Last active

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist
View lmer vs JAGS comparison
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
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.