Warnings (?) from a varying-intercept, varying-slope linear model with a categorical predictor

[Xiao He]

Dear all, 

I recently started trying out Stan via Rstan. I built a model like below based on this guide.

data {
 int<lower=0> N;                     //total # of observations
 int<lower=0> NS;                    //total # of subjects
 int<lower=0> NI;                    //total # of experimental items
 int<lower=1, upper=NS> subject[N];  //Subject ID
 int<lower=1, upper=NI> item[N];     //Item ID
 vector[N] cond;                     //Categorical predictor (3 levels)
 vector[N] rt;                       //Outcome variable (reaction time)
 vector[2] zero;                     //A vector of 0's
}


parameters {
 vector[2] beta;              //intercept and slope
 vector[2] u[NS];             //subject random effects
 vector[2] w[NI];             //item random effects
 real<lower=0> sigma_e;       //residual sd
 vector<lower=0>[2] sigma_u;  //subject sd
 vector<lower=0>[2] sigma_w;  //item sd
 corr_matrix[2] omega_u;      //correlation matrix for subject random effects
 corr_matrix[2] omega_w;      //correlation matrix for item random effects
}

transformed parameters{
 cov_matrix[2] Sigma_u;       //VCV matrix for subject random effects
 cov_matrix[2] Sigma_w;       //VCV matrix for item random effects
 for(r in 1:2){
   for(c in 1:2){
     Sigma_u[r, c] <- sigma_u[r] * sigma_u[c] * omega_u[r, c];
     Sigma_w[r, c] <- sigma_w[r] * sigma_w[c] * omega_w[r, c];
   }
 }
}


model{
 real mu[N];
 for(j in 1:NS){
   u[j] ~ multi_normal(zero, Sigma_u);
 }
 for(k in 1:NI){
   w[k] ~ multi_normal(zero, Sigma_w);
 }
 
 sigma_u ~ gamma(1.5, 1.0E-4);
 sigma_w ~ gamma(1.5, 1.0E-4);
 sigma_e ~ gamma(1.5, 1.0E-4);

 omega_u ~lkj_corr(2.0);
 omega_w ~lkj_corr(2.0);


 for(i in 1:N){
   mu[i] <- (beta[1] + u[subject[i], 1] + w[item[i], 1]) +
            (beta[2] + u[subject[i], 2] + w[item[i], 2])*cond[i];
 }

 rt ~ normal(mu, sigma_e);
}

When I ran the model, I got messages shown below. I wonder what may have caused it. Also, I have a 3-level categorical predictor. I wonder if it needs to be manually dummy coded. The session info is at the end of the post, and files containing the Stan model code, the dataset, and the R code to run the Stan code are attached to this post. Thank you!

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_u: validate transformed params N4stan5agrad3varE: Sigma_uSigma_u is not positive definite. Sigma_u(0,0) is 3.89667e-14:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_u: validate transformed params N4stan5agrad3varE: Sigma_uSigma_u is not positive definite. Sigma_u(0,0) is 1.41065e-14:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_w: validate transformed params N4stan5agrad3varE: Sigma_wSigma_w is not symmetric. Sigma_w[1,2] is -inf:0, but Sigma_w[2,1] element is -inf:0If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:stan::prob::lkj_corr_log(N4stan5agrad3varE): yy is not positive definite. y(0,0) is 1:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_w: validate transformed params N4stan5agrad3varE: Sigma_wSigma_w is not positive definite. Sigma_w(0,0) is 734.841:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_w: validate transformed params N4stan5agrad3varE: Sigma_wSigma_w is not positive definite. Sigma_w(0,0) is 5.00582:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Iteration:   1 / 1000 [  0%]  (Warmup)Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_u: validate transformed params N4stan5agrad3varE: Sigma_uSigma_u is not positive definite. Sigma_u(0,0) is 0:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_u: validate transformed params N4stan5agrad3varE: Sigma_uSigma_u is not positive definite. Sigma_u(0,0) is 9.96031e-88:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Invalid value of Sigma_w: validate transformed params N4stan5agrad3varE: Sigma_wSigma_w is not positive definite. Sigma_w(0,0) is 3.55041e+09:0.If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Iteration: 100 / 1000 [ 10%]  (Warmup)

Iteration: 200 / 1000 [ 20%]  (Warmup)

Iteration: 300 / 1000 [ 30%]  (Warmup)

Iteration: 400 / 1000 [ 40%]  (Warmup)

Iteration: 500 / 1000 [ 50%]  (Warmup)

Iteration: 501 / 1000 [ 50%]  (Sampling)

Iteration: 600 / 1000 [ 60%]  (Sampling)

Iteration: 700 / 1000 [ 70%]  (Sampling)

Iteration: 800 / 1000 [ 80%]  (Sampling)

Iteration: 900 / 1000 [ 90%]  (Sampling)

Iteration: 1000 / 1000 [100%]  (Sampling)

#  Elapsed Time: 20.6092 seconds (Warm-up)

#                13.1574 seconds (Sampling)

#                33.7665 seconds (Total)

sessionInfo()

R version 3.1.2 (2014-10-31)

Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:

[1] splines   stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

 [1] rstan_2.5.0       inline_0.3.13     psych_1.4.5       MCMCglmm_2.21     ape_3.1-4        

 [6] coda_0.16-1       lattice_0.20-29   multcomp_1.3-6    TH.data_1.0-3     survival_2.37-7  

[11] mvtnorm_1.0-0     lme4_1.1-7        Rcpp_0.11.3       Matrix_1.1-4      WRS_0.24         

[16] robustbase_0.91-1 akima_0.5-11      MASS_7.3-35      

loaded via a namespace (and not attached):

 [1] codetools_0.2-9  colorspace_1.2-4 corpcor_1.6.6    DEoptimR_1.0-1   digest_0.6.4    

 [6] ggplot2_1.0.0    grid_3.1.2       gss_2.1-3        gtable_0.1.2     minqa_1.2.3     

[11] munsell_0.4.2    nlme_3.1-118     nloptr_1.0.0     plyr_1.8.1       proto_0.3-10    

[16] reshape2_1.4     sandwich_2.3-2   scales_0.2.4     stats4_3.1.2     stringr_0.6.2   

[21] tensorA_0.36     zoo_1.7-11      

[Bob Carpenter]

What makes you think there's a problem? It's OK to get a few
of those messages during warmup. You can try to crank down the
initial step size lower (say 0.001), which should avoid some of
these messages during warmup. The problem is that the discretization
of the Hamiltonian needs to find the right step size before
things become arithmetically stable.

What did the output look like? Can you get it to fit simulated
data?

- 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.
> <model.stan.txt><dat.txt><Rcode.txt>

[Xiao He]

Thanks for the quick reply. The red font of the message made me a bit paranoid...that's why I thought there were problems. I printed the output, and based on the output, I think I made the mistake of coding the 3-level predictor as one variable (the 3 levels were coded as 1, 2, and 3). I thought Stan would automatically dummy code it as R would. This might explain why there were only two betas in my initial output, as shown below. I dummy coded the predictor, and re-ran the model. The new results (as shown by the output at the bottom of this reply) seem reasonable.  

Incorrect coding

> print(fit)
Inference for Stan model: model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.


              mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
beta[1]       5.78    0.00  0.05    5.69    5.75    5.78    5.81    5.87   632 1.00
beta[2]       0.06    0.00  0.01    0.04    0.06    0.06    0.07    0.08  4000 1.00

Dummy coded predictor

> print(fit)
Inference for Stan model: model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.


                mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
beta[1]         5.86    0.00  0.05    5.77    5.83    5.86    5.89    5.95   229 1.02
beta[2]         0.01    0.00  0.02   -0.03    0.00    0.01    0.03    0.06  1471 1.01
beta[3]         0.13    0.00  0.02    0.08    0.11    0.13    0.14    0.17  1268 1.00

[Shravan Vasishth]

Please also see our improved code (based on Stan 2.5 manual, so make sure you have Stan 2.5), which will be incorporated in the tutorial you linked to. This runs much faster (we can also fit 8x8 variance covariance matrices this way with very large datasets by psychology standards, and get sensible results): 

https://github.com/vasishth/StanJAGSexamples/blob/master/FrankEtAlCogSci2015/FTVCogSci.Stan