Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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 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 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!

share|improve this question
2  
Cross posting is typically discouraged, stackoverflow.com/q/12714591/604456 –  Andy W Oct 3 '12 at 19:01
    
This question unfortunatelly fits both sites equally.. –  Curious Oct 3 '12 at 19:25
1  
Here's a recent discussion regarding cross-posting of equivalent questions on different stack exchange sites: meta.stackexchange.com/questions/149257. Your question seems more appropriate to this site than to SO to me. –  smillig Oct 3 '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... –  Curious Oct 3 '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 7 '12 at 16:44
add comment

1 Answer

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.

share|improve this answer
    
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? –  Curious Oct 3 '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 3 '12 at 20:36
    
1) posterior but what? I understand "posterior" as adjective 2) they are not exact zero, please see updated question. –  Curious Oct 3 '12 at 21:29
    
So is 4) what you suggest? –  Curious Oct 3 '12 at 21:29
    
Try fixing those two parameters. –  tedddd Oct 5 '12 at 20:10
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.