Value-added model coding question (based on Mariano, McCaffrey, and Lockwood 2010)

[Christopher Candelaria]

Hello all,

My question is about the efficiency of my code. I am not sure if I am leveraging Stan properly to take advantage of all its speed to estimate a value-added model, especially in the presence of missing data in the dependent variable vector. 

Here is a brief explanation of what I am doing:

I have a few research projects that are using Stan to perform value-added estimation using the generalized persistence model that was developed by Louis Mariano, Dan McCaffrey, and J.R. Lockwood. The model was originally coded in WinBUGS, but I have translated it into Stan code. As I understand, very few people have attempted to estimate their model because it was prohibitively slow using WinBUGS. I started out using WinBUGS, but I could not estimate the model with a large data set. I have since transitioned to Stan and have had success. 

The model leverages the longitudinal nature of student test scores; however, in many cases, student test scores are missing. The model is layered so that one can account for the full history of teachers that a student has. Here is an example:

y_1stgrade = mu_1 + theta_11[grade1teacher] + error_1

y_2ndgrade = mu_2 + theta_12[grade1teacher] + theta_22[grade2teacher] + error_2

y_3rdgrade = mu_3 + theta_13[grade1teacher] + theta_23[grade2teacher] + theta33[grade3teacher] + error_3

In this model, y is the test score vector for a given student, mu is a vector of year fixed effects, and the theta variables are random effects or the effects that teachers had on student test scores. What is unique about this model is that we can obtain the effect of the teacher in the year he or she teachers as well as the effect that he or she has when the students is in a future grade. 

Overall, the model above and my explanation is a gross simplification of everything the model can do; I just wanted to give a brief summary. The original paper can be found using the following doi: 10.3102/1076998609346967

Back to my question...

I would appreciate any general feedback about my code, but I think I found a reason why my model tends to run slowly; specially, I am not sure if I am handling missing test score data values efficiently. Specifically, I have been filling in missing test score data values with parameter values as necessary. Please see lines 153 to 167 for specific details. I borrowed some code from Ben Goodrich to do this, but I am not sure if there is a better way of handling missing values given that Stan has been updated regularly over the past year. I will also say that I am trying to estimate this model using a large, district data set. 

Any feedback would be greatly appreciated.

Thanks,

Chris

P.S. I am also happy to share any other details that may be needed. FYI, I am using CmdStan given the complexity of the data. 

[Bob Carpenter]

The missing data handling is OK.

The covariance matrices are very slow. Models are much faster (and more
robust arithmetically) using Cholesky factors. We don't have the (inverse) Wishart
implemented for Cholesky factors because we use the LKJ prior on correlation
matrices and scale it. I think you'll find that more flexible in terms of prior
(choose degree of concentration for correlation matrix around the unit vector and
flexible scale modeling). It'll also be a whole lot faster.

You also want to vectorize the multi_normal calls.

for(i in 2:Ntch[5]) {
theta5[i] ~ multi_normal(rep_vector(0,ih_l[5]), Gamma5);
}

is more efficient as

theta5[2:Ntch[5]] ~ multi_normal(rep_vector(0,ih_l[5]), Gamma5);

because you only need to solve (invert and determinant) Gamma5 once.

I find using Gamma5 and Gamma05 as variables very confusing!

There are a lot of these --- it'll make a huge difference in speed.

Also vectorize this:

for(i in 1:N) {
...
y ~ multi_normal(u[i], Sigma);

Also things like rep_vector(0,ih_l[5]) can be declared as transformed
data --- it's a bit more efficient than creating and destroying it each
time it's used (but with just data, it's not too bad).

The loops like these can be made more efficient

for (i in 1:N) {
u[i,1] <- mu[1] + theta1[tchid[i,1],1];

as

u[, 1] <- mu[1] + theta1[tchid[, 1], 1];

but that's not going to work until I fix the includes for our next release.
But you can do this:

u[1:N, 1] <- mu[1] + theta1[tchid[1:N, 1], 1];

Or, you can transpose the represntation of u and tchid and just do

u[1] <- mu[1] + theta1[tchid[1], 1];

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

[Ben Goodrich]

On Saturday, January 30, 2016 at 6:42:21 PM UTC-5, Christopher Candelaria wrote:

Back to my question...

I would appreciate any general feedback about my code, but I think I found a reason why my model tends to run slowly; specially, I am not sure if I am handling missing test score data values efficiently. Specifically, I have been filling in missing test score data values with parameter values as necessary. Please see lines 153 to 167 for specific details. I borrowed some code from Ben Goodrich to do this, but I am not sure if there is a better way of handling missing values given that Stan has been updated regularly over the past year. I will also say that I am trying to estimate this model using a large, district data set. 

If you have multivariate normal priors on unknown parameters, I pretty much always do a multivariate Matt trick and it pretty much always seems to be fine. In other words, instead of having

theta ~ multi_normal(zeros, Gamma);

in the model block, put

theta <- L * z_theta;

in the transformed parameters block where L is a Cholesky factor such that L * L' = Gamma but you often don't need to form Gamma until the generated quantities block if at all and put

z_theta ~ normal(0,1); // implies theta ~ multi_normal_cholesky(zeros, L);

in the model block. Then you can compose L as the product of a diagonal matrix of standard deviations and the Cholesky factor of a correlation matrix that has an lkj_corr_cholesky() prior. You can do whatever you want for the prior on the standard deviations and it usually doesn't matter too much.

Ben

[Ben Goodrich]

On Saturday, January 30, 2016 at 6:42:21 PM UTC-5, Christopher Candelaria wrote:

Back to my question...

I would appreciate any general feedback about my code, but I think I found a reason why my model tends to run slowly; specially, I am not sure if I am handling missing test score data values efficiently.

Also, it looks as if you forgot to translate this line

theta6[1] ~ normal(0,0.001);

from BUGS. In Stan, the second argument to normal() is the standard deviation, not the precision.

Ben
 

[Jonah Gabry]

On Sunday, January 31, 2016 at 6:44:05 PM UTC-5, Ben Goodrich wrote:

Also, it looks as if you forgot to translate this line

theta6[1] ~ normal(0,0.001);

from BUGS. In Stan, the second argument to normal() is the standard deviation, not the precision.

It might not matter a whole lot, but instead of translating it to the equivalent normal(0, sqrt(1000)) a much smaller standard deviation is probably better.

[Christopher Candelaria]

Thanks so much, Bob. I added a few comments below and attached a new file. 

On Sunday, January 31, 2016 at 2:35:18 PM UTC-8, Bob Carpenter wrote:

The missing data handling is OK.

The covariance matrices are very slow. Models are much faster (and more
robust arithmetically) using Cholesky factors.  We don't have the (inverse) Wishart
implemented for Cholesky factors because we use the LKJ prior on correlation
matrices and scale it.  I think you'll find that more flexible in terms of prior
(choose degree of concentration for correlation matrix around the unit vector and
flexible scale modeling).  It'll also be a whole lot faster.

I am definitely interested in rethinking the current approach to the covariance matrices. I will look more into using Cholesky factors. I might seek some advice about this after I attempt to implement it. 

 

You also want to vectorize the multi_normal calls.  

  for(i in 2:Ntch[5]) {
        theta5[i] ~ multi_normal(rep_vector(0,ih_l[5]), Gamma5);
  }

is more efficient as

  theta5[2:Ntch[5]] ~ multi_normal(rep_vector(0,ih_l[5]), Gamma5);

because you only need to solve (invert and determinant) Gamma5 once.

This is a great idea; I have already changed my code to reflect this. Is this type of vectorizing now possible because of the new R indexing capabilities?

 

I find using Gamma5 and Gamma05 as variables very confusing!

There are a lot of these --- it'll make a huge difference in speed.

Yes, I agree; I need to better name the variables. I am using the "0" to indicate the covariance matrices for the one "non-existent" teacher in each grade that gets assigned missing student links. By the way, when you mention there are a lot of these, are you referring to multi-normal calls or is there something fundamentally wrong with the Gamma0# and Gamma#, where # stands for the teacher effects covariance matrix for grade #. Sorry, I just want to make sure I understood your statement correctly. 

 

Also vectorize this:

 for(i in 1:N) {
   ...
   y ~ multi_normal(u[i], Sigma);

I am not sure how to vectorize this statement, given that I am dealing with missing data for each student in the same loop. I feel as if I am missing something obvious. As a reminder, here is the full code block:

   /* Fill in missing test score data values with parameter values as necessary
    * Code borrowed from Ben Goodrich's post on Stan users list:
    * https://groups.google.com/d/msg/stan-users/y9BalO9ILF0/XmpNOwIVsEcJ
    */
    count <- 1;
    for(i in 1:N) {
        vector[Nyear] y;
        for(j in 1:Nyear) {
            if(R[i,j] == 1)
                y[j] <- Z[i,j];
            else {
                y[j] <- Z_miss[count];
                count <- count + 1;
            }
        }      
       /*
        * Specify the likelihood for the test scores: In transformed parameters, we
        * specified the expected values conditional on theta.
        * The conditional distribution is normal with covariance matrix Sigma.
        * Note that "y", as it appears here, is student i's vector of test scores.
        */
        y ~ multi_normal(u[i], Sigma);
    }

  

Also things like rep_vector(0,ih_l[5]) can be declared as transformed
data --- it's a bit more efficient than creating and destroying it each
time it's used (but with just data, it's not too bad).

This is good to know. I will change this. I haven't done it yet in the attached file.  

The loops like these can be made more efficient

 for (i in 1:N) {    
   u[i,1] <- mu[1] + theta1[tchid[i,1],1];

as

  u[, 1] <- mu[1] + theta1[tchid[, 1], 1];

but that's not going to work until I fix the includes for our next release.
But you can do this:

  u[1:N, 1] <- mu[1] + theta1[tchid[1:N, 1], 1];

Or, you can transpose the represntation of u and tchid and just do

  u[1] <- mu[1] + theta1[tchid[1], 1];

Thanks, this is really helpful. I approached it as follows:

    u[1:N,1] <- mu[1] + theta1[tchid[i,1],1];
    u[1:N,2] <- mu[2] + theta1[tchid[i,1],2] + theta2[tchid[i,2],1];
    u[1:N,3] <- mu[3] + theta1[tchid[i,1],3] + theta2[tchid[i,2],2] +
                      theta3[tchid[i,3],1];
    u[1:N,4] <- mu[4] + theta1[tchid[i,1],4] + theta2[tchid[i,2],3] +
                      theta3[tchid[i,3],2] + theta4[tchid[i,4],1];
    u[1:N,5] <- mu[5] + theta1[tchid[i,1],5] + theta2[tchid[i,2],4] +
                      theta3[tchid[i,3],3] + theta4[tchid[i,4],2] +
                      theta5[tchid[i,5],1];
    u[1:N,6] <- mu[6] + theta1[tchid[i,1],6] + theta2[tchid[i,2],5] +
                      theta3[tchid[i,3],4] + theta4[tchid[i,4],3] +
                      theta5[tchid[i,5],2] + theta6[tchid[i,6]];

 Finally, just to check, is this acceptable:

   /*
    * Specify the priors for the grade means, Normal 0, variance 50
    * Assume normal 0, variance 50 for grade means.
    * NOTE: In Stan, specify the standard deviation.
    */
    mu[1:Nyear] ~ normal(mu_p[1:Nyear], sqrt(50));

This block is lines 142 to 147 in the attached file...the original code was:

    for(i in 1:Nyear) {
        mu[i] ~ normal(mu_p[i], sqrt(50));
    }  

[Christopher Candelaria]

I do want to move in this direction. In the next few days, I will try to experiment with the Matt trick. I'm not an expert with this, but I do want to learn. I think it will greatly improve the code and the efficiency of the model. Thanks!

[Christopher Candelaria]

I completely missed this. I have changed the 2nd argument to a standard deviation now. Thanks for catching it. 

[Christopher Candelaria]

As Bob and Ben suggested, I'm interested in using Cholesky factors of covariance matrix that is scaled from the an LKJ prior on a correlation matrix. I have been reading up on this, but I want to make sure I implement it properly.

I want to provide a minimally informative prior on each of the Gamma matrices that I have in my model. Mariano, McCaffrey, and Lockwood assumed that the Gamma matrices for each grade were distributed Wishart with K_g +1 degrees of freedom (where K_g is the size of the matrix for a given grade) and a scale matrix where all the off-diagonal elements were 0.8 and the variance elements had the following pattern: 0.25 x 0.7^(t-g). The basic idea is that the variance will decline over time.

To provide an example, consider the Gamma matrix for third grade teachers in my model, which is called Gamma3 in my code. The matrix is of size 6, as I have test scores from grades 3 to 8. For this set of teachers, g = 3 and t takes values 3 to 8. The idea behind the scale matrix is to reflect prior research suggesting that the variance of teacher effects is largest in the year teachers teach and that the variance declines over time (i.e., in the years when teachers' students move on to subsequent grades). So, the first variance element would be 0.25, and the subsequent ones would be 0.175, 0.1225, 0.08575, etc.  

Is there a practical way to code this type of structure using the LKJ prior. I understand the properties of the LJK shape parameter (e.g., if eta >1, then the modal correlation matrix would be the identity matrix with the concentration increasing the larger eta is). I'm not sure this would be a problem; I would have to do some simulations to check. Nevertheless, I would like to be able to be minimally informative with the variances, using the structure described above.

Is the approach I want to use similar to what is described on page 157 on the Stan Reference 2.9.0 manual (which I pasted below)? I see how the covariance matrix is generated from the LJK correlation matrix in the transformed parameters section, but I am not sure how to deal with the variances I described above. Would I draw from a Cauchy distribution? 

And then, after I scale the matrix, I would then do a Cholesky decomposition and then use something like theta ~ multi_normal_cholesky(zeros, L)...better yet, I would use the Matt trick, as Ben suggested.

Thanks for the help!

parameters {
    vector[K] mu; // topic mean
    corr_matrix[K] Omega; // correlation matrix
    vector<lower=0>[K] sigma; // scales
    vector[K] eta[M]; // logit topic dist for doc m
    simplex[V] phi[K]; // word dist for topic k
}

transformed parameters {
    ... eta as above ...
    cov_matrix[K] Sigma; // covariance matrix
    for (m in 1:K)
        Sigma[m,m] <- sigma[m] * sigma[m] * Omega[m,m];
    for (m in 1:(K-1)) {
        for (n in (m+1):K) {
            Sigma[m,n] <- sigma[m] * sigma[n] * Omega[m,n];
            Sigma[n,m] <- Sigma[m,n];
        }
    }
}

model {
    mu ~ normal(0,5); // vectorized, diffuse
    Omega ~ lkj_corr(2.0); // regularize to unit correlation
    sigma ~ cauchy(0,5); // half-Cauchy due to constraint
    ... words sampled as above ...
}

On Sunday, January 31, 2016 at 3:41:30 PM UTC-8, Ben Goodrich wrote:

[Bob Carpenter]

The prior on scales (sqrt variances) is independent of the LKJ prior on the
correlation matrices. Just create the vector of scale parameters

vector<lower=0>[G] sigma;

and you can put whatever kind of prior you want on it, such as:

for (g in 1:G)
sigma[g] ~ ...some prior that depends on g...

I wasn't sure where the t came in for 0.25 * 0.7^(t - g).

And don't build up a covariance matrix --- you should only need
its Cholesky factor, which you can get with diag_pre_multiply(sigma, L_Omega)
where L_Omega is the Cholesky factor of the correlation matrix.

- Bob

[Christopher Candelaria]

I have been trying to optimizing my code using all the suggestions given, but I am running into problems using the new R-style indexing. (By the way, I am compiling the program using cmdstan-2.9.0.)

For example, I want to do the following:

    theta1[2:Ntch[1]] <- diag_pre_multiply(scale1, L_Gamma1)*z_theta1[2:Ntch[1]];                    
    theta2[2:Ntch[2]] <- diag_pre_multiply(scale2, L_Gamma2)*z_theta2[2:Ntch[2]];
    theta3[2:Ntch[3]] <- diag_pre_multiply(scale3, L_Gamma3)*z_theta3[2:Ntch[3]];
    theta4[2:Ntch[4]] <- diag_pre_multiply(scale4, L_Gamma4)*z_theta4[2:Ntch[4]];
    theta5[2:Ntch[5]] <- diag_pre_multiply(scale5, L_Gamma5)*z_theta5[2:Ntch[5]];

but I am getting the following compiling error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for: 

  matrix * vector[]

Instead, I have written the code inefficiently as:

    for(i in 2:Ntch[1]) {
        theta1[i] <- diag_pre_multiply(scale1, L_Gamma1)*z_theta1[i];
    }

    for(i in 2:Ntch[2]) {
        theta2[i] <- diag_pre_multiply(scale2, L_Gamma2)*z_theta2[i];
    }    

    for(i in 2:Ntch[3]) {
        theta3[i] <- diag_pre_multiply(scale3, L_Gamma3)*z_theta3[i];
    }        

    for(i in 2:Ntch[4]) {
        theta4[i] <- diag_pre_multiply(scale4, L_Gamma4)*z_theta4[i];

    }    

    for(i in 2:Ntch[5]) {

        theta5[i] <- diag_pre_multiply(scale5, L_Gamma5)*z_theta5[i];
    }    

Similar problem but with real+real[ ] being ill formed when I do the following:

    u[1:N,1] <- mu[1] + theta1[tchid[1:N,1],1];
    u[1:N,2] <- mu[2] + theta1[tchid[1:N,1],2] + theta2[tchid[1:N,2],1];
    u[1:N,3] <- mu[3] + theta1[tchid[1:N,1],3] + theta2[tchid[1:N,2],2] +
                      theta3[tchid[1:N,3],1];
...etc.

Finally, the last block of code that is causing issues is:

    z_theta1[2:Ntch[1]] ~ normal(0,1); // implies theta1[2:Ntch[1]] ~ multi_normal_cholesky(zeros, L);
    z_theta2[2:Ntch[2]] ~ normal(0,1);
    z_theta3[2:Ntch[3]] ~ normal(0,1);
    z_theta4[2:Ntch[4]] ~ normal(0,1);
    z_theta5[2:Ntch[5]] ~ normal(0,1);

I have attached the revised version of my program. All the suggestions thus far have been extremely helpful. I assume I must use "for" loops? Perhaps I am missing something simple.

Thanks!

On Sunday, January 31, 2016 at 2:35:18 PM UTC-8, Bob Carpenter wrote:

[Daniel Lee]

Hi Chris,

Just a general comment -- when posting and asking about syntax errors, it helps to provide the full error message. It's a lot easier to help you interpret the message when the context is available.

More below.

On Wed, Feb 10, 2016 at 5:19 AM, Christopher Candelaria <chris.ca...@gmail.com> wrote:

I have been trying to optimizing my code using all the suggestions given, but I am running into problems using the new R-style indexing. 

For example, I want to do the following:

    theta1[2:Ntch[1]] <- diag_pre_multiply(scale1, L_Gamma1)*z_theta1[2:Ntch[1]];                    
    theta2[2:Ntch[2]] <- diag_pre_multiply(scale2, L_Gamma2)*z_theta2[2:Ntch[2]];
    theta3[2:Ntch[3]] <- diag_pre_multiply(scale3, L_Gamma3)*z_theta3[2:Ntch[3]];
    theta4[2:Ntch[4]] <- diag_pre_multiply(scale4, L_Gamma4)*z_theta4[2:Ntch[4]];
    theta5[2:Ntch[5]] <- diag_pre_multiply(scale5, L_Gamma5)*z_theta5[2:Ntch[5]];

but I am getting the following compiling error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for: 

  matrix * vector[]

The error message is telling you that Stan doesn't know how to multiply a matrix (matrix) by an array of vectors (vector[]).

- diag_pre_multiply(scale1, L_Gamma1) results in a matrix.

- z_theta1[2:Ntch[1]] results in an array of vectors. The size of the array is (Ntch[1] - 2 + 1).

We haven't vectorized our multiplication operator to broadcast the matrix on the left to multiply by each of the vectors in the array on the right to result in a array of vectors. Vectorizing all our operators in this way is on our todo list, but it'll be a while before you see it in Stan.

Instead, I have written the code inefficiently as:

    for(i in 2:Ntch[1]) {
        theta1[i] <- diag_pre_multiply(scale1, L_Gamma1)*z_theta1[i];
    }

    for(i in 2:Ntch[2]) {
        theta2[i] <- diag_pre_multiply(scale2, L_Gamma2)*z_theta2[i];
    }    

    for(i in 2:Ntch[3]) {
        theta3[i] <- diag_pre_multiply(scale3, L_Gamma3)*z_theta3[i];
    }        

    for(i in 2:Ntch[4]) {
        theta4[i] <- diag_pre_multiply(scale4, L_Gamma4)*z_theta4[i];
    }    

    for(i in 2:Ntch[5]) {
        theta5[i] <- diag_pre_multiply(scale5, L_Gamma5)*z_theta5[i];
    }    

You should cache the diag_pre_multiply(scale1, L_Gamma1) outside the for loop.

 

Similar problem but with real+real[ ] being ill formed when I do the following:

    u[1:N,1] <- mu[1] + theta1[tchid[1:N,1],1];
    u[1:N,2] <- mu[2] + theta1[tchid[1:N,1],2] + theta2[tchid[1:N,2],1];
    u[1:N,3] <- mu[3] + theta1[tchid[1:N,1],3] + theta2[tchid[1:N,2],2] +
                      theta3[tchid[1:N,3],1];
...etc.

Sorry -- not following.

 

Finally, the last block of code that is causing issues is:

    z_theta1[2:Ntch[1]] ~ normal(0,1); // implies theta1[2:Ntch[1]] ~ multi_normal_cholesky(zeros, L);
    z_theta2[2:Ntch[2]] ~ normal(0,1);
    z_theta3[2:Ntch[3]] ~ normal(0,1);
    z_theta4[2:Ntch[4]] ~ normal(0,1);
    z_theta5[2:Ntch[5]] ~ normal(0,1);

You'll have to use a loop. We've vectorized the normal distribution to allow for arrays of reals, vectors, and row vectors, but we haven't recursively vectorized it so it will automatically recurse through arrays of those structures. If you're worried about computation time for those operations, no need. It won't affect the runtime in a noticeable way by looping here.

Daniel

[Bob Carpenter]

You can also solve some of these problems by redeclaring
the variables. For instance, coding that array of vectors
as a matrix (just make sure the dimensions are in the right order).

And if you have a matrix u on the left of ~, you can use

to_vector(u) ~ ...

which will be more efficient than

for (n in 1:N)
u[n] ~ ...

- Bob