Crazy divergence in mixed model - chains start well but go crazy later

Question

I'm trying to compute a mixed model using jags (R2jags) and I got very strange divergence. The chains started so well, very well according to the results expected (also see lmer output of the same model below). But at certain point, the chains went crazy. Just look at the traceplot for delta_tau variable - the chains start so well, but then the green chain goes crazy, then blue and finally red...

Any ideas why this happens? Can't be in initial values, because the chains started so well. Maybe the priors? Why is the system unstable?

enter image description here

EDIT: Variables gamma_tau and delta_tau don't fall to exact zero, as you can see on these zoomed-in figures:

enter image description here

This is the jags model:

model {

# likelihood
for (i in 1:N) {
    logInd[i] ~ dnorm(mu[i], eps_tau)
    mu[i] <- alpha[crit[i]] + (beta[crit[i]] + delta[species[i]])*year[i] + gamma[species[i]] # ekviv mix1b/c podle me
}

# priors
eps_tau ~ dgamma(1.0E-3, 1.0E-3) 

for (j in 1:no_crit) {
    alpha[j] ~ dnorm(0, 0.0001)
    beta[j] ~ dnorm(0, 0.0001) 
}

for (k in 1:no_species) {
    gamma[k] ~ dnorm(0, gamma_tau)
    delta[k] ~ dnorm(0, delta_tau)
}

gamma_tau ~ dgamma(1.0E-3, 1.0E-3) 
delta_tau ~ dgamma(1.0E-3, 1.0E-3)
}

Code used to run jags (using R2jags):

no_crit = length(levels(crit))

win.data = list(logInd = mydata$logInd, crit = (as.integer(crit)), 
 year = mydata$Year, species = (as.integer(mydata$Taxon)),
 N = nrow(mydata), no_crit = no_crit, no_species = length(levels(mydata$Taxon))
)

inits = function () { list(
    alpha = rnorm(no_crit, 0, 10000),
    beta = rnorm(no_crit, 0, 10000)
)}  

params = c("alpha", "beta", "eps_tau", "gamma_tau", "delta_tau")

# ni: 1000 -> .. sec
ni <- 20000
nt <- 8
nb <- 8000
nc <- 3

out <- R2jags::jags(win.data, inits, params, "model.txt",
    nc, ni, nb, nt,  
    working.directory = paste(getwd(), "/tmp_bugs/", sep = "")
)
R2jags::traceplot(out, mfrow = c(4, 2))

Here is output from the equivalent lmer model:

> summary(lmer(logInd ~ 0 + crit_i + Year:crit_i + (1 + Year|Taxon), data = datai2))
Linear mixed model fit by REML 
Formula: logInd ~ 0 + crit_i + Year:crit_i + (1 + Year | Taxon) 
   Data: datai2 
  AIC  BIC logLik deviance REMLdev
 8558 8630  -4267     8495    8534
Random effects:
 Groups   Name        Variance   Std.Dev.   Corr  
 Taxon    (Intercept) 1.1682e-12 1.0808e-06       
          Year        5.3860e-07 7.3389e-04 0.000 
 Residual             8.7038e-01 9.3294e-01       
Number of obs: 2987, groups: Taxon, 103

Fixed effects:
               Estimate Std. Error t value
crit_iA      29.0539403  8.8116915   3.297
crit_iF       0.1848404  6.0286726   0.031
crit_iU      12.3405800 10.3326242   1.194
crit_iW       5.3248537  9.7416915   0.547
crit_iA:Year -0.0122717  0.0044174  -2.778
crit_iF:Year  0.0022365  0.0030222   0.740
crit_iU:Year -0.0038701  0.0051799  -0.747
crit_iW:Year -0.0003054  0.0048836  -0.063

Correlation of Fixed Effects:
            crit_A crit_F crit_U crit_W cr_A:Y cr_F:Y cr_U:Y
crit_iF      0.000                                          
crit_iU      0.000  0.000                                   
crit_iW      0.000  0.000  0.000                            
crit_iA:Yer -0.999  0.000  0.000  0.000                     
crit_iF:Yer  0.000 -0.999  0.000  0.000  0.000              
crit_iU:Yer  0.000  0.000 -0.999  0.000  0.000  0.000       
crit_iW:Yer  0.000  0.000  0.000 -0.999  0.000  0.000  0.000

Thanks in advance!

Cross posting is typically discouraged, http://stackoverflow.com/q/12714591/604456 — Andy W, Oct 03 '12 at 19:01
Here's a recent discussion regarding cross-posting of equivalent questions on different stack exchange sites: http://meta.stackexchange.com/questions/149257. Your question seems more appropriate to this site than to SO to me. — smillig, Oct 03 '12 at 19:28
It might be that `lmer` doesn't show correct results and the chains don't turn crazy, but turn right! Considering also the [suspicious differences between lmer and lme...](http://stats.stackexchange.com/q/38524/5509) — Tomas, Oct 03 '12 at 21:32
I didn't have time to read carefully, but my quick suggestion: fit a model from simulated data where you know the true parameters (take the parameters from lmer, for instance). Also, try a uniform prior on the variance, rather than a gamma. — Manoel Galdino, Oct 07 '12 at 16:44

tedddd · Answer 1 · 2012-10-03T20:40:50.933

2

I can't speak to this specific class of models, but it appears you sampler has found another (probably higher) mode. At this mode the variances of gamma (this is the species random effect?) and delta, are 0. When this effect isn't present, there is an identification problem, which is what your correlations are telling you.

Try adjusting the hyperparameters for your gamma_tau and delta_tau distributions to make them tighter away from zero.

edited Oct 03 '12 at 20:40

answered Oct 03 '12 at 20:07

tedddd

21
2

Hi Ted, thanks for answer. Few questions on your answer: 1) what you mean with "probably higher mode"? What do you mean with higher? More probable? 2) gamma is a species random effect for intercept. What you probably meant is that for the second mode, gamma_tau goes to lower values? 3) what do you mean with make it thighter, shall I use `delta_tau ~ dgamma(9.01, 0.01)` instead of `dgamma(1.0E-3, 1.0E-3)`? 4) so you basically suggest to adjust priors in favor of the second mode? – Tomas Oct 03 '12 at 20:24
1

1) The posterior in this region is much higher than the posterior evaluated at the earlier parts of your chain. 2) From your traceplot, gamma_tau and and delta_tau are zero, which means there are no species random effects, which is causing an identification problem. 3) I'm not an expert in JAGs syntax, but I think that this would tighter your prior, but also move it closer to zero. Try a mean preserving tightening. – tedddd Oct 03 '12 at 20:36
1) posterior but what? I understand "posterior" as adjective 2) they are not exact zero, please see updated question. – Tomas Oct 03 '12 at 21:29
So is 4) what you suggest? – Tomas Oct 03 '12 at 21:29
Try fixing those two parameters. – tedddd Oct 05 '12 at 20:10

Crazy divergence in mixed model - chains start well but go crazy later

1 Answers1