Skip to content

@vasishth /lmer vs JAGS comparison
Last active

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.