Covariance Matrix not Fitting Correctly in Varying Slopes Model

[mrdevlar]

Covariance Matrix not Fitting Correctly in Varying Slopes Model

[rmcelreath]

Seems like Stan changed something about how warmup and/or sampling works. If you set control=list(adapt_delta=0.99), the spike goes away with modest samples. So that suggests a change to warmup algorithm.

I don't understand why the results should differ, however, once the adaptation is adjusted.

I can't get the same random data from the provided seed, in any event. I can't imagine anything changed with the way the seed is used. So I must have used a different seed for the example shown in 13.4 and 13.5. I'll see if I figure it out.

[mrdevlar]

Thanks the adaptive_delta seems to have eliminated the spike. I initially expected that the mass at 0 was shifting the distribution's MAP, but it seems with the adaptive_delta it is still at:

               Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
Rho[1,1]       1.00   0.00       1.00       1.00  6000  NaN
Rho[1,2]      -0.18   0.33      -0.71       0.33  2211 1.00
Rho[2,1]      -0.18   0.33      -0.71       0.33  2211 1.00
Rho[2,2]       1.00   0.00       1.00       1.00  6000 1.00

Rather bizarre. Just to make sure I am understanding this correctly, I should expect Rho[1,2], Rho[2,1] to be -0.7 yes? As the original Rho matrix we defined earlier:

sigmas <- c(sigma_a,sigma_b) # standard deviations
Rho <- matrix( c(1,rho,rho,1) , nrow=2 ) # correlation matrix

Returns:

     [,1] [,2]
[1,]  1.0 -0.7
[2,] -0.7  1.0

Right? Or am I missing something?

[rmcelreath]

I think the provided seed in the text won't reproduce Figure 13.4. That's an issue. I should figure out which random seed will reproduce that figure.

But we shouldn't expect any particular simulation to recover the "true" correlation of -0.7. Many simulated sets of data will produce a milder empirical correlation than that, and then the finite evidence will lead to regularization towards a less strong correlation.

[paul-buerkner]

I think the Stan team recently changed the default of adapt_delta from 0.95 to 0.8, which may explain the convergence problems. Apart from that, adapt_delta = 0.8 should work ok when using the non-centered parameterizations.

When trying to get reliable estimates of random effects correlations, n = 20 (for the cafes) will likely not suffice. Also lkjcorr(2) pulls the estimates somewhat towards zero. As long as rho = -0.7 is in the 95% credible interval (when using Rho ~ lkjcorr(1)) I wouldn't worry too much.

[lexmart]

set.seed(4) instead of set.seed(5) gives me roughly the same plot as figure 13.4 (1st edition)