Help with hierarchical linear model

[David Westergaard]

Hi everyone,

I am trying to develop a model with four nominal variables (age, gender, analysis method, disease type), and the interaction between some of these. The setup is a kind of repeated measures. Some of the participants have been measured 100+ times, other have been measured once. The response value is numerical. On a small dataset the model seems to perform as I would have expected, however on the full data set, with ~400k observations, the model is simply to slow. After 12 hours the model has not even gone through 500 iterations in any of the four chains run in parallel, whereas in the JAGS the roughly equivalent model took just under a day with ~10,000 iterations. I have tried to follow the guidelines for code optimisation in the manual plus what I have found in this forum.

What I would like to model is,

  y = a0 + a_S + a_age + a_sex + a_sex*age + a_analysis + a_type + a_age*type + a_sex*type + a_age*sex*type

a_S is the coefficient for each subject, modelled in as suggested by Kruschke (2010) for repeated measure experiments. I have precomputed all of the interactions in R, since I could not find a clever way of assigning each coefficient in Stan its own prior without a nested loop.

I have added a gist to the model file, and any comments on the model would be most welcome, both in terms of optimization but also scientific interpretation.

https://gist.github.com/demodw/4fe7e274bcbb76c0353a

Thanks!

David

[Jonah]

The real Stan folks can give you a more complete answer, but I noticed a few things that you can change that might help things out: 

1. You can get rid of all the uniform(0,10) sampling statements and just declare the parameters with lower=0 and upper=10 constraints. Not sure how much this will help performance, but it certainly won't hurt. 

2. You can use the "Matt trick" for all the normal priors on the coefficients. For example, for a0, you would do 

parameters {

  real a0_raw ;

  real<lower=0, upper=10> a0sigma ; 

}

transformed parameters {

 real a0 ;

 a0 <- a0_raw * a0_sigma ; 

}

model {

 a0_raw ~ normal(0,1) ;

}

Making the sampling statements all for standard normals often helps quite a bit. 

[Jonah Sol Gabry]

In my example of the Matt trick I forgot to declare a0_raw and a0 as reals. I edited my original post to correct that. 

On Friday, January 30, 2015, Jonah <jga...@gmail.com> wrote:

The real Stan folks can give you a more complete answer, but I noticed a few things that you can change that might help things out: 

1. You can get rid of all the uniform(0,10) sampling statements and just declare the parameters with lower=0 and upper=10 constraints. Not sure how much this will help performance, but it certainly won't hurt. 

2. You can use the "Matt trick" for all the normal priors on the coefficients. For example, for a0, you would do 

parameters {

  a0_raw ;

}

transformed parameters {

 a0 ;

 a0 <- a0_raw * a0_sigma ; 

}

model {

 a0_raw ~ normal(0,1) ;

}

Making the sampling statements all for standard normals often helps quite a bit. 

On Friday, January 30, 2015 at 5:45:14 AM UTC-5, David Westergaard wrote:

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/s1uqTRUN12c/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

We're all "real Stan folks"!

If you declare y as a vector, then you can define
the standardized version in one line as

vector[Ntotal] zy;
zy <- (y - mean(y)) / sd(y);

We should really implement a z_score() or standardize()
function to do this in one function call.

We go over which priors we recommend for regressions
in the manual regression chapter. Short story is we
prefer softer constraints to hard intervals like you have.
(Warning: We haven't even updated all of our own examples to
follow our own advice yet!)

Is there a reason you're not giving the hierarchical
priors intercepts. (Then the Matt trick, aka non-centered
parameterization, requires adding the mean back in.)

Without actually going into details of all the indices
on the predictors and effects, the basic structure looks
OK.

- Bob

> On Jan 30, 2015, at 11:37 AM, Jonah <jga...@gmail.com> wrote:
>
> The real Stan folks can give you a more complete answer, but I noticed a few things that you can change that might help things out:
>
> 1. You can get rid of all the uniform(0,10) sampling statements and just declare the parameters with lower=0 and upper=10 constraints. Not sure how much this will help performance, but it certainly won't hurt.
>
> 2. You can use the "Matt trick" for all the normal priors on the coefficients. For example, for a0, you would do
>
> parameters {
> a0_raw ;
> }
> transformed parameters {
> a0 ;
> a0 <- a0_raw * a0_sigma ;
> }
> model {
> a0_raw ~ normal(0,1) ;
> }
>
> Making the sampling statements all for standard normals often helps quite a bit.
>
>
>
> On Friday, January 30, 2015 at 5:45:14 AM UTC-5, David Westergaard wrote:

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

[David Westergaard]

I managed to implement the "Matt trick", though I had to convert my data type from real to vectors. I am not really sure what impact this has on performance, but hopefully it is better than ~18 hours for 500 iterations.

Could you point me to the place on the reference manual where you discuss priors in regression?

With regard to no intercept on the hierarchical priors, are you referring to the mean value of 0 on all coefficients? I actually started out by having hyperpriors on the coefficient mean, but I removed them since I was doing some debugging. Then plan is definitely to add them back, though.

This is the "new" model with the "Matt trick": https://gist.github.com/demodw/87e42a6d583fea4651c4

Cheers,

David

[Rob Trangucci]

Section 5.2 in the manual discusses prior a for regression. 

--

Rob

[Bob Carpenter]

On Jan 30, 2015, at 1:59 PM, David Westergaard <da...@harsk.dk> wrote:
>
> I managed to implement the "Matt trick", though I had to convert my data type from real to vectors. I am not really sure what impact this has on performance, but hopefully it is better than ~18 hours for 500 iterations.

It's not the time per iteration that matters, but the time
per effective sample. The non-centered parameterization will
mainly give you more effective samples per iteration, though
it can also speed up the iterations by requiring fewer leapfrog
steps.

You may not even need 500 iterations depending on how well
the model mixes with your data.

> Could you point me to the place on the reference manual where you discuss priors in regression?

There are several sections in the regression chapter.

Tightening priors (even if it doesn't have a noticeable impact
on the posterior) can be a big help with computation. HMC
really likes to stretch its legs and explore when you give
it room to run.

> With regard to no intercept on the hierarchical priors, are you referring to the mean value of 0 on all coefficients? I actually started out by having hyperpriors on the coefficient mean, but I removed them since I was doing some debugging. Then plan is definitely to add them back, though.

Yes, that's what I meant.

- Bob

P.S. Rob's working on implementing marginal maximum likelihood, so we
should have a more general version of bglme soon.

[Andrew Gelman]

Matt trick yes.  But I say no on the <lower=0,upper=10> thing.  I much prefer soft constraints.

On Jan 30, 2015, at 11:37 AM, Jonah <jga...@gmail.com> wrote:

The real Stan folks can give you a more complete answer, but I noticed a few things that you can change that might help things out: 

1. You can get rid of all the uniform(0,10) sampling statements and just declare the parameters with lower=0 and upper=10 constraints. Not sure how much this will help performance, but it certainly won't hurt. 

2. You can use the "Matt trick" for all the normal priors on the coefficients. For example, for a0, you would do 

parameters {

  a0_raw ;

}

transformed parameters {

 a0 ;

 a0 <- a0_raw * a0_sigma ; 

}

model {

 a0_raw ~ normal(0,1) ;

}

Making the sampling statements all for standard normals often helps quite a bit. 

On Friday, January 30, 2015 at 5:45:14 AM UTC-5, David Westergaard wrote:

[David Westergaard]

Thanks a lot for the help everyone. Applying the "Matt trick" really made things go faster. 

After 5,000 iteration some of the R-hat coefficients are still between 1.3-1.4. It might be that I need better initialisation values, though. Right now I am just letting Stan set some at random.

This is the "final" code with priors on the coefficient intercept: https://gist.github.com/demodw/88552f402d82669bfef6

David 

[Bob Carpenter]

My guess is that with so many varying intercepts, you’re running
into identifiability issues:

mu[i] <- a0_raw + aS_raw[subj[i]] + a1_raw[age[i]]
+ a2_raw[sex[i]] + + a3_raw[analysis[i]] + a4_raw[type[i]]
+ a1a2_raw[sexage[i]] + a1a4_raw[agetype[i]] + a2a4_raw[sextype[i]]
+ a1a2a4_raw[agesextype[i]];

If agesextype is an interaction term, it’ll only make things
worse.

There’s a chapter in the manual on problematic posteriors with some illustrations
of what goes wrong when parameters are unidentified or only weakly identified.

I’m not quite sure how to address the problem here, though.
Anyone else have any ideas?

- Bob

[David Westergaard]

I am not quite sure i understand the problem. I am well aware that some of the intercepts will have minimal influence on the measurement value, and that is exactly the case. I have included some, plus their interactions, to show these have no effect. E.g. that the measurement value does not vary regardless of age and gender. 

If this is not the “correct" way of doing it, then I am open for suggestions. Would it be better to have a group wide prior on each nominal variable?

Cheers,

David

[Bob Carpenter]

Your model’s OK theoretically, but it presents
computational issues due to weak identifiability.

You can subtract a constant from a0_raw and add it
to any of the other a… vars, such as aS_raw and get exactly
the same likelihood. Therefore the likelihood’s not identified.

Adding normal(0,1) priors technically gives you identification,
but it’s weak. What are the scales of the parameters you’re getting?
It may help to standardize them.

There’s an example of just how this works and
how it slows things down versus a properly identified model in
the problematic posteriors chapter of the manual (it also shows
why a truly unidentified model won’t work in Stan even though it
may appear to be working [even though it isn’t] in BUGS/JAGS).

Also, adding means to the coefficients won’t work here as they
will also be unidentified.

- Bob

[Rob Trangucci]

Your model isn't well-identified as it's written, but it isn't because of a0_raw, but rather because of your inclusion of distinct means for each of your random effects. You should have one non-zero mean for all of your random effects, and no intercept *or* a global intercept and all zero-mean random effects. My suggestion would be to change the block:

|

 a0 <- a0_raw * a0sigma + a0_raw_mean;
  aS <- aS_raw * aSsigma + aS_raw_mean;
  a1 <- a1_raw * a1sigma + a1_raw_mean;
  a2 <- a2_raw * a2sigma + a2_raw_mean;
  a3 <- a3_raw * a3sigma + a3_raw_mean;
  a4 <- a4_raw * a4sigma + a4_raw_mean;
 
  a1a2 <- a1a2_raw * a1a2sigma + a1a2_raw_mean;
  a1a4 <- a1a4_raw * a1a4sigma + a1a4_raw_mean;
  a2a4 <- a2a4_raw * a2a4sigma + a1a2_raw_mean;
  a1a2a4 <- a1a2a4_raw * a1a2a4sigma + a1a2a4_raw_mean;
|

to:

|
 a0 <- a0_raw * a0sigma + a0_raw_mean;
  aS <- aS_raw * aSsigma;
  a1 <- a1_raw * a1sigma;
  a2 <- a2_raw * a2sigma;
  a3 <- a3_raw * a3sigma;
  a4 <- a4_raw * a4sigma;
 
  a1a2 <- a1a2_raw * a1a2sigma;
  a1a4 <- a1a4_raw * a1a4sigma;
  a2a4 <- a2a4_raw * a2a4sigma;
  a1a2a4 <- a1a2a4_raw * a1a2a4sigma;
|

while eliminating lines 127 through 135 and 81 through 89.

Also, for clarity's sake, you should be putting your upper bounds on your sigma parameters in the parameter block, like so:

|
real<lower=0,upper=100> a0sigma;
|

Rob

[David Westergaard]

I am a bit unsure if this is the way to go. As I mentioned, what I am really interested in with this model is to estimate the effect from the random effects on the measurement. E.g., is there an negative, zero or non-zero effect from age on the measurement value. If I add a random intercept to the baseline value (a0), I am not really sure if the suggested model is describing what I am actually interested in. On the other hand, the commonly accepted knowledge about these coefficients is actually that they should have close to zero effect on the measurement. 

What it make more sense if I split up the model into three, one for each phenotype? Currently, I have combined three (related) phenotypes (type) in the model, which leads to the three way interaction term (age*sex*type). 

Would it help to re-evaluate the priors on the standard deviation for each random effect? I can imagine there are more suitable priors for e.g. gender rather than just a uniform one.

On a somewhat unrelated note, is there a smart "trick" in STAN for converting these standardised nominal coefficients back to their original scale?

Cheers,

David

[Bob Carpenter]

> On Feb 6, 2015, at 7:55 AM, David Westergaard <da...@harsk.dk> wrote:
>
> I am a bit unsure if this is the way to go. As I mentioned, what I am really interested in with this model is to estimate the effect from the random effects on the measurement. E.g., is there an negative, zero or non-zero effect from age on the measurement value. If I add a random intercept to the baseline value (a0), I am not really sure if the suggested model is describing what I am actually interested in. On the other hand, the commonly accepted knowledge about these coefficients is actually that they should have close to zero effect on the measurement.
>
> What it make more sense if I split up the model into three, one for each phenotype? Currently, I have combined three (related) phenotypes (type) in the model, which leads to the three way interaction term (age*sex*type).

Usually you want to keep as many predictors in a model as you
can fit.

> Would it help to re-evaluate the priors on the standard deviation for each random effect? I can imagine there are more suitable priors for e.g. gender rather than just a uniform one.
>
> On a somewhat unrelated note, is there a smart "trick" in STAN for converting these standardised nominal coefficients back to their original scale?

You can define a function if you want, and apply it in the
generated quantities block. But there's nothing built-in for
this.

- Bob