How would you parametrize a markov chain model and its hierarchical extension?

[matus.s...@privatdemail.net]

I have a Markov chain model with four states A, B, C and D. The transition probabilities are given by the 4x4 transition matrix M. The data are given as chains of states e.g. AAABBBACCCAA .... How would you parametrize this model?  How would you parametrize a hierarchical extension of the model? For hierarchical extension assume that I have N subjects and for each subject I have M_n chains.

Idea 1: Estimate transition probabilities directly. I would then use inverse logit to transform the transition probabilities and pool the transformed vars with a hierarchical gaussian prior.

Idea 2: In the actual data, the states repeat very often. Same-state chains of length 100 and more are not unusual. As a consequence the values on the diagonal of the transition matrix will be near 1 and other elements will be near 0. My alternative idea is to estimate the transition probabilities between different consecutive states. I would estimate the length of a chain of same-state transitions with a poisson model. Again, there is a multilevel generalization of poisson model in chap. 15 of gelman's DAUR&MHM. I would adapt this.

I haven't coded anything yet. I'm just interested in your opinion. Which idea is better or do you see alternative formulations?

My hunch is that idea 2 is easier to estimate with STAN and also the parameters are easier to interpret. However, I've never seen idea 2 in the literature, so I'm not very confident.

[Dan Stowell]

Hi -

Idea 2 is called the "semi-Markov model". I hope that helps your
literature search :)
Egotistically enough, I'll point you to a short note I wrote recently
about what to call them:
http://rspb.royalsocietypublishing.org/content/281/1792/20141370.e-letters

In response to the question about choosing between the MM and SMM, I'd
say choose whichever fits closest to your beliefs about the underlying
behaviour. Note that a non-hidden Markov model will always generate
"holding times" of geometrically-distributed length. If a priori you
don't think that fits well to your situation, then that's a good
argument to use your SMM.

Best
Dan

> --
> You received this message because you are subscribed to the Google Groups
> "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.



--
http://www.mcld.co.uk

[matus.s...@privatdemail.net]

Hi Dan, thanks, that's very helpful to know what the model is called.

I'm still looking for some feedback on hierarchical pooling of the transition probabilities. Maybe, someone has experience with that?

m.

[Luc Coffeng]

If you define transition matrix M per individual (M_i) and use M_i for all chains of states for a given individual, you can use a dirichlet hyperprior for columns of transition probabilities in your transition matrix M_i. However, this would only pool information for each column separately, and not account for possible correlation between column within a given matrix M_i. I don't know whether there are priors for whole transition matrices - perhaps this is something you can find literature.

[Bob Carpenter]

The simplest thing to do would be what Luc suggested --- use a Dirichlet
prior and give it some kind of hyperprior (typically reparameterizing
in terms of a simplex mean alpha/sum(alpha) and a prior on the total
count sum(alpha). Then each state would draw its transition probabilities
from Dirichlet(alpha).

You could enrich this using a logistic multinormal. Just use
a multivariate normal, then softmax transform the draws to get
your simplexes for the transitions. This was one of the first
pieces of stats advice Andrew ever gave me, long before we
started working together.

- Bob

[Luc Coffeng]

Thanks Bob for clarifying what I meant! And of course, the trusted multivariate normal distribution is *the* way to capture the correlation between the columns. Care should be taken though to restrict the support of the alphas to positive numbers. For identifiability, fix the mean of the multinormal hyperprior to 0.0, fix sd to 1.0, and only estimate the cholesky decomposition of the correlation matrix.

The pseudo-code for a N*N transition matrix would then be:
parameters:
corr_cholesky[N] L;
vector<lower=0>[N] alpha_sum;  // column-specific sample sizes for Dirichlet prior for transition probabilities
vector<lower=0>[N*N] alpha_u;  // normalized untransformed alpha

model block:
vector[N*N] alpha_raw;
matrix[N,N] alpha; // each column is a vector of alpha parameters of the dirichlet prior for a column

alpha_u ~ normal(0, 1);
L ~ lkj_prior(...);  // if you want something else than a uniform prior

alpha_raw <- u %*% L;  // alpha_raw ~ multinormal(0, LL')T[0,]

for (n in 1:N) {  // not sure the indexing is right here, check that values are organized into columns
  alpha[n] <- to.column(softmax(n-th set of N elements from alpha_raw) .* alpha_sum[n]);
}

[Luc Coffeng]

As I'm riding my bike home, I realise my previous post was (partly) hogwash.

For identifiability, only fix the sd of the multinormal hyperprior for alpha_u; the mean vector is a parameter to be estimated. Further, for identifiability, you may fix the first element of each of the N sets of length N to the value 1.

Unfortunately, I don't see a way to implement the mu + u*L trick here to simulate the multinormal random alpha_raw, as the results of this trick needs to be positive (for the softmax). Will a call to multinormal_cholesky be very ineffecient?

The correct pseudo code then becomes (somebody please correct me if I'm wrong):

parameters:
corr_cholesky[N] L;
vector<lower=0>[N] alpha_sum;  // column-specific sample sizes for Dirichlet prior for transition probabilities

vector<lower=0>[N*N] alpha_u_mu;   // mean untransformed alpha (each column of the transition matrix has it's own N-vector)
vector<lower=0>[N*N)] alpha_raw;  // untransformed alpha

model block:
vector[N*(N-1)] alpha_raw;

matrix[N,N] alpha; // each column is a vector of alpha parameters of the dirichlet prior for a column

alpha_u ~ normal(0, 1);
L ~ lkj_prior(...);  // if you want something else than a uniform prior

alpha_raw ~ multinormal_cholesky(alpha_u_mu, L);  // I don't see a way to comobine the support of alpha_raw
                                                  // with the mu + u*L trick.

[Luc Coffeng]

Duh, correction:

On Thursday, 5 March 2015 17:23:51 UTC+1, Luc Coffeng wrote:

As I'm riding my bike home, I realise my previous post was (partly) hogwash.

For identifiability, only fix the sd of the multinormal hyperprior for alpha_u; the mean vector is a parameter to be estimated. Further, for identifiability, you may fix the first element of each of the N sets of length N to the value 1. <<-- ignore last sentence (grey); fixing the sd should provide soft identifiability.

[Bob Carpenter]

More precisely, here's what I was thinking:

vector[M] theta_raw[M];
for (m in 1:M)
theta_raw[m] ~ multi_normal...

simplex[M] theta[M];
for (m i 1:M)
theta[m] <- softmax(theta_raw[m]);

That lets you have a prior with covariance on theta_raw, then
it gets transformed with softmax to a proper transition matrix.

As you bring up, there's the issue of efficiently parameterizing
that multi_normal, and there's also the issue of identifiability
using a size-M vector.

What we really need to do is add versions of softmax that
involve only (M-1) parameters and treat the Mth as 0. That
provides the necessary identification, but then it makes
the prior hard to interpret. We've discussed this before among
me, Michael, and Andrew, but without much consensus. The
(M-1) version is better identified, so should sample better,
the M version is better in terms of efficiency.

For now, to get the (M-1) version, you could instead take

vector[M-1] theta_raw_raw[M];
for (m in 1:M)
theta_raw_raw ~ multi_normal ...

vector[M] theta_raw[M];
for (m in 1:M) {
for (n in 1:(M-1))
theta_raw[m,n] <- theta_raw_raw[m,m];
theta_raw[m,M] <- 0;
}

simplex[M] theta[M];
for (m i 1:M)
theta[m] <- softmax(theta_raw[m]);

Ugly, but it gives you the K-1 parameterization.

- Bob

[Luc Coffeng]

Thanks for clarifying! Regarding the following piece:

vector[M] theta_raw[M];
for (m in 1:M)
  theta_raw[m] ~ multi_normal...

This assumes that for a given individual, the columns of the final transition matrix M are i.i.d. following the multinormal distribution. Not sure whether this is valid for this particular application. It might well be that the transition probability M[1,2] is correlated with transition probability M[2,3]. To capture that, one would have to define theta_raw as a vector of length m^2 and evaluate that with a M^2-dimensional multinormal distribution, and then do the softmax transform on each of the m sets of m-length theta_raws.

(Theta_raw needs to be constrained to be positive.)

[Bob Carpenter]

> On Mar 6, 2015, at 7:06 PM, Luc Coffeng <lucco...@gmail.com> wrote:
>
> Thanks for clarifying! Regarding the following piece:
>
> vector[M] theta_raw[M];
> for (m in 1:M)
> theta_raw[m] ~ multi_normal...
>
> This assumes that for a given individual, the columns

I always get confused because the choice of row vs. column
is arbitrary.

What matters is that for each state m there is a simplex
of transition probabilities to each possible next state 1,...,M.

What I was suggesting is that you model the states hierarchically
with a multivariate logistic prior. That allows for modeling
covariance among the transitions. For instance, transitions to m1
may be correlated with transitions to m2.

> of the final transition matrix M are i.i.d. following the multinormal distribution.

Not quite --- it's multivariate logistic. You could
also do multivariate probit by using the probit link
function (cumulative normal Phi).

> Not sure whether this is valid for this particular application. It might well be that the transition probability M[1,2] is correlated with transition probability M[2,3].

Without replicates, there'll be no way to estimate that. What
I was suggesting has M replicates, one for each state.

Of course, if you have bunches of problems with the same HMM
state structure then you could do it.

Then I think you could generalize to matrix normal distributions,
which I believe Rob's about to add to Stan:

http://en.wikipedia.org/wiki/Matrix_normal_distribution

> To capture that, one would have to define theta_raw as a vector of length m^2 and evaluate that with a M^2-dimensional multinormal distribution, and then do the softmax transform on each of the m sets of m-length theta_raws.
>
> (Theta_raw needs to be constrained to be positive.)

The point was that theta_raw is unconstrained and each
set of outgoing transition probabilities theta_raw[m] can
be normalized to a simplex softmax(theta_raw[m]). That induces
a hierarchical structure on the outgoing transitions that allows
correlation between the targets. You could do the same kind of
thing with the Dirichlet, which gives you mean and variance, but
not covariance.

- Bob

[Luc Coffeng]

Thanks for clarifying!

And yes, I do believe the original poster mentioned that the data included multiple state chains per individual.

[matus.s...@privatdemail.net]

Bob, Luc,

thanks for help. My original goal was to pool transition probabilities across subjects. I haven't thoght about correlation between the columns of the transition matrix (each column=outgoing state). But once I have a simple model running, I may look into that.

[Bob Carpenter]

Oh, then you can do the same thing, but with the
whole matrix. You'll need to flatten it into a vector
if you want to have covariance. Or you can just
pool per transition without allowing covariance across
states. Or you can use a Dirichlet for each transition
probability simplex if you don't care about covariance.

- Bob

> On Mar 7, 2015, at 3:52 AM, matus.s...@privatdemail.net wrote:
>
> Bob, Luc,
>
> thanks for help. My original goal was to pool transition probabilities across subjects. I haven't thoght about correlation between the columns of the transition matrix (each column=outgoing state). But once I have a simple model running, I may look into that.
>
>