A political scientist wrote in with a question that actually comes up a lot, having to do with hard-to-estimate group-level variance and correlation parameters in multilevel models. The short answer is, when these things are hard to estimate, it’s often because there’s not much information about them in the data, and sometimes you can get away with just setting these parameters to zero and making your life easier.

Anyway, here’s his question:

I’m trying to estimate an overdispersed event count model with a correlated random intercept/random slope in BUGS. I’ve implemented both the multivariate normal and the scaled-inverse-Wishart approaches from 17.1 of your book in an overdispersed Poisson setup like your police stops model. The random intercept works fine, but the variance of the random slope is pretty low. As a result, BUGS gets stuck at zero all the time for the variance of the random slope even in the SIW which has the extra scale parameters, which means that the density of the variance (even after tons of iterations) has a big point mass at 0. (Also, estimated correlation between the random intercept and slope has an extremely wide, somewhat negatively skewed distribution (95% CI is often ~[-9, .9]).) Do I need to add further redundant parameterization per 19.5? I’m not sure how to do this. The only way I’ve been able to handle this issue is to put a uniform prior that does not include zero (i.e. sigma.b ~ dunif(.01, 20) on the random slope in the simpler model. This keeps the Gibbs sampler from getting stuck as much but might offend the vague priors people out there. Also, I’d prefer to use the SIW and have multiple varying coefficients if possible (that’s the next step in theory), which doesn’t allow me to manipulate the prior over rho directly.

My reply:

If it makes you feel any better, I struggle with implementing these varying-intercept, varying-slope models myself. Can you use lmer() with this model? I’m not sure if lmer() does overdispesed Poisson, but it might (it would be called “quasipoisson”, I think; I’m not sure. lmer() doesn’t always work either but at least it’s fast.

Regarding Bugs having convergence problems, yeah, a redundant multiplicative parameterization can help. Beyond this, if the variance of the slopes is low, you could simply fit a varying-intercept model. Remember that the group-level variance is the variance of the slopes, after accounting for any group-level predictors that have been interacted with your individual-level predictor that corresponds to the “slope”. So if the residual group-level variance is low, you don’t have to feel so bad about removing it from the model.

The other thing to remember about varying-intercept, varying-slopes is that your zero point should be reasonable. You could center your x variable or else have zero make some sort of sense, so that the “intercept” is interpretable.

The wide uncertainty on the correlation might be a computational problem or it might simply be a reflection of you having very little information on this parameter.

He then followed up:

I re-estimated it in lmer (quasi-poisson) and Stata’s xtmepoisson (not overdispersed) and both find a strong negative correlation between the random effects and slopes plus significant group-level variance, so I’m increasingly convinced that something is not kosher with my model. Also, running it as a negative binomial in WinBUGS with no correlation between the random effects yields very similar results to the lmer/xtempoisson. The predictor with the random slope is already centered to zero so I guess the only other step is redundant multiplicative parameterization, but adding that to the SIW seems really ugly — do you think it’s workable?

[Side note: There’s a kind of interesting interpretive question with random effects that have low group-level variance but affect the statistical significance of the predicted effect within groups. In my (misestimated) model, even though group-level variance is low, the random effects that are estimated change whether the 95% credible interval includes zero for specific groups.]

Also, before I wear out my welcome, do you know about any issues with using R2WinBUGS on 64 bit Windows boxes? We have high-end shared machines for the social sciences here and R2WinBUGS always times out on both when I call BUGS using the same data and script that works on my wife’s Windows box (“Error in while (regexpr(“\r”), info) > 0) { : argument is of length zero. Error in file(file,”r”): unable to open connection”). Besides some sort of permissions problem, 64 bit was the only other hypothesis the tech guy here had.

And I replied:

If the unexplained variance of the slopes is low, then it can happen that the correlation parameter is difficult to estimate. Basically, the likelihood on the correlation is really flat so the maximum could be anywhere–even at a place such as -0.98. This is consistent with your Bugs findings also.

I agree that all sorts of things can change whether the 95% interval includes zero. But ultimately that shouldn’t be a key decision point (even though I know it is).

I don’t know about R2WinBUGS on 64 bit Windows boxes.