Longitudinal Growth Study - Computation Time Lengthy...

[Allan Debertin]

Hi Stan Users,

I'm thinking of using STAN, with 'rstan' package instead of using 'lme4' for a longitudinal 'growth' study design, because eventually I would like to include missing information for the 'X2' variable below and STAN is less computational demanding than JAGS or BUGS. The computation time for the model seems little slow. I've started with simple models and built up. The more partial pooling variables I add, the much more computationally demanding the model becomes. I am assuming slopes are fixed in the linear regression,but I may code a model where slopes are variable among individuals, increasing the computation time further.

I wonder if there is a more efficient way of coding this example or if I am coding the model incorrectly? That being said, I have approximately 3000 individuals, so I think the slow computation time may be just the number of parameters that STAN needs to estimate.

Following Gelman and Hill 2006 the multilevel model is:

Y[i,j] ~ Normal(alpha[k[i,j] + Beta1[X1[i,j] + Beta2[X2[i,j] +Beta3(X1[i,j]*X2[i,j,), sigma^2), for i = 1,...,n subjects, for j = 1 to 2,m gender.

The 'lmer' function in 'lme4' package would code the above as:

Model1 <- lmer(Y ~ X1 + X2 + X1:X2 + (1 | subject) + (1 | gender) 

The Stan code I am using is:

data {

  int<lower=0> N; //no. of observations

  int<lower=0> n_subject; //no. of subjects

  int<lower=0> n_gender; //no. of genders

  int<lower=1, upper=n_subject> subject[N]; // label for each subject or individual

  int<lower=1, upper=n_sex> gender[N]; // grouping gender

  vector[N] u1; //fixed or predictator variable

  vector[N] u2; //fixed or predicator variable 2 

  vector[N] y; //Y 'growth' variable

}

parameters {

  real b1;

  real b2;

  real b3;

  vector[n_subject] a1;

  vector[n_gender] a2;

  real<lower=0, upper=100> sigma_a1;

  real<lower=0, upper=100> sigma_a2;

  real<lower=0, upper=100> sigma_y;

}

transformed parameters{

  vector[N] y_hat;

  vector[N] inter; interaction term

  

  for (i in 1:N)

    inter <- u1*u2

    y_hat[i] <- b1*u1[i] + b2*u2[i] + b3*inter[i] + a1[subject[i]] + a2[gender[i]];

}  

model {

b1 ~ normal(0, 100);

b2 ~ normal(0, 100);

b3 ~ normal(0, 100);    

a1 ~ normal(0, sigma_a1);

a2 ~ normal(0, sigma_a2);

growth ~ normal(y_hat, sigma_y);

}

Cheers,

Allan

[Andrew Gelman]

Don’t do real<lower=0, upper=100>.  Simply real<lower=0> is better.

Then add a prior for these parameters if you want.

Also, don’t do normal(0,100); better to do normal(0,10) or even normal(0,1).  I’m assuming your predictors have ben scaled reasonably.

A

--
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

>
> for (i in 1:N)
> inter <- u1*u2
> y_hat[i] <- b1*u1[i] + b2*u2[i] + b3*inter[i] + a1[subject[i]] + a2[gender[i]];

Can you send something that compiles?

u1 and u2 are both vectors, so you can't multiply
them like this. And you're missing semicolons and
brackets on the loop.

I'm not trying to be picky, I just don't know what
you intend with this.

The goal is to pull shared computations out and use
as much vectorization as you can, such as

y_hat <- b1 * u1 + b2 * u2 + b3 * inter;
for (n in 1:N)
y_hat[n] <- y_hat[n] + a1[subject[n]] + a2[gender[n]];

You might be able to go even further and compute the
n_subject x n_gender table of sums and double-index it
with a1_plus_a2[subject[n], gender[n]], but that's a pretty
extreme optimization and you probably shouldn't start there.

You can also define y_hat as a local variable in the model.
Unless you want to use it for inference later.

- Bob

> On Nov 25, 2015, at 7:03 PM, Allan Debertin <alland...@gmail.com> wrote:
>

[Allan Debertin]

Thanks Andrew and Bob, 

Also, don’t do normal(0,100); better to do normal(0,10) or even normal(0,1).  I’m assuming your predictors have been scaled reasonably.

The scaling for the predictors are reasonable, the coefficients for b1, b2, b3 are below 1, so I'll set the prior for them to normal(0,1). 

Can you send something that compiles? 

u1 and u2 are both vectors, so you can't multiply 
them like this.  And you're missing semicolons and 
brackets on the loop. 

Bob, do you want the stanfit file (200MB)? or just an example of some of the outputs for b1, b2 and b3 and error terms? I'll add the semicolons and brackets to the code. The goal of the u1*u2 is to describe an interaction term in the linear regression. I'm following the ARM examples, specifically Ch. 13, section 13.1 "radon_inter_vary.stan"  actually  I think I need to write

inter <-u1.*u2;

and then the loop for this multiplication isn't necessary. 

For this below

  y_hat <- b1 * u1 + b2 * u2 + b3 * inter; 
  for (n in 1:N) 
    y_hat[n] <- y_hat[n] + a1[subject[n]] + a2[gender[n]]; 

With this coding, would the b1, b2, b3 coefficients still be available as outputs from the stanfit? I may need to compare the coefficients among different predictor variables to see which have larger effects. 

Helpful comments and questions as always,

Allan

[Bob Carpenter]

I was just talking about having a model that compiles.
The first message had ill-formed Stan code like

for (i in 1:N)
inter <- u1*u2
y_hat[i] <- b1*u1[i] + b2*u2[i] + b3*inter[i] + a1[subject[i]] + a2[gender[i]];

This does an elementwise product:

> inter <-u1.*u2;

Stan saves only parameters, transformed parameters, and
generated quantities. With RStan you can select a subset of
those to save.

- Bob

[Allan Debertin]

Ah, right generated quantities, thanks Bob.

You are correct, the model wasn't compiling. 

Next step is to build a model that will work with missing values, the BUGS examples with missing data look helpful. I'll post that model when I have something that compiles. 

[Ryan Batt]

Allan,

I'm pretty new to Stan too. I was recently very surprised (and ultimately pleased) at how much of a speed-up I got from using non-centered parameterization. I had a multilevel modeel too, but maybe a bit more deeply nested. So YMMV.

By non-centered, I mean try something like

transformed parameters {
  b1 <- b1_mu + b1_raw * b1_sd;
}

parameters { 
  real b1_mu;
  real<lower=0> b1_sd;
  real b1_raw;
}

model { 
  b1_raw ~ normal(0, 1); // implies b1 ~ normal(b1_mu, b1_sd)
}

​

Then you just use b1 in your model as you normally would. Also, you can just drop the b1_mu if you want to use a mean of 0.

​

Such a procedure was a truly massive speed-up for me. It was shocking! Good luck.

-Ryan

[Michael Betancourt]

> Such a procedure was a truly massive speed-up for me. It was shocking! Good luck.

For some discussion why see http://arxiv.org/abs/1312.0906.

[Allan Debertin]

It turns out that the slowdown in computation was caused by the mislabeled data. I was using for "n_subject" was not numbered correctly, so instead of having  i = 1,...,10 subjects I had 8 subjects, with the label for subjects 7 and 8 missing. 

I tried the centering method you suggested, but YMMV. The computation time was worse with centering. However see below, I changed the entire modeling structure, I did some more reading on longitudinal modeling and I needed a multilevel linear model.

[Bob Carpenter]

This looks problematic:

for (j in 1:n_subject)
B_temp ~ multi_normal(transpose(row(B_hat, j)), Sigma_b);

That models B_temp n_subject times, each with a different mean
parameter corresponding to a row of B_hat.

I'm also not sure why you have a parameter named B_temp --- it doesn't
look like a temp as it's a saved parameter -- maybe it should be called mu_B?

- Bob


> On Dec 13, 2015, at 3:52 PM, Allan Debertin <alland...@gmail.com> wrote:
>
> Hi Stan Users and everyone that posted,
>
> I've done some more reading on longitudinal growth studies and came to the conclusion that I should be using a multilevel linear model to add grouping variables. I'd recommend reading "Applied Longitudinal Data Analysis: Modeling Change and Event Occurence" Singer and Willet 2003 for an introduction on longitudinal studies, particularly chapter 1 through 4.
>
> They also have example code for each of their chapters in various chapters http://www.ats.ucla.edu/Stat/r/examples/alda/default.htm
>
> I'm going to use a simpler model structure for an example below, but I'll try to follow it up with coding that includes an interaction with a second predictor variable.
>
> So following Gelman and Hill the new model structure is:
>
>
> Y[i,j] ~ Normal(alpha[j[i] + Beta[TIMEj[i], sigma_y^2), for i = 1,...,n subjects
> alpha[j] ~ Normal(gamma00 + gamma01m, sigma_alpha^2), for j = 1,...,J gender
> Beta[j] ~ Normal(gamma10 + gamma11m, sigma_Beta^2), for j = 1,...,J gender

>
>
> The 'lmer' function in 'lme4' package would code the above as:
>

> Model2 <- lmer(Y ~ TIME + (1 | gender/subject))
>
> The Stan code I am using:
>
> data {
> int<lower=0> N; //observations
> int<lower=0> n_subject; //number of individuals
> int<lower=0> n_gender; //number of gender
> int<lower=1, upper=n_subject> subject[N]; // # individuals
> int<lower=1, upper=n_gender> gender[N]; // # genders
> vector[N] time; //age of subject fixed or predictor variable
> vector[N] y; //response "growth" variable
> }
> parameters {
> // Level 1 parameters
> real<lower=0> sigma_bo; // sigma_alpha^2 from above
> real<lower=0> sigma_b1; // sigma_Beta^2 from above
> real<lower=0> sigma_y; // sigma_y^2
>
> // Level 2 parameters
> real g00;
> real g01;
> real g10;
> real g11;
>
> //covariance matrix
> real<lower=-1, upper=1> rho; // correlation term
> vector[2] B_temp;

>
>
> }
> transformed parameters{
> vector[N] y_hat;

> vector[n_subject] b0; //alpha coefficient
> vector[n_subject] b1; //BETA coefficient
> matrix[n_subject,2] B_hat;
> matrix[2,2] Sigma_b;
> matrix[n_subject,2] B;
>
> Sigma_b[1,1] <- pow(sigma_bo, 2);
> Sigma_b[2,2] <- pow(sigma_b1, 2);
> Sigma_b[1,2] <- rho * sigma_bo * sigma_b1;
> Sigma_b[2,1] <- Sigma_b[1,2];
>
> for (j in 1:n_subject){
> B_hat[j,1] <- g00 + g01*gender[j];
> B_hat[j,2] <- g10 + g11*gender[j];
> B[j,1] <- B_temp[1];
> B[j,2] <- B_temp[2];
> b0[j] <- B[j,1];
> b1[j] <- B[j,2];
> }
>
> for( i in 1:N){
> y_hat[i] <- b0[subject[i]] + b1[subject[i]]*time[i];
> }
> }
> model {
> g00 ~ normal(0, 1);
> g01 ~ normal(0, 1);
> g10 ~ normal(0, 1);
> g11 ~ normal(0, 1);
> rho ~ uniform(-1, 1);
> for (j in 1:n_subject){
> B_temp ~ multi_normal(transpose(row(B_hat,j)),Sigma_b);
> }
> y ~ normal(y_hat, sigma_y);

[Allan Debertin]

Hi everyone,

I finally had time to revisit the Longitudinal Growth Study using Stan. With some discussions, I ended up using a random-intercept for the random effect of individual in the study, instead of a random slope, random intercept and covariance structures. We were interested in the effect of gender/sex in the study, as well as the interactions between age and two covariates, we choose to include it as a factor. For our study, the R-code in 'lme4' package would be:

Model1 <-lmer(Y~F.sex+Age*C1*C2+(1|subject)

For STAN, the code would be:

data {

  int N; //number of data points

  real Y[N]; //growth

  real age[N]; //age

  real Fsex2[N]; //factor of gender coded as 1 for being sex "2" and 0 for not being sex "2"

  real C1[N]; //covariate 1

  real C2[N]; //covariate 2

  int subject_index[N];   //subject id

  int N_subject_index;

}

transformed data{

 real age_C1[N];  

 real age_C2[N]; 

 real C1_C2[N];    

 real age_C1_C2[N];

 for(i in 1:N){

 age_C1[i] <- age[i]*C1[i];

 age_C2[i] <- age[i]*C2[i];

 C1_C2[i] <- C1[i]*C2[i];

 age_C1_C2[i] <- age[i]*C1[i]*C2[i];

 }

}

parameters {

  real Intercept;   // Grand intercept

  real beta_age;    // could use a vector of betas instead, i.e. vector[9] Beta, but its easier to pick out each coefficient this way.

  real beta_Fsex2; 

  real beta_C1;

  real beta_C2;

  real beta_ageC1; // age_chum interaction

  real beta_ageC2; // age_npgo interaction

  real beta_C1C2;

  real beta_ageC1C2;

  real<lower=0> sigma;

  real vary_subject_index[N_subject_index];

  real<lower=0> sigma_subject_index;

} 

model {

  real vary[N];

  real glm[N];

    

  //priors

  Intercept ~ normal(0,10);  

  beta_age ~ normal(0,10);
  beta_Fsex2 ~ normal(0,10);

  beta_C1 ~ normal(0,10);

  beta_C2 ~ normal(0,10);

  beta_ageC1 ~ normal(0,10);

  beta_ageC2 ~ normal(0,10);

  beta_C1C2 ~ normal(0,10);

  beta_ageC1C2 ~ normal(0,10);

  sigma_subject_index ~ lognormal(0,1);  

  sigma ~ lognormal(0,1); 

    

  // Varying or Random effects   

  for ( j in 1:N_subject_index ) 

    vary_subject_index[j] ~ normal( 0 , sigma_subject_index );  

  

  // Fixed effects   

  for ( i in 1:N ) {      

    vary[i] <- vary_subject_index[subject_index[i]];      

    glm[i] <- vary[i] + Intercept 

    + beta_age * age[i] 

    + beta_Fsex2 * Fsex2[i]           

    + beta_C1*C1[i]

    + beta_C2*C2[i]

    + beta_ageC1*age_C1[i]

    + beta_agec2*age_C2[i] 

    + beta_C1C2*C1_C2[i]

    + beta_ageC1C2*age_C1_C2[i];   

  }  

  Y ~ normal(glm, sigma);

}

generated quantities{

  real vary[N];

  real glm[N];

  vector [N] log_lik;

  real dev;

  dev <- 0;

  for ( i in 1:N ) {  

    vary[i] <- vary_subject_index[subject_index[i]];      

    glm[i] <- vary[i] + Intercept 

     vary[i] <- vary_subject_index[subject_index[i]];      

    glm[i] <- vary[i] + Intercept 

    + beta_age * age[i] 

    + beta_Fsex2 * Fsex2[i]           

    + beta_C1*C1[i]

    + beta_C2*C2[i]

    + beta_ageC1*age_C1[i]

    + beta_agec2*age_C2[i] 

    + beta_C1C2*C1_C2[i]

    + beta_ageC1C2*age_C1_C2[i]; 

    log_lik[i] <- normal_log( growth[i] , glm[i] , sigma );

    dev <- dev + (-2)* normal_log(growth[i] , glm[i] , sigma);

  }

}

[Allan Debertin]

I should also mention that my data had 3500 (J) individuals and about 13000 data points (i.e. N). So using STAN  and HMC is likely a little overkill for this analysis. I ran 500 iterations and used 3 chains, I got covergence with this smaller amount iterations. Without using parallel processors (I am still having problems the new parallel processing commands), the model ran for 2 hours! Takes about a second to run in lme4, but kinda of expected given the difference in methods.

In the end, got the same estimate cofficients and estimated Std. Error (lme4) and Std. Deviation (MCMC chain in STAN)

[Bob Carpenter]

You can make all that a bit faster.

If you convert your data types from arrays to vectors, you
can replace this

for(i in 1:N) {
age_C1[i] <- age[i]*C1[i];

with

age_C1 <- age .* C1;

This is more to prevent errors (e.g., bounds checks built in,
no indexes to fiddle) than for efficiency.


The big win will be expressing the linear predictor as
a vectorized operation. First, glm and everything else needs
to be defined as a vector rather than array. Then you can
replace

for ( i in 1:N ) {
vary[i] <- vary_subject_index[subject_index[i]];
glm[i] <- vary[i] + Intercept
+ beta_age * age[i]
+ beta_Fsex2 * Fsex2[i]
+ beta_C1*C1[i]
+ beta_C2*C2[i]
+ beta_ageC1*age_C1[i]
+ beta_agec2*age_C2[i]
+ beta_C1C2*C1_C2[i]
+ beta_ageC1C2*age_C1_C2[i];
}

with

glm <- Intercept
+ vary_subject_index[subject_index]
+ beta_age * age
+ beta_Fsex2 * Fsex2
+ beta_C1 * C1
+ beta_C2 * C2
+ beta_ageC1 * age_C1
+ beta_agec2 *age_C2
+ beta_C1C2 * C1_C2
+ beta_ageC1C2 * age_C1_C2;

You could even just put this all right inside normal, but defining
the variable and assigning isn't that expensive (cheap copy).

That still has a bunch of multiply-by-zero you don't need in
the Fsex2 case, which is a boolean indicator, but there's not
a great way to speed that up without writing a function that either
inserts aconstant 0 if Fsex2[n] is 0 and inserts beta_Fsex2 otherwise.
No useless multiplications by 1 or 0 that way.

vector fill_if(int[] include, real beta) {
vector[num_elements(vector)] v;
real z;
z <- 0;
for (i in 1:num_elements(vector))
if (include[i])
v[i] <- beta;
else
v[i] <- z;
}
}

and then replace

+ beta_Fsex2 * Fsex2

with

+ fill_if(Fsex2, beta_Fsex2)

It's goofy to define the zero and reuse it, but it's a bit more
efficient with memory than using v[i] <- 0, which creates a new
autodiff variable for each instance.

Finally, you want to use a non-centered parameterization and
replace this:

for ( j in 1:N_subject_index )
vary_subject_index[j] ~ normal( 0 , sigma_subject_index );

by making vary_subject_index a vector and using

vary_subject_index_raw ~ normal(0, 1);

and then replacing the line above

+ vary_subject_index[subject_index]

with

+ sigma_subject_index * var_subject_index_raw[subject_index]

- Bob

> On Feb 19, 2016, at 1:43 PM, Allan Debertin <alland...@gmail.com> wrote:
>
> I should also mention that my data had 3500 (J) individuals and about 13000 data points (i.e. N). So using STAN and HMC is likely a little overkill for this analysis. I ran 500 iterations and used 3 chains, I got covergence with this smaller amount iterations. Without using parallel processors (I am still having problems the new parallel processing commands), the model ran for 2 hours! Takes about a second to run in lme4, but kinda of expected given the difference in methods.
>
> In the end, got the same estimate cofficients and estimated Std. Error (lme4) and Std. Deviation (MCMC chain in STAN)
>

[Allan Debertin]

Hi Bob,

Thanks for advice, I was able to figure out a few of the edits you mentioned, make sense.

Question about this bit:

That still has a bunch of multiply-by-zero you don't need in
the Fsex2 case, which is a boolean indicator, but there's not
a great way to speed that up without writing a function that either
inserts aconstant 0 if Fsex2[n] is 0 and inserts beta_Fsex2 otherwise.
No useless multiplications by 1 or 0 that way.

  vector fill_if(int[] include, real beta) {
    vector[num_elements(vector)] v;
    real z;
    z <- 0;
    for (i in 1:num_elements(vector))
      if (include[i])
        v[i] <- beta;
      else
        v[i] <- z;
    }
  }

This code would be a put in the stan function, parameters or model sections? I'm guessing its a STAN function.

and then replace

         + beta_Fsex2 * Fsex2

with

         + fill_if(Fsex2, beta_Fsex2)

Okay Got this part, makes sense.

[Bob Carpenter]

Yes, that funciton would go in the functions block.
I haven't tested it or even tried to compile it. I'd
only do this kind of thing if you were running the model multiple
times and the speed was very important. You might not even
notice the difference here with the fill_if() function because
there's so much other multiplication going on that it can't be
more than 5% or 10% of the total compute 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.

> To post to this group, send email to stan-...@googlegroups.com.