Very inefficient hierarchical model - trouble vectorizing multi-level IRT with correlated residuals

[Shauna Sweet]

Hello Bob and all,

I feel like something of a broken record... Because I keep running into similar problems: inefficiency!

I have successfully fit a multi-level Rasch model with regressions at multiple levels, and it worked extremely well. I employed all of the suggestions I have received thus far, including vectorizing my code, separating out the calculation of means as transformed parameters versus including them in sampling statements, using the if_else construction with design variables instead of relying on matrix calculations, and adjusting my priors.

My next step is to fit a 2 PL model. My concern is that while I have code that runs, it is frustratingly inefficient. All of the examples that I have seen and reviewed on this site, in assorted documentation, etc. parameterize a 2 PL model by specifying some version of vector[2] xi[K]; where xi is the vector of alpha-beta pairs over all of K items.

When this specification is extended into multiple levels, unfortunately, the model becomes... Really bogged down. My code successfully runs, but extremely slowly, and I'm looking for any suggestions as to how I can vectorize my code.

The problem that I seem to be running into is that I'm permitting the alphas and betas to be correlated. Right now, each alpha beta pair is drawn from a normal distribution with some (group-specific) mean vector and (group-specific) correlation. Unfortunately, when I attempt to vectorize this portion over the code so that I don't have to loop over items, I run into errors with the types of data structures permitted with multi_normal_cholesky. Also, because xi is a vector of length two in an array the length of all the items, I am also having difficulty vectorizing at level 1, and so I'm looping over observations.

Please help? I'm wondering if folks can suggest something - or if my weeks of working on this just means I'm too close to the problem and can't see the forest through the trees.

My model code is attached, and also pasted below.

Thanks so much,
Shauna



model<- "

data {
   int<lower=1> N;            // number of responses

   int<lower=1> J;            // number of students
   int<lower=1> K;            // number of items
   int<lower=1> P;            // number of families
   int<lower=1> M;            // number of item models

   int peeps[N];                  // student giving response n
   int<lower=1> item[N];         // item for response n
   int<lower=0,upper=1> score[N];   // correctness for response n

   int<lower=0,upper=J> pid[J];   // Person ID number
   int<lower=0,upper=K> iid[K];   // Item ID number
   int<lower=0,upper=P> fid[P];   // family ID number 
   int<lower=0,upper=M> mid[M];   // item model ID

   int<lower=1,upper=P> parent[K];         //indexes items to families
   int<lower=1,upper=M> imodel[P];    //indexes families to item models

   int<lower=0,upper=1> X1[P];                   
   int<lower=0,upper=1> X2[P];
}

transformed data{

   real noval;
   noval<-0;

}

parameters {

   vector[J] theta;                  // abilities

   vector[2] mum[M];
   vector[2] mup[P];

   real bxx1a;
   real bxx1b;
   real bxx2a;
   real bxx2b;

   vector<lower=0>[2] taup[P];             // vector for family-level alpha/beta residual sds

   vector<lower=0>[2] taui;              // vector for alpha/beta residual sds
   cholesky_factor_corr[2] L_Omega;    // cholesky factorization of residual correlation matrix

   vector[2] xi[K];              // alpha/beta pair vectors

}

transformed parameters {

   vector[2] mupm[P];             //parent means

   for (p in 1:P){

      mupm[p,1]<- mum[imodel[p],1] + if_else(X1[p], bxx1a, noval) + if_else(X2[p], bxx2a, noval);
      mupm[p,2]<- mum[imodel[p],2] + if_else(X1[p], bxx1b, noval) + if_else(X2[p], bxx2b, noval);

   }

}

model {

   mum[1,1] ~ normal(0,10);
   mum[2,1] ~ normal(0,10);
   mum[3,1] ~ normal(0,10);
   mum[1,2] ~ normal(0,10);
   mum[2,2] ~ normal(0,10);
   mum[3,2] ~ normal(0,10);

   bxx1a ~ normal(0,10);
   bxx1b ~ normal(0,10);
   bxx2a ~ normal(0,10);
   bxx2b ~ normal(0,10);

   taup[1] ~ inv_gamma(.1,.1);
   taup[2] ~ inv_gamma(.1,.1);

   taui[1] ~ inv_gamma(.1,.1);
   taui[2] ~ inv_gamma(.1,.1);

   L_Omega ~ lkj_corr_cholesky(4);

   theta ~ normal(0, 1);

   for (p in 1:P){
      mup[p,1] ~ normal(mupm[p,1], taup[p,1]);
      mup[p,2] ~ normal(mupm[p,2], taup[p,2]);
   }

   for (k in 1:K){
      xi[k] ~ multi_normal_cholesky(mup[parent[k]], diag_pre_multiply(taui,L_Omega));
   }

   for (n in 1:N){
      score[n] ~ bernoulli_logit(xi[item[n],1] * (theta[peeps] - xi[item[n],2]));
   }

}

generated quantities {

   corr_matrix[2] Omega;

   vector[K] alpha;
   vector[K] beta;

   vector[N] log_lik;
   real log_lik0;
   real dev;

   Omega <- multiply_lower_tri_self_transpose(L_Omega);

   dev <- 0;

   for (n in 1:N){
      log_lik[n] <- bernoulli_logit_log(score[n], xi[item[n],1] * (theta[peeps[n]] - xi[item[n],2]));
      dev <- dev + (-2 * log_lik[n]);
   }

   log_lik0 <- sum(log_lik); // global log-likelihood

   for (k in 1:K) {
      alpha[k] <- exp(xi[k,1]);
      beta[k] <- xi[k,2];
   }

}

"

[Bob Carpenter]

There are two kinds of efficiency --- the time per iteration
and the effective sample size per iteration. All that vectorization,
etc. do is deal with the first problem. To deal with the second you
need to look at reparameterizations like the non-centered transform.
You can see in the manual how to do that with a multivariate normal.

The main slow thing in your code is this:

for (k in 1:K){
xi[k] ~ multi_normal_cholesky(mup[parent[k]], diag_pre_multiply(taui,L_Omega));
}

Two big issues here. First, you don't ever want to recompute
the same thing twice. Compute it once, assign it to a local
variable, and reuse. Here, you compute diag_pre_multiply(taui,L_Omega)
each time in the loop. What you want to do is this:

xi ~ multi_normal_cholesky(mup[parent], diag_pre_multiply(taui, L_Omega));

That only computes the product once and only factors the matrix once.
It should be MUCH faster. You can vectorize the rest of the code
the same way, but it won't make as much difference because it
doesn't involve a matrix factorization in the inner loop.

- 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.
> For more options, visit https://groups.google.com/d/optout.
> <Correlated Residuals - Hierarchical IRT Model - SJS.R>

[Shauna Sweet]

Hi Bob and all -

Thank you so much! I made that simple change in the code (I attempted to remove other loops but that didn't work), and I also adjusted some of my code to ensure that I was actually specifying a prior for all of my parameters (I had to loop across k values of taup). The time to work through 1500 draws was shortened by a factor of 5 to approximately 20 minutes. I didn't think that the efficiency gains would be *that* significant?!

Regarding the second type of efficiency - I will write back when I work through some of the reparameterizing. I assume that the second kind of efficiency will improve my estimates of alpha, beta, and Omega? Right now, my effective sample size even for 1500 draws is only 3 for some items!

Beyond reparameterization, is there anything I can do, i.e. adjust my priors, to improve precision of the estimates? All other estimates look good, and I'm able to reproduce simulated values.

Thanks again!

Shauna

--
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/UjnUqlOi7wc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Shauna Sweet]

Hello Bob and all,

Below is the newly revised code that is currently running, with the changes highlighted in boldface type. In short, I reparameterized xi ~ multi_normal_cholesky(mup[

parent], diag_pre_multiply(taui,L_Omega)); according to the examples in the manual.

The questions that I have are:

(1) would a reparameterization of the priors on taui and taup also help with efficiency? I do see that it's possible to reparameterize an inv_gamma draw using the gamma distribution, and I'm not sure if that would be helpful here?

(2) is the reparameterization that I've used below where I can hope to gain the most in terms of information? Unlike the examples in the manual, I don't know anything about Sigma. Both tau (the scale) and L_Omega (the correlation) are unknown. That said, where I seemed to be having difficulty with rejection was at the item-parameter level, so I would assume draws of taui.

#######################################################
##### MODEL
#######################################################

   vector[2] xi_r[K];              // alpha/beta pair vectors

}

transformed parameters {

   vector[2] mupm[P];             //parent means

   vector[2] xi[K];

   for (p in 1:P){

      mupm[p,1]<- mum[imodel[p],1] + if_else(X1[p], bxx1a, noval) + if_else(X2[p], bxx2a, noval);
      mupm[p,2]<- mum[imodel[p],2] + if_else(X1[p], bxx1b, noval) + if_else(X2[p], bxx2b, noval);

   }

   for (k in 1:K){

      xi[k] <- mup[parent[k]] + taui .* (L_Omega * xi_r[k]);
   }
}

model {
  

   for (k in 1:K){

      xi_r[k]~normal(0,1);
      //xi ~ multi_normal_cholesky(mup[parent], diag_pre_multiply(taui,L_Omega));
   }
 
   mum[1,1] ~ normal(0,3);
   mum[2,1] ~ normal(0,3);
   mum[3,1] ~ normal(0,3);
   mum[1,2] ~ normal(0,3);
   mum[2,2] ~ normal(0,3);
   mum[3,2] ~ normal(0,3);

   bxx1a ~ normal(0,3);
   bxx1b ~ normal(0,3);
   bxx2a ~ normal(0,3);
   bxx2b ~ normal(0,3);

   for (p in 1:P){

      taup[p,1] ~ inv_gamma(.1,.1);
      taup[p,2] ~ inv_gamma(.1,.1);

   }

   taui[1] ~ inv_gamma(.1,.1);
   taui[2] ~ inv_gamma(.1,.1);

   L_Omega ~ lkj_corr_cholesky(4);

   theta ~ normal(0, 1);

   mup[1] ~ normal(mupm[1], taup[1]);
   mup[2] ~ normal(mupm[2], taup[2]);

   for (n in 1:N){

      score[n] ~ bernoulli_logit(xi[item[n],1] * (theta[peeps[n]] - xi[item[n],2]));

[Bob Carpenter]

You can replace this:

for (k in 1:K){
xi_r[k]~normal(0,1);

with just

to_vector(xi_r) ~ normal(0, 1);

for a bit more savings. It requires making xi_r
a matrix.

Same with this:

> for (p in 1:P){
> taup[p,1] ~ inv_gamma(.1,.1);
> taup[p,2] ~ inv_gamma(.1,.1);
> }

which can just be

to_vector(taup) ~ inv_gamma(0.1, 0.1);

But that would require defining taup as
and xi_r as matrices.

BUT, we don't recommend these inverse gamma priors.

You can also replace these

> mup[1] ~ normal(mupm[1], taup[1]);
> mup[2] ~ normal(mupm[2], taup[2]);

with

mup ~ normal(mupm, taup);

You can also vectorize this:

> for (n in 1:N){
> score[n] ~ bernoulli_logit(xi[item[n],1] * (theta[peeps[n]] - xi[item[n],2]));
> }

score ~ bernoulli_logit(xi[item, 1] .* (theta[peeps] - xi[item, 2]));

[Shauna Sweet]

Hello all,

Thank you so much!

I have worked to make all of the suggested changes, but turning one thing into a matrix snowballed, and now I'm lost inside my own code. I apologize because my intent is not to use thee list to do troubleshooting, but rather to have specific questions - but I guess here my question is: do I need to make ALL of these objects matrices?

If I understand the to_vector function, I need my alpha-beta pairs to be in each column, so my matrices should all have two rows and K, P, etc. columns, correct?

I apologize for being quite lost. The code compiles - meaning now all of my objects are "legal", but the dimensions are off in the following statement and I can't seem to get things righted: xi <- mup[,parent] + taui * (L_Omega * xi_r); ... Any suggestions would be greatly appreciated.

Is there space somewhere for a model example that shows progressive efforts to add levels, parameters, and then optimize the code? Are you looking for such examples? I would be happy to clean up and annotate my code, etc. if you feel it is appropriate. I've received a lot of help from folks in this group.

S

model<- "

data {
int<lower=1> N;            // number of responses

int<lower=1> J;            // number of students
int<lower=1> K;            // number of items
int<lower=1> P;            // number of families
int<lower=1> M;            // number of item models

int peeps[N];                  // student giving response n
int<lower=1> item[N];         // item for response n
int<lower=0,upper=1> score[N];   // correctness for response n

int<lower=0,upper=J> pid[J];   // Person ID number
int<lower=0,upper=K> iid[K];   // Item ID number
int<lower=0,upper=P> fid[P];   // family ID number 
int<lower=0,upper=M> mid[M];   // item model ID

int<lower=1,upper=P> parent[K];         //indexes items to families
int<lower=1,upper=M> imodel[P];    //indexes families to item models

int<lower=0,upper=1> X1[P];                   
int<lower=0,upper=1> X2[P];
}

transformed data{

real noval;
noval<-0;

}

parameters {

vector[J] theta;                  // abilities

vector[2] mum[M];

matrix[P,2] mup;
//vector[2] mup[P];

real bxx1a;
real bxx1b;
real bxx2a;
real bxx2b;

matrix<lower=0>[2,K] taup;
//vector<lower=0>[2] taup[P];             // vector for family-level alpha/beta residual sds

matrix<lower=0>[1,2] taui;
//vector<lower=0>[2] taui;              // vector for alpha/beta residual sds

cholesky_factor_corr[2] L_Omega;    // cholesky factorization of residual correlation matrix

matrix[2,K] xi_r;
//vector[2] xi_r[K];              // alpha/beta pair vectors

}

transformed parameters {

matrix[2,P] mupm;
matrix[2,K] xi;
//vector[2] mupm[P];             //parent means
//vector[2] xi[K];

for (p in 1:P){

mupm[1,p]<- mum[imodel[p],1] + if_else(X1[p], bxx1a, noval) + if_else(X2[p], bxx2a, noval);
mupm[2,p]<- mum[imodel[p],2] + if_else(X1[p], bxx1b, noval) + if_else(X2[p], bxx2b, noval);

}


xi <- mup[,parent] + taui * (L_Omega * xi_r);


//   for (k in 1:K){
//      xi[k] <- mup[parent[k]] + taui .* (L_Omega * xi_r[k]);
//   }

}

model {

to_vector(xi_r) ~ normal(0, 1);

//for (k in 1:K){
//  xi_r[k]~normal(0,1);

  //xi ~ multi_normal_cholesky(mup[parent], diag_pre_multiply(taui,L_Omega));

//}

mum[1,1] ~ normal(0,3);
mum[2,1] ~ normal(0,3);
mum[3,1] ~ normal(0,3);
mum[1,2] ~ normal(0,3);
mum[2,2] ~ normal(0,3);
mum[3,2] ~ normal(0,3);

bxx1a ~ normal(0,3);
bxx1b ~ normal(0,3);
bxx2a ~ normal(0,3);
bxx2b ~ normal(0,3);

to_vector(taup) ~ cauchy(0,5);

//for (p in 1:P){
//taup[p,1] ~ inv_gamma(.1,.1);
//taup[p,2] ~ inv_gamma(.1,.1);
//}

to_vector(taui) ~ cauchy(0,5);
//taui[1] ~ inv_gamma(.1,.1);
//taui[2] ~ inv_gamma(.1,.1);

L_Omega ~ lkj_corr_cholesky(4);

theta ~ normal(0, 1);

to_vector(mup) ~ normal(to_vector(mupm),to_vector(taup));

//mup[1] ~ normal(mupm[1], taup[1]);
//mup[2] ~ normal(mupm[2], taup[2]);

score ~ bernoulli_logit(xi[item,1] .* (theta[peeps] - xi[item,2]));

//for (n in 1:N){
//score[n] ~ bernoulli_logit(xi[item[n],1] * (theta[peeps[n]] - xi[item[n],2]));
//}

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/UjnUqlOi7wc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

Sorry you're lost.

My suggestion is to back out the changes and not try
to optimize until you're more comfortable with Stan's
data structures. And then to do optimization one item
at a time --- don't try to do them all at once.

The key is that all the argument types need to be appropriate
for the functions/operators to which they're applied. And
it all has to be consistent. There's a chapter on arrays
vs. matrix/vector types in the manual and I'd suggest starting
there to get a better feeling for how Stan's data structures work.
All of the functions are indexed with their argument types.

Feel free to post examples to this mailing list.
The problem we're running into now is that we have so much
doc and so many mailing list posts that people are having
a hard time finding relevant content and they're overwhelmed
by the hundreds of pages we already have.

- Bob

P.S. I'll do this one for you as an example. You say you're having
problems with:

> xi <- mup[,parent] + taui * (L_Omega * xi_r);

>

The thing to do is break down all the data types here
and see what they look like:

> matrix[2,K] xi;

> matrix[P,2] mup;
> matrix<lower=0>[1,2] taui;
> cholesky_factor_corr[2] L_Omega;
> matrix[2,K] xi_r;

And then look at the expression types bottom-up:

L_Omega * xi_r : mat[2,P]

taui * (L_Omega * xi_r) : mat[1,P]

mup[,parent] : vector[P]

mup[,parent] + taui * (L_Omega * xi_r): ***ill-formed***

There's your problem. You're trying to add a vector[P]
to a mat[1,P], which don't match.

You almost never want a matrix with a dimension of size 1.