Growth Mixture Model implemented in Stan

[W Robert Long]

Hi all

I would like to know if anyone has implemented a Growth Mixture Model
(GMM) in Stan ? More generally I would like to know the experts opinions
on the suitability of Stan for estimating GMMs

To avoid confusion, here I consider a GMM to be an extension of a Latent
Growth Curve model (LGCM) in the following sense. In LGCM, repeated
measures of some variable are modelled as a function of a small number
of latent variables, for example an intercept and slope. Such a model
can also be considered as a mixed effects / multilevel model where the
random effects (random intercepts and random slopes in this case)
correspond to the latent growth variables.

GMM can then be viewed as an extension of LGCM which allows for
differences in latent growth variables across unobserved subpopulations,
modelled as latent classes.

Individuals are classified by estimating their posterior probabilities
of class membership. The number of latent classes is generally specified
in advance - models with different numbers of classes are compared using
statistics such as BIC (as well as expert a priori considerations) to
determine the optimum number of classes.

Covariates may enter the model as predictors of the latent growth
variables, and/or as predictors of the latent class variable.

In more complex models, the latent class variable may also predict a
distal outcome variable.

The variance-covariance matrix of the latent growth variables is often
an area of modelling focus.

The variance-covariance matrix of the measurement-level residual
variances is also an area of modelling focus - for example they are
sometimes all constrained to be equal, sometimes freely estimated, and
sometimes assumed to be autoregressive.

Convergence is often a problem with these kinds of finite mixture
models. The variance of the latent growth variables is sometimes
constrained to zero in order to achieve convergence. MPlus is often
used, using maximum likelihood (EM algorithm) with thousands of
different random starting values to overcome the problem of local
solutions.

I wondered whether a fully Bayesian approach in Stan with
suitable/flexible choices of priors would be a good alternative ?

In a nutshell - is Stan a good choice for fitting these kinds of models
? Has anyone has any success ? If so I would love to hear about the
results, particularly compared to using Mplus.

Thanks
RL

[Ben Goodrich]

I don't know of anyone doing a model like this, but I can't see any reason off the top of my head why it would be a problem. It would be even better after the next Stan release when we have lower-triangular matrices as a variable type.

Ben

[Bob Carpenter]

Cholesky factors for correlation matrixes may have to wait.
We're getting bogged down in ODE solvers and getting the
basic 2.0 out the door at the moment.

But it's the next thing I plan to add to the language when
I free up some time.

- Bob

> --
> 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/groups/opt_out.

[Bob Carpenter]

You probably want to take a look at the mixture modeling
chapter in the Stan manual -- it's a bit tricky how to do it.
No problem in theory as long as there are preset numbers of
classes (or you can approximate unbounded with a big fixed number).

Stan's set up to let you change priors easily, so moving
from independent to time series is no big deal. It can have
a big impact on sampling, so parameterization and efficiency
become important considerations.

As to particular priors, that's going to depend on your data
and how much of it there is. The less data there is, the
more careful you have to be with priors to control computation
(and control the posterior for inference).

- Bob

[Sierra]

(I am happy to be included in reference materials, etc., no restrictions)

Hello Stan Users,

I am also trying to estimate a growth mixture model in Stan, but I am not sure if I am doing it correctly. I have read the mixture model chapter from the manual. If I am doing it correctly, than I am almost surely doing it poorly. So far the model compiles and begins sampling, but on my machine it has taken about 6 hours to get to 10% of 2000 iterations [Warmup] for Chain 1. In a previous attempt to let the model run overnight, just a single chain, 200 iterations, the effective sample size for k(# of components) was very small relative to the other parameters, but I don't remember the specifics before my computer froze.

The model I am trying to estimate is just as Robert described. In this case the repeated measures are described as a function of three potentially correlated latent variables, vector "eta" ("eta[,1]"= intercept, "eta[,2]"= slope, "eta[,3]"= quadratic). There is only a single estimated error variance "sigma_y", so the residual variance is equal at each time point. The only parameters that vary by class are the means on the latent variables "mu_eta". The parameters of interest are the class-specific mean vectors for eta, the covariance matrix among the etas ("Psi"), the probabilities of belonging to each class "pi", and the residual error variance. 

I am attaching example fake data for a two-class quadratic growth mixture model, with 5 repeated observations "T" for 2000 individuals "N" (the data has 6 columns, the last column is the true latent class assignment). I was able to get this model to work as a latent growth curve model in the log-scale (without bringing in the classes), but in translating to a mixture model it seems to be going awry. 

Any feedback is appreciated. Thank you.

My code setup (in R):

y<-read.table('gmmq1.dat')

y<-y[,1:5]

x<- c(0, 1, 2, 3, 4)

q<- c(0, 1, 4, 9, 16)

Omega <- structure(c(5, 0, 0, 0, 5, 0, 0, 0, 5), .Dim = c(3, 3)) 

gmm_dat<- list(N = nrow(y), T = ncol(y), x = x, q = q, y = y, Omega=Omega, K=2)

gmm_fit <- stan(file = 'gmm-q.stan', data = gmm_dat, verbose = FALSE)

Stan code:

data {

  int<lower=0>  N;// number of ppl

  int<lower=0>  T;// number of timepoints

  real x[T];      // X (time)

  real q[T];      // q (time^2)

  real y[N,T];    // Y matrix of T obs of N ppl

  int <lower=1> K;// number of classes

  cov_matrix[3] Omega;

}

transformed data {

}

parameters {

  simplex[K] pi;

  vector[3] eta[N];  // eta (int, slope, quad) for each person

  vector[3] mu_eta[K];  // mean (int, slope, quad) for etas

  cov_matrix[3] Psi; // cov matrix among etas

  real<lower=0> sigmasq_y; // y residual error variance 

}

transformed parameters {

  real<lower=0> sigma_y;

  vector[K] log_q_z[N];

  sigma_y <- sqrt(sigmasq_y);

  for (n in 1:N){

    log_q_z[n]<-log(pi);

    for (t in 1:T)

      for(k in 1:K)

         log_q_z[n,k]<-log_q_z[n,k]+normal_log(y[n,t], eta[n, 1] +eta[n, 2]*x[t] +eta[n,3]*q[t], sigma_y);

  }

} 

model {

   real tt[K];

   mu_eta[K] ~ normal(0, 100);

   sigmasq_y ~ inv_gamma(0.001, 0.001); //y residual error variance drawn from prior

   Psi~ inv_wishart(3, Omega); //cov matrix among etas drawn from prior  

   for (n in 1:N)

     for (k in 1:K)

      tt[k]<-log(pi[k])+multi_normal_log(eta[n],mu_eta[k],Psi);

   for (n in 1:N){ 

       increment_log_prob(log_sum_exp(log_q_z[n]));

       increment_log_prob(log_sum_exp(tt));

  } 

}

[Bob Carpenter]

I can make several suggestions on model and implementation.

First, don't use any of these priors:

mu_eta[K] ~ normal(0, 100);
sigmasq_y ~ inv_gamma(0.001, 0.001); //y residual error variance drawn from prior
Psi~ inv_wishart(3, Omega); //cov matrix among etas drawn from prior

They're all problematic for both computational and statistical
reasons. The fat normal and inverse gamma are just too fat.
the inverse Wishart usually doesn't make any sense. These priors
just get used by BUGS/JAGS becaus they're conjugate. We discuss
recommended alternatives in the manual chapter on regression.

We recommend a Cholesky-parameterized, scaled correlation matrix with
an LKJ prior for Psi, and at least weakly informative priors for the
others. There are references as to why in the manual, too.

Computationally, you probably don't want to define all those variables
as transformed parameters, because then they get written out at each
iteration.

I'm not sure what you're trying to do, but this doesn't make sense:

> for (n in 1:N)
> for (k in 1:K)
> tt[k]<-log(pi[k])+multi_normal_log(eta[n],mu_eta[k],Psi);

It redefines tt N times, then uses the last one (n == N) going forward.

I'm lost as to what this one's doing:

> for (n in 1:N){
> log_q_z[n]<-log(pi);
> for (t in 1:T)
> for(k in 1:K)

> log_q_z[n,k] <- log_q_z[n,k]+normal_log(y[n,t], eta[n, 1] +eta[n, 2]*x[t] +eta[n,3]*q[t], sigma_y);

The variables are:

// x : T x 1
// pi : K x 1
// eta: : N x 3
// q : T x 1
// sigma_y : 1

so there's nothing in that loop other that varies by k other than log_q_z.
If that's really what you want (I don't really have time to figure it out),
then you're far better of using vectorization on the normal_log over t
outside of the loop over K --- as is, you're doing way too much extra
work.

So that'd be something like this:

{
vector[K] log_pi;
log_pi <- log(pi);

for (n in 1:N) {

real lp;
lp <- normal_log(y[n],
eta[n,1] + eta[n,2] * x + eta[n,3] * q,
sigma_y);

for (k in 1:K)

log_q_z[n,k] <- log_pi[k] + lp;
}
}

But I'm not convinced you're doing what you want to do here.

- Bob

> --
> 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.
> <gmmq1.dat>

[Sierra]

Thank you Bob.

Okay, I read the regression chapter and I think I understand why to make the changes to the priors and parameterization (though I haven't Cholesky-factored yet), thank you.

As you can see from my original code, I am quite confused about summing out over classes K while incrementing the log probability, etc. A part of my confusion is what to sum over because the model is hierarchical but not everything varies over classes... and I am new around here.

Only the vector on the means varies over classes.... is that the only portion where I need to sum over classes? The core model is simply something like:

    eta[n]~MVN[3](mu_eta[k],Sigma_eta)

    y[n]~N(X*eta,sigma_y)

Right now my code isn't working at all (ill-formed expression for the multiplication in the normal log.... I have tried this umpteen ways), but something like:

data {
  int<lower=0>  N;// number of ppl
  int<lower=0>  T;// number of timepoints

  real x[T];      // X (time loadings)
  real q[T];

  real y[N,T];    // Y matrix of T obs of N ppl
  int <lower=1> K;// number of classes
}

parameters {
  simplex[K] pi;           // class weights
  corr_matrix[3] Psi;      // prior correlation
  vector<lower=0> [3] tau; // prior scale

  vector[3] mu_eta[K];     // mean (int,slope,quad) for etas

  vector[3] eta[N];        // etas (int,slope,quad) for each person
  real<lower=0> sigmasq_y; // y error variance

}
transformed parameters {
  real<lower=0> sigma_y;

 
  sigma_y<-sqrt(sigmasq_y);
}
model {
  matrix[3,3] Sigma_eta;
  real ps[K];
  Sigma_eta<-quad_form_diag(Psi,tau);
  mu_eta[K] ~ normal(0, 10);
  tau~cauchy(0,2.5);
  Psi~lkj_corr(2);

  for (n in 1:N){

    for (k in 1:K){
      ps[k]<-log(pi[k])+ multi_normal_log(eta[n],mu_eta[k],Sigma_eta);
    }
    increment_log_prob(normal_log(y[n], eta[n,1]+eta[n,2]*x+eta[n,3]*q,sigma_y));
    increment_log_prob(log_sum_exp(ps));
  }
}

Is this closer to summing over the classes and incrementing the log probability correctly? Thank you again.

[Bob Carpenter]

I'm not sure what your model is. The first thing to do
is write the model down in math. Then if it has discrete
parameters, marginalize them out.

None of that depends on Stan.

When you have that, we can help you translate to Stan.
As is, I don't understand your model well enough to help
and trying to infer it from a broken Stan program isn't
so easy.

I can tell you that the problem in your multiplication is
that you can only multiply vectors, not arrays. There's a
chapter in the manual on arrays vs. vectors you might find
useful in this regard.

- Bob

[Sierra]

Hi Bob,

Sorry for the radio silence, I drafted a reply but apparently never sent it. 

This is the model in math, I hope this is clear. 

My updated code attempt to sum over the classes compiles without errors, runs for a day or so, but comes nowhere near convergence... 

(I don't remember the details of non-convergence unfortunately, they were included in the lost draft from almost a month ago.)

r setup:

#example with GMM data

y<-read.table('gmmq1.dat')
y<-y[,1:5]


x<- c(0, 1, 2, 3, 4)
q<- c(0, 1, 4, 9, 16)

lcm_dat<- list(N = nrow(y), T = ncol(y),
               x = x, q=q, y = y, K=2)
gmm_fit <- stan(file = 'gmm-q.stan', data = lcm_dat, verbose = FALSE)

gmm-q.Stan file:

# simple GMM

data {
  int<lower=0>  N;// number of ppl
  int<lower=0>  T;// number of timepoints

  vector[T] x;      // X (time loadings)
  vector[T] q;
  matrix[N,T] y;    // Y matrix of T obs of N ppl

  int <lower=1> K;// number of classes
}

transformed data {
}
parameters {

[Bob Carpenter]

It might help to put a prior on sigma (you can do it directly
or by putting a prior on sigma^2 and transforming). You'd have
to see if too large sigma was causing problems.

Otherwise, you have to look and see what's correlated or not
fitting well and try to adjust priors, etc.

There are two big optimizations for what you're doing.

First, you want to use a Cholesky-factor parameterization of
the prior --- there's a description in the manual chapter on
regression in a section called something like "multivariate priors".

Second, it would be much more efficient to pull this out of the loop

> increment_log_prob(normal_log(y[n], eta[n,1]+eta[n,2]*x+eta[n,3]*q,sigma_y));

and write this outside the loop:

y ~ normal(col(eta,1) + col(eta,2) * x + col(eta,3) * q, sigma_y);

- Bob

> <GMM.pdf><gmmq1.dat>