Multivariate multilevel regression with second level predictors

[Mauricio Garnier-Villarreal]

Hi Stan masters

I have been working on a random change point model, that is working and we are finishing up the paper with it. 

Now we are working on the next extension, where at the same time I am estimating a linear growth for a second variable (biomarker), all for the same subjects. I can model both longitudinal models simultaneously. The challenge is how we want to relate the 2 longitudinal trends. We want to use the linear growth of the biomarker as a second level predictor of the random change point. This would be simple enough for me to code if the parameters were independent, but the parameters are multivariate.  

Here is my code so far, but I think I amo doing something wrong. I tried to follow the coding examples from Gelman and Hill (2006), but in those examples all the random parameters have a second order predictor, in my case only one has it, and the book example use the inverse wishart for the correlation while I am trying to use the lkj since the wishart bias the correlations.

Any help would be greatly appreaciate

Thanks

data{

int n; 

int parmt;

int nsubj;

int nvars;    

vector[nvars] y[n]; // arrray of order y[n,nvars]

int id[n];

real x[n]; 

vector[parmt] pm[nvars]; 

cov_matrix[parmt] cov_vb;

cov_matrix[nvars] cov_vm;

vector[nvars] mm;

int n_bio; // new

real x_bio[n_bio]; // new

real y_bio[n_bio]; // new

int id_bio[n_bio]; //new

}

transformed data{ // cholesky decompisition of prior covariances

cholesky_factor_cov[parmt,parmt] Lvb;

cholesky_factor_cov[nvars,nvars] Lvm;

Lvb <- cholesky_decompose(cov_vb);

Lvm <- cholesky_decompose(cov_vm);

}

parameters {

vector[parmt] b[nvars,nsubj]; // array of order b[nvars,nsubj,parmt]

corr_matrix[parmt] corr_b[nvars]; // correlation between bs, one cor matrix for each variable

corr_matrix[nvars] corr_y; // correlation between y variables

vector<lower=0>[parmt] sd_b[nvars];

vector<lower=0>[nvars] sd_y; //

real<lower=0> var_bio; // new

vector[2] b_bio[nsubj]; // new array of order b[nsubj,2]

vector[2] mub_bio; // new

real<lower=0> varb_bio[2]; // new

vector[2] gm[nvars];//

vector[parmt] mub[nvars]; // array of order mub[nvars,parmt]

}

transformed parameters {

vector[nvars] mu[n]; // vector of order mu[n,nvars]

cholesky_factor_cov[nvars,nvars] Ly;

cholesky_factor_cov[parmt,parmt] Lb[nvars];

vector[n_bio] mu_bio; // new

vector[nvars] difs[nsubj];

vector[parmt] mub_hat[nsubj,nvars]; // mub_hat[nsubj,nvars,parmt]

for(k in 1:nvars){

Lb[k] <- diag_matrix(sd_b[k]) * cholesky_decompose(corr_b[k]);}

Ly <- diag_matrix(sd_y) * cholesky_decompose(corr_y);

for(j in 1:n_bio){

mu_bio[j] <- b_bio[id_bio[j],1] + b_bio[id_bio[j],2]*x_bio[j];}

for (i in 1:n){ // this loops run the piecewise regression

for (k in 1:nvars){

mu[i,k] <- b[k,id[i],1] + b[k,id[i],2] * x[i];

if (x[i] > b[k,id[i],4])

mu[i,k] <- mu[i,k] + (b[k,id[i],3] * (x[i] - b[k,id[i],4] ));

}}

for(j in 1:nsubj){

for(k in 1:nvars){

mub_hat[j,k,1] <- mub[k,1]; 

mub_hat[j,k,2] <- mub[k,2]; 

mub_hat[j,k,3] <- mub[k,3]; 

mub_hat[j,k,4] <- gm[k,1] + gm[k,2] * b_bio[j,2]; 

}}

for(i in 1:nsubj){

for(j in 1:nvars){

difs[i,j] <- b[j,i,3]-b[j,i,2];}}

}

model{

for(k in 1:nvars){

gm[k] ~ normal(0, 100);

sd_b[k] ~ cauchy(0, 2.5);}

sd_y ~ cauchy(0, 2.5);

varb_bio ~ normal(0,5);

var_bio ~ normal(0,5);

corr_y ~ lkj_corr(1);

for(k in 1:nvars){

corr_b[k] ~ lkj_corr(1);}

for (j in 1:nsubj){

for (k in 1:nvars){

b[k,j] ~ multi_normal_cholesky(mub_hat[k], Lb[k]);

mub_hat[k] ~ multi_normal_cholesky(mub[k], Lvb);

}}

for (i in 1:n){

y[i] ~ multi_normal_cholesky(mu[i], Ly);

mu[i] ~ multi_normal_cholesky(mm, Lvm);}

y_bio ~ normal(mu_bio, var_bio);

for(j in 1:nsubj){

b_bio[j] ~ normal(mub_bio,varb_bio);}

mub_bio ~ normal(0, 10);

}

[Bob Carpenter]

I think this has a problem (in addition to quirky formatting!):

for (j in 1:nsubj){
for (k in 1:nvars){
b[k,j] ~ multi_normal_cholesky(mub_hat[k], Lb[k]);
mub_hat[k] ~ multi_normal_cholesky(mub[k], Lvb);
}}

Shouldn't that mub_hat[k] be in a loop over k --- as is, that
same sampling is getting repeated nsubj times.

I'm not an expert in multivariate modeling, but this also looks
suspicious:

for (i in 1:n){
y[i] ~ multi_normal_cholesky(mu[i], Ly);
mu[i] ~ multi_normal_cholesky(mm, Lvm);}

Do you really want that compound multi-normal behavior?

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

[Mauricio Garnier-Villarreal]

Thanks for your comments 

About the quirky formatting, I call that individual expression. Are there formatting guidelines? 

About this double loop, I code it originally when the multi_normal distribution wasnt completely vectorized, its great that now I can take out one of the loops

 for (j in 1:nsubj){ 
    for (k in 1:nvars){ 
      b[k,j] ~ multi_normal_cholesky(mub_hat[k], Lb[k]); 
      mub_hat[k] ~ multi_normal_cholesky(mub[k], Lvb); 
    }} 

About the compound multi-normal behavior. This is the only way I know to include priors for the expected values, should this be excluded? Examples in JAGS include this structure, is this not neccesary or recommended in Stan?

 for (i in 1:n){ 
    y[i] ~ multi_normal_cholesky(mu[i], Ly); 
    mu[i] ~ multi_normal_cholesky(mm, Lvm);} 

I am also not clear if this way of coding is right to predict a random parameter, 

for(j in 1:nsubj){

for(k in 1:nvars){

mub_hat[j,k,1] <- mub[k,1]; 

mub_hat[j,k,2] <- mub[k,2]; 

mub_hat[j,k,3] <- mub[k,3]; 

mub_hat[j,k,4] <- gm[k,1] + gm[k,2] * b_bio[j,2]; 

}}

for (k in 1:nvars){

b[k] ~ multi_normal_cholesky(mub_hat[k], Lb[k]);

mub_hat[k] ~ multi_normal_cholesky(mub[k], Lvb);

}

Thank you

[Bob Carpenter]

> On Jul 7, 2015, at 4:48 PM, Mauricio Garnier-Villarreal <mauga...@gmail.com> wrote:
>
> Thanks for your comments
>
> About the quirky formatting, I call that individual expression. Are there formatting guidelines?

Stan's not whitespace sensitive on purpose.
I don't think I've ever seen anyone end all the
open curly braces in a block like you do. No
space between "for" and "(" or between ")" and "{"
seems to be common in R programmers, but I don't
think I've ever seen it in C or Java code.

There are formatting guidelines, which are basically
what we use for our own C++ code and basically follow
Google's C++ formatting guidelines (the only real difference
being the outer blocks).

> About this double loop, I code it originally when the multi_normal distribution wasnt completely vectorized, its great that now I can take out one of the loops

The biggest upside is that it'll be a lot faster because
the matrix only needs to be factored once.

>
> for (j in 1:nsubj){
> for (k in 1:nvars){
> b[k,j] ~ multi_normal_cholesky(mub_hat[k], Lb[k]);
> mub_hat[k] ~ multi_normal_cholesky(mub[k], Lvb);
> }}
>
>
> About the compound multi-normal behavior. This is the only way I know to include priors for the expected values, should this be excluded? Examples in JAGS include this structure,

Do you have an example in mind? I think I was mostly confused
because I think that mub_hat should be on the outside --- as is,

> mub_hat[k] ~ multi_normal_cholesky(mub[k], Lvb);

is getting called j times. This would probably be caught by JAGS
when the data was read in if nvars > 1.

> is this not neccesary or recommended in Stan?
>
> for (i in 1:n){
> y[i] ~ multi_normal_cholesky(mu[i], Ly);
> mu[i] ~ multi_normal_cholesky(mm, Lvm);}

If mu[i] is the prior here, do you have an informative prior Lvm?
Don't use it just for something like a diagonal as it's hugely
inefficient.

>
>
> I am also not clear if this way of coding is right to predict a random parameter,

I'm not sure what you mean here. You can certainly declare
predictions as parameters and then just sample as with the
data. It's a little redundant compared to JAGS or BUGS letting
you just input NA values where you want predictions.

You can also do prediction more efficiently in the generated quantities
block. That way you get proper MC draws rather than relying on MCMC.

- Bob