Markov chains, hysteresis and logistic regression.

Logistic regression models the response of a binary variable to a covariate. The binary data produced by a two state Markov process coupled to an external covariate can be fit using logistic regression. Under what circumstances, if any, is this a valid approach?

At each time step $i$ , the non-stationary Markov process is governed by a transition matrix

(1) $\begin{equation*}P_i=\left(\begin{array}{cc} p_i & 1-p_i \\ 1-q_i & q _i\end{array}\right)\end{equation*}$

where $0 \leq p_i,q_i \leq 1$ . $P$ can be coupled to the external covariate $x$ using logistic functions, which preserve the stochastic matrix property. For instance, choose:

(2) $\begin{eqnarray*} p_i &=& \frac{1}{1+e^{2\left(1+2x_i\right)}}\\ q_i &=& \frac{1}{1+e^{2\left(1-2x_i\right)}} \end{eqnarray*}$

Also, assume that $x$ varies harmonically with timescale (angular period) $T$ i.e.

(3) $\begin{equation*} x_i= \sin \left( \frac{i}{T} \right) \end{equation*}$

The graph below the expected response (a number between 1 & 2) as a function of $x$ obtained by simulating 200 Markov chains of length $10^5$ using Equations (1)-(3).

The striking feature is the appearance of an hysteresis loop. Hysteresis, where the mean response has separate branches for increasing and decreasing $x$ , becomes pronounced when $x$ varies rapidly (small $T$ ). On the other hand, when $x$ varies slowly, the mean response collapses to a single branch.

What’s going on?

Starting from some initial state, a markov chain approaches it’s equilibrium distribution after a mixing time. The mixing time is longer when there are bottlenecks e.g.when $q \simeq 1$ in (1). When the covariate varies on a time-scale which is long compared to the mixing time, the system remains close to equilibrium at all times. Hysteresis is absent in this case. Conversely, when the covariate varies rapidly the transition matrix changes too quickly for the system to reach equilibrium. This leads to hysteresis.

Logistic regression is valid as long as the mixing time is short compared to the timescale over which the covariate varies. However, if the mixing time is long, multiple branches of the response are present. This aspect of Markov chain dynamics is missing from ordinary logistic regression. This conclusion generalises to the N-state/multinomial regression case.

Markov chain simulation

A good way simulate non-stationary Markov chains in R is to pass a list of transition matrices to a C++ function (MCsim1) using Rcpp. For example:

Rcpp::sourceCpp('~/markov.cpp')
x <- sin(1:100000/4)
T.ts <- lapply(x,function(d) Tmat(1/(1+exp(-2*(1-2*d))),1/(1+exp(-2*(1+2*d)))) )
simulation <- MCsim1(T.ts,start=1)

where “~/markov.cpp” contains:

// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadilloExtensions/sample.h>

using namespace Rcpp ;

// [[Rcpp::export]]

IntegerVector MCsim1(Rcpp::List t,int start) {

  NumericMatrix t1 = t[1];
  int nrow = t1.nrow();
  int n = t.length();
  NumericMatrix s;
  IntegerVector states = seq_len(nrow)-1;
  IntegerVector sim(n);
  sim(0) = start;
  for (int j = 1; j < n; j++) {
    NumericMatrix s=t[j];
    NumericVector probs = s(sim(j-1), _);
    IntegerVector newstate = RcppArmadillo::sample(states, 1, false, probs);
    sim(j) = newstate(0);
  }
  for (int j = 0; j < n; j++) sim(j) = sim(j)+1;
  return sim;
}

What’s going on?

Markov chain simulation

Leave a Comment Cancel reply