lmer vs JAGS comparison

library(MASS)
new.df <- function(cond1.rt=600, effect.size=10, sdev=40,
                   sdev.int.subj=10, sdev.slp.subj=10, 
                   rho.u=0.6,
                   nsubj=10,
                   sdev.int.items=10, sdev.slp.items=10, 
                   rho.w=0.6,
                   nitems=10) {

  ncond <- 2

  subj <- rep(1:nsubj, each=nitems*ncond)
  item <- rep(1:nitems, nsubj, each=ncond)
  ## needs to be generalized:
  cond <- rep(0:1, nsubj*nitems)
  err  <- rnorm(nsubj*nitems*ncond, 0, sdev)
  d    <- data.frame(subj=subj, item=item, cond=cond+1, err=err)

  Sigma.u<-matrix(c(sdev.int.subj^2,
                    rho.u*sdev.int.subj*sdev.slp.subj,
                    rho.u*sdev.int.subj*sdev.slp.subj,
                    sdev.slp.subj^2),nrow=2)

  Sigma.w<-matrix(c(sdev.int.items^2,
                    rho.u*sdev.int.items*sdev.slp.items,
                    rho.u*sdev.int.items*sdev.slp.items,
                    sdev.slp.items^2),nrow=2)

  
  
  # Adding random intercepts and slopes for subjects:
  ## first col. has adjustment for intercept, 
  ## second col. has adjustment for slope
  subj.rand.effs<-mvrnorm(n=nsubj,rep(0,ncond),Sigma.u)
  
  item.rand.effs<-mvrnorm(n=nitems,rep(0,ncond),Sigma.w)
  
   re.int.subj <- subj.rand.effs[,1]
  d$re.int.subj <- rep(re.int.subj, each=nitems*ncond)
  re.slp.subj   <- subj.rand.effs[,2]

  ## sanity check:
  print(cor(re.int.subj,re.slp.subj))
  d$re.slp.subj <- rep(re.slp.subj, each=nitems*ncond) * (cond - 0.5)

  # Adding random intercepts and slopes for items:
  re.int.item <- item.rand.effs[,1]
  d$re.int.item <- rep(re.int.item, nsubj, each=ncond)
  re.slp.item <- item.rand.effs[,2]
  d$re.slp.item <- rep(re.slp.item, nsubj, each=ncond) * (cond - 0.5)

  ## sanity check:
  print(cor(re.int.item,re.slp.item))


  d$rt <- (cond1.rt + cond*effect.size
           + d$re.int.subj + d$re.slp.subj
           + d$re.int.item + d$re.slp.item
           + d$err)

  d
}


d<-new.df(nsubj=10,nitems=10,rho.u=0.6,rho.w=0.6)
colnames(d)

d<-d[,c(1,2,3,9)]
d$c1<-ifelse(d$cond==1,-0.5,0.5)

library(lme4)
m1<-lmer(rt~c1+(1+c1|subj)+(1+c1|item),d)
summary(m1)

cat("
    # Fixing data to be used in model definition
    data
{
    zero.u[1] <- 0
    zero.u[2] <- 0
    zero.w[1] <- 0
    zero.w[2] <- 0
}
# Then define model
    model
{
# Intercept and slope for each person, including random effects
    for( i in 1:I )
{
    u[i,1:2] ~ dmnorm(zero.u,invSigma.u)
}
# Intercepts and slope by item
    for( k in 1:K)
{
    w[k,1:2] ~ dmnorm(zero.w,invSigma.w)
}
    
# Define model for each observational unit
    for( j in 1:N )
{
    mu[j] <- ( beta[1] + u[subj[j],1] + w[item[j],1]) +
    ( beta[2] + u[subj[j],2] + w[item[j],2] ) * ( c1[j] )
    rt[j] ~ dnorm( mu[j], tau.e )
    pred[j] ~ dnorm( mu[j], tau.e )
}
    
    minimum <- min(pred)
    maximum <- max(pred)
    mean <- mean(pred)
        
# Fixed intercept and slope (uninformative)
    beta[1] ~ dnorm(0.0,1.0E-5)
    beta[2] ~ dnorm(0.0,1.0E-5)
    
# Residual variance
    tau.e <- pow(sigma.e,-2)
    sigma.e  ~ dunif(0,15)

# SUBJECT:    
# variance-covariance matrix of subject ranefs
    invSigma.u  ~ dwish( R.u , 2 )
    R.u[1,1] <- pow(sigma.a,2)
    R.u[2,2] <- pow(sigma.b,2)
    R.u[1,2] <- rho.u*sigma.a*sigma.b
    R.u[2,1] <- R.u[1,2]
    Sigma.u <- inverse(invSigma.u)
    
# var intercepts, var slopes, following Chung et al 2013 recommendations.
## Not used because this systematically underestimates the variance components
#    tau.a ~  dgamma(1.5, pow(10,-4))   
#    tau.b ~ dgamma(1.5, pow(10,-4))
#    sigma.a <- pow(tau.a,-1/2)
#    sigma.b <- pow(tau.b,-1/2)

# gives better estimates:
     sigma.a ~ dunif(0,15)
     sigma.b ~ dunif(0,15)
     tau.a <- 1/pow(sigma.a,2)
     tau.b <- 1/pow(sigma.b,2)


# correlation, following Chung et al 2013 recommendations.
#    rho.u <- rho2.u*2-1
#    rho2.u ~ dbeta(1.5,1.5)
     rho.u ~ dunif(-1,1)   

# ITEMS:
# variance-covariance matrix of item ranefs
    invSigma.w ~ dwish( R.w, 2 )
    R.w[1,1] <- pow(sigma.c,2)
    R.w[2,2] <- pow(sigma.d,2)
    R.w[1,2] <- rho.w*sigma.c*sigma.d
    R.w[2,1] <- R.w[1,2]
    Sigma.w <- inverse(invSigma.w)
    
# var intercepts, var slopes
#    tau.c ~ dgamma(1.5, pow(10,-4))
#    tau.d ~ dgamma(1.5, pow(10,-4))
#    sigma.c <- pow(tau.c,-1/2)
#    sigma.d <- pow(tau.d,-1/2)
# gives better estimates:
     sigma.c ~ dunif(0,15)
     sigma.d ~ dunif(0,15)
     tau.c <- 1/pow(sigma.c,2)
     tau.d <- 1/pow(sigma.d,2)

# correlation
#    rho.w <- rho2.w*2-1
#    rho2.w ~ dbeta(1.5,1.5)
     rho.w ~ dunif(-1,1)   
}",file="simulateddata.jag")

dat <- list( subj = sort(as.integer(factor(d$subj) )),
                      item = sort(as.integer(factor(d$item) )),
                      rt = d$rt,
                      c1 = d$c1,
                      N = nrow(d),
                      I = length( unique(d$subj) ),
                      K = length( unique(d$item) ))  


library(rjags)

sim.mod <- jags.model(
  file = "simulateddata.jag",
  data = dat,
  n.chains = 4,
  n.adapt = 20000 ,quiet=T) 

track.variables<-c("beta","sigma.e","sigma.a","sigma.b","sigma.c","sigma.d","rho.u","rho.w")
update(sim.mod,n.iter=40000)
sim.res <- coda.samples(sim.mod,
                            var = track.variables,
                            n.iter = 80000,
                            thin = 20 )

gelman.diag(sim.res) 

plot(sim.res)

res<-round(cbind(summary(sim.res)$statistics[,1:2],
      summary(sim.res)$quantiles[,c(1,3,5)]),digits=3)

## compare results:
res
           Mean    SD    2.5%     50%   97.5%
beta[1] 605.956 7.143 591.435 606.035 620.039
beta[2]  14.643 5.834   2.878  14.749  26.125
rho.u     0.088 0.549  -0.910   0.124   0.953  <- underestimate
rho.w     0.070 0.556  -0.923   0.100   0.953  <- underestimate
sigma.a   8.835 4.069   0.748   9.353  14.734
sigma.b   7.888 4.190   0.580   8.002  14.603
sigma.c   8.584 4.150   0.705   9.023  14.734
sigma.d   8.035 4.218   0.473   8.253  14.658
sigma.e  14.989 0.011  14.959  14.992  15.000

summary(m1)
Random effects:
 Groups   Name        Variance Std.Dev. Corr
 subj     (Intercept) 3.27e+02 18.080       
          c1          3.50e+01  5.913   1.00 <- overestimate (failure to estimate)
 item     (Intercept) 1.08e+02 10.380        
          c1          8.12e-02  0.285   1.00 <- overestimate (failure to estimate)
 Residual             1.77e+03 42.032       
Number of obs: 200, groups: subj, 10; item, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)   606.41       7.23    83.9
c1             14.75       6.23     2.4