[R-sig-ME] random effects CIs

[Greg Dropkin]

Dear list

apologies if these questions have already been answered here - if so just
point me to them.

Am I implementing these 4 approaches to CIs for the random effects
correctly (for Dyestuff), and if so, which is best (or is "best"
meaningless)? bootMer use.u=FALSE seems to be recommended in the lme4
manual, but gives quite different results to the other 3 which are roughly
comparable. Is there a better way?

thanks

Greg Dropkin

--

library(lme4)
library(lattice)
library(boot)
library(R2jags)

#1
#as in Ch 1 of Douglas Bates' book on lme4:
fm01<-lmer(Yield~1+(1|Batch),data=Dyestuff)
dotplot(ranef(fm01,condVar=TRUE))

attr(ranef(fm01,condVar=TRUE)$Batch,"postVar")
#all the same
sd1<-attr(ranef(fm01,condVar=TRUE)$Batch,"postVar")[[1]]^0.5
sd1
#the estimates are ranef(fm01), and the CIs shown in the dotplot are +/-
sd1*qnorm(0.975)
#e.g. for Batch A:
ranef(fm01)$Batch[1,]
ranef(fm01)$Batch[1,]+c(-1,1)*sd1*qnorm(0.975)

#2
#as recommended in the lme4 manual
#bootMer with use.u=FALSE
g<-function(.) ranef(.)$Batch[1,]
bootg<-bootMer(fm01,g,seed=1,nsim=1000,use.u=FALSE,.progress="txt",PBarg=list(style=3))
boot.ci(bootg,type=c("norm","basic","perc"))

#3
#bootMer with use.u=TRUE
bootgu<-bootMer(fm01,g,seed=1,nsim=1000,use.u=TRUE,.progress="txt",PBarg=list(style=3))
boot.ci(bootgu,type=c("norm","basic","perc"))

#4
#R2jags
I<-1:30
J<-1+floor((I-1)/5)

zY<-with(Dyestuff,(Yield-mean(Yield))/sd(Yield))
sdY<-with(Dyestuff,sd(Yield))

mod1.data<-list("zY","I","J")
mod1.parameters<-c("Con.adj","a.adj","mu.a","sigma.a","sig")

mod1.mod<-function()
{
Con.adj<-Con+mean(a[])
for (j in 1:6)
{
a[j]~dnorm(mu.a,1/sigma.a^2)
a.adj[j]<-a[j]-mean(a[])
}

for (i in I)
{
zY[i]~dnorm(lp[i],1/sig^2)
lp[i]<-Con+a[J[i]]
}
Con~dnorm(0,0.000001)
mu.a~dnorm(0,0.0001)
sigma.a~dunif(0,100)
sig~dunif(0,100)
}

set.seed(1)
mod1.mod.fit<-jags(data=mod1.data, inits=NULL,
parameters.to.save=mod1.parameters, n.chains=3, n.iter=50000,
n.burnin=10000, model.file=mod1.mod)
mod1.mod.fit
mat1<-mod1.mod.fit$BUGSoutput$sims.matrix

mean(mat1[,1]+mat1[,2])*sdY
quantile(mat1[,1]+mat1[,2],c(0.025,0.975))*sdY

#comparison

cim<-matrix(0,8,3)
rownames(cim)<-c("conditional sd","bootMer use.u=F norm","bootMer use.u=F
basic","bootMer use.u=F perc","bootMer use.u=T norm","bootMer use.u=T
basic","bootMer use.u=T perc","jags")
colnames(cim)<-c("estimate","low","high")
cim[1:7,1]<-ranef(fm01)$Batch[1,]
cim[1,2:3]<-ranef(fm01)$Batch[1,]+c(-1,1)*sd1*qnorm(0.975)
cim[2,2:3]<-boot.ci(bootg,type=c("norm","basic","perc"))$norm[,2:3]
cim[3,2:3]<-boot.ci(bootg,type=c("norm","basic","perc"))$basic[,4:5]
cim[4,2:3]<-boot.ci(bootg,type=c("norm","basic","perc"))$perc[,4:5]
cim[5,2:3]<-boot.ci(bootgu,type=c("norm","basic","perc"))$norm[,2:3]
cim[6,2:3]<-boot.ci(bootgu,type=c("norm","basic","perc"))$basic[,4:5]
cim[7,2:3]<-boot.ci(bootgu,type=c("norm","basic","perc"))$perc[,4:5]
cim[8,1]<-mean(mat1[,1]+mat1[,2])*sdY
cim[8,2:3]<-quantile(mat1[,1]+mat1[,2],c(0.025,0.975))*sdY

cim