Multivariate ordered logistic model not converging

[Kevin Rush]

I am trying to model ordinal data from a survey as a multivariate ordered logistic model. I am finding that as I increase the number of questions in the data (from 10 items to 16 items), the chains do not converge. The output for all of the chains in the 10 item dataset show convergence, but for the 16 item dataset, only 1 of the chains converge. 

The underlying latent parameter for the questions are modelled as a multivariate normal distribution, and each integer response as an ordered logistic. I am modelling the prior for the cutpoints for each item as a fixed-scaled cumulative sum of a symmetric dirichlet. For identifiability, I set the first cutpoint to 0. I am running on CmdStan (v2.4), Linux Ubuntu.

Do you have any suggestions for choices of priors or scales that might help convergence? My eventual goal is to include all the data, which is 56 questions.

[Xulong Wang]

Hi Kevin,

I am also modeling a ordered logistic with Stan. What brought you to set the cut points' prior as fixed-scaled cumulative sum of a dirichlet? How different are the results comparing with a non-informative prior on cut points by not specifying anything?

My model might also has the identificability issue given that intercept is included. The cut points and intercepts might run into each other. But I am not sure because I am modeling the intercept as a random effect. As shown in my model below, do I also need to specify priors for muAlpha and sigmaAlpha?

---
code <- '  // STAN: hierarchical ordered logistic model
data {
  int<lower=0> N;  // number of objects
  int<lower=1> K;  // number of categories
  int<lower=1> L;  // number of groups
  int<lower=1,upper=K> y[N];  // response
  vector<lower=0,upper=1>[N] sex;  // predictor
  vector<lower=1>[N] age;  // predictor
  int<lower=1,upper=L> group[N];  // group index
}
parameters {
  vector[L] alpha;  // intercept
  real muAlpha;  // mu for intercept
  real<lower=0> sigmaAlpha;  // sigma for intercept
  real beta1;  // slope for sex
  real beta2;  // slope for age
  ordered[K-1] c;  // cutpoints for ordered logistic model
}
transformed parameters {
  vector[N] y_hat;
  for (n in 1:N)
    y_hat[n] <- alpha[group[n]] + beta1 * sex[n] + beta2 * age[n];
}
model {
  alpha ~ normal(muAlpha, sigmaAlpha);  // intercept
  /* priors for muAlpha, sigmaAlpha, c can be further specified */

  for (n in 1:N)
    y[n] ~ ordered_logistic(y_hat[n], c);
}
generated quantities {
  vector[N] loglik;
  for (n in 1:N)
    loglik[n] <- ordered_logistic_log(y[n], y_hat[n], c);
}
'
---

[Kevin Rush]

I got the idea for setting the cut point prior as fixed-scaled cumulative sum of a dirichlet from this thread:

https://groups.google.com/d/msg/stan-users/2YMo3SUzHTg/jgeoyZfUSBsJ

When the cutpoints have a non-informative prior, I found it was extremely slow to converge (over 10 times as long) versus setting a more informative prior. However, looking at the older thread, the poster seems to have convergence issues as well in some of the chains with this approach.

I'm definitely not an expert, but I have found that constraining the prior helps with non-identifiability issues and convergence speed. You might want to put a broad prior, even a normal(0,10) or cauchy(0,10) (for longer tails) on muAlpha and sigmaAlpha.

[Xulong Wang]

Thanks Kevin for your kind response. Put priors could possibly affect the posterior in a bad way if the priors are not appropriate. I agree with you it is helpful to put a broad prior to speed up the convergence process. Also, put distinct priors to parameters can solve the non-identifiability issue if I am understanding the sampling process correctly. I will give all of them broad norm in my model.

Somehow, I feel there should be a general guidance on putting priors to the common models, given this is such important and obvious question. Let me search the literature more.

[Bob Carpenter]

We don't (yet) have any default recommendations for priors on
cutpoints.

The main issue you have to deal with is making sure you
don't have multiplicative (scale) or additive (location)
identifiability issues with the cutpoints and a regression.

If you give the regression a mean and the cutpoints float,
then you have a classic non-identifiability problem. You
can pin the value of one of the cutpoints to identify the
location. Then if you still need to worry about identifying
the scales.

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

[Xulong Wang]

Thanks Bob, I am trying to fix the cut points' first element as 0 by first declaring an ordered vector, then to substrate the first element.

However, the compilation does not go through with the message: variable "vector" does not exist. Please see the code below. Something wrong happens in the transformed parameter block.

Then, for a categorical predictor, how can I manually fix the coefficient for a given group as 0?

code <- '  // STAN: hierarchical ordered logistic model
data {
  int<lower=0> N;  // number of objects
  int<lower=1> K;  // number of categories
  int<lower=1> L;  // number of groups
  int<lower=1,upper=K> y[N];  // response

  matrix[N,2] x;  // predictor matrix

  int<lower=1,upper=L> group[N];  // group index
}

parameters {
  vector[L] alpha;  // intercept
  real muAlpha;  // mu for intercept
  real<lower=0> sigmaAlpha;  // sigma for intercept

  vector[2] beta;  // predictor coefficients
  ordered[K-1] cu;  // cut points
}

transformed parameters {
  vector[K-1] c;  // cut points
  for (k in 1:(K-1))
    c[k] <- cu[k] - cu[1];

  vector[N] y_hat;
  for (n in 1:N)

    y_hat[n] <- alpha[group[n]] + x[n] * beta;
}

model {
  muAlpha ~ normal(0, 5);
  beta ~ normal(0, 5);

  alpha ~ normal(muAlpha, sigmaAlpha);  // intercept

  for (n in 1:N)
    y[n] ~ ordered_logistic(y_hat[n], c);
}

generated quantities {
  vector[N] loglik;
  for (n in 1:N)
    loglik[n] <- ordered_logistic_log(y[n], y_hat[n], c);
}
'

[Daniel Lee]

For future reference, it helps to paste the whole error message. The error message provides context and makes it much easier to find the error.

In Stan, you have to declare all variables at the top of its scope. In this case, inside transformed parameters, y_hat needs to be moved up a few lines.

[Xulong Wang]

Thanks Daniel, the model is now somehow complete to me, and the computation time is acceptable.

Bob, I have now fixed the first cut point as 0 for the identificability issue with the intercept. The factorial predictor issue was handled by simply ignoring the reference group. Could you please take a look at the model below before my serious run. The model in Stan is also attached if it is more convenient to you.

---
data {
  int<lower=1> N;  // Sample number
  int<lower=1> K;  // AD categories
  int<lower=1> L;  // Family number
  int<lower=1,upper=K> Ad[N];  // AD status
  int<lower=1,upper=L> Family[N];  // Family index
  int<lower=1,upper=2> Sex[N];  // M:F
  int<lower=1,upper=5> Apoe[N];  // 33:22:23:24:34
  vector[N] Age;
}

transformed data {
  vector[N] sexF;  // Relative to M
  vector[N] apoe22;  // Relative to 33
  vector[N] apoe23;
  vector[N] apoe24;
  vector[N] apoe34;

  for (i in 1:N) {
    sexF[i]  <- Sex[i]  == 2;
    apoe22[i] <- Apoe[i] == 2;
    apoe23[i] <- Apoe[i] == 3;
    apoe24[i] <- Apoe[i] == 4;
    apoe34[i] <- Apoe[i] == 5;

  }
}

parameters {
  vector[L] alpha;  // intercept
  real muAlpha;  // mu for intercept
  real<lower=0> sigmaAlpha;  // sigma for intercept

  real betaAge;  // Age
  real betaSex;  // Female
  real betaApoe22;  // Relative to 33
  real betaApoe23;
  real betaApoe24;
  real betaApoe34;

  ordered[K-1] cu;  // cut points
}

transformed parameters {
  vector[K-1] c;  // cut points

  vector[N] yHat;

  for (k in 1:(K-1))
    c[k] <- cu[k] - cu[1];

  for (n in 1:N)

    yHat[n] <- alpha[Family[n]] + Age[n] * betaAge +
               Sex[n] * betaSex +
               apoe22[n] * betaApoe22 +
               apoe23[n] * betaApoe23 +
               apoe24[n] * betaApoe24 +
               apoe34[n] * betaApoe34;

}

model {
  muAlpha ~ normal(0, 5);

  betaAge ~ normal(0, 5);
  betaSex ~ normal(0, 5);
  betaApoe22 ~ normal(0, 5);
  betaApoe23 ~ normal(0, 5);
  betaApoe24 ~ normal(0, 5);
  betaApoe34 ~ normal(0, 5);

  alpha ~ normal(muAlpha, sigmaAlpha);  // intercept

  for (n in 1:N)

    Ad[n] ~ ordered_logistic(yHat[n], c);
}

generated quantities {
  vector[N] logLik;

  for (n in 1:N)

    logLik[n] <- ordered_logistic_log(Ad[n], yHat[n], c);

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

Three suggestions.

1. You don't need to save yHat as a transformed parameter if you don't
need it in the output. You can make it a local variable in the model.
Then you'd have to recompute it in the generated quantities block, but
you can also do that as a local variable. Then you don't output yHat.

2. You don't want to try to set c[k] to cu[k] - cu[1]. That doesn't
identify the cu because they can all shift and you get the same behavior.
Instead, do something like declaring cu as pos_ordered[K-2], then setting
c[1] to 0, and all c[k] to cu[k-1].

3. You can use the non-centered parameterization for alpha. Define
alpha_raw ~ normal(0,1) and then define alpha <- muAlpha + sigma_alpha * alpha_raw.
This is much faster when there's not much data to estimate alpha with.
It may not matter here if alpha is well identified by the model.

- Bob


On Oct 18, 2014, at 3:48 AM, Xulong Wang <emai...@gmail.com> wrote:

> ...
> parameters {
> ...

> <orderedLogitV2.stan>

[Xulong Wang]

Hi Bob, thanks for all the suggestions. All well taken. I will keep you updated how the project goes with Stan.

Best,

Xulong

[Xulong Wang]

Bob, 

A related question with your first suggestion. Local variables in the generated quantities block are actually all saved in the output, at least in my case.  

However, seems there is an argument in the stan function, par, which will declare the parameters you want to save. Want to confirm with you the "par" argument in the stan or sampling function does not affect the sampling process, but only select parameters to be saved. Am I correct?

X.

[Daniel Lee]

Hi X,

On Wednesday, November 5, 2014, Xulong Wang <emai...@gmail.com> wrote:

A related question with your first suggestion. Local variables in the generated quantities block are actually all saved in the output, at least in my case.  

None of your generated quantities are within a local block. Try putting something like this inside your generated quantities block, including the addition curly braces:

{

   real x;

   x <- 1;

}

 

x should be saved to output.

However, seems there is an argument in the stan function, par, which will declare the parameters you want to save. Want to confirm with you the "par" argument in the stan or sampling function does not affect the sampling process, but only select parameters to be saved. Am I correct?

This is for RStan. We'll let some of the RStan users let you know what the behavior is.

[Bob Carpenter]

The variables with scope over the blocks are saved. That
means in the transformed data, transformed parameters, and generated
quantities block, the initial sequence of variable declarations are
the ones that are saved.

For ones you don't want to save, you can make a local block
using curly braces. For example, in

generated quantities {
real gq_scope;
{
real local_scope;
... use either var ...
}
... only use gq_scope ...
}

the variable "gq_scope" scopes over the entire generated quantities
block, whereas the variable local_scope declared inside curly braces
can only be used inside those braces and won't be saved.

- Bob

[Bob Carpenter]

Sorry for the duplication --- I hadn't seen this. One correction
below.

> On Nov 5, 2014, at 8:40 AM, Daniel Lee <bea...@alum.mit.edu> wrote:
>
> None of your generated quantities are within a local block. Try putting something like this inside your generated quantities block, including the addition curly braces:
> {
> real x;
> x <- 1;
> }
>
> x should be saved to output.

That should be:

x WILL NOT be saved to output.

- Bob

[Xulong Wang]

Thanks Bob, this is clear.