[R-sig-ME] random effects CIs

[Greg Dropkin]

Hi Phillip

Thanks. Yes, I was aware of your points but what surprised me was that the
R2jags CI was closer to use.u = true.

Greg

> Hi Greg,
>
> I haven't checked your code for correctness, but the pattern you've
> found is typical. I don't know if there is a "best" answer, but broader
> confidence intervals are more conservative and so should have better
> coverage, if that's what you're going for.  The biggest difference for
> use.u=FALSE is how the random effects are handled, from the Details for
> ?bootMer:
>
>>  If ‘use.u’ is ‘FALSE’ and ‘type’ is ‘"parametric"’, each
>>  simulation generates new values of both the “_spherical_”
>>  random effects u and the i.i.d. errors epsilon, using
>>  ‘rnorm()’ with parameters corresponding to the fitted model
>>  ‘x’.
>>
>>  If ‘use.u’ is ‘TRUE’ and ‘type=="parametric"’, only the
>>  i.i.d. errors (or, for GLMMs, response values drawn from the
>>  appropriate distributions) are resampled, with the values of
>>  u staying fixed at their estimated values.
>
> In other words, use.u=FALSE generates new levels of the random effects,
> so in the case of Dyestuff, it simulates new, unseen Batches based on
> the estimates of the between-Batch variance. In the use.u=TRUE case, it
> re-uses the observed Batches and just simulates the observation-level /
> residual error.  So, the use.u=FALSE case introduces "new" variation /
> models additional variation and has broader confidence intervals. If I'm
> not mistaken, this is related to the difference between prediction and
> confidence intervals.
>
> Best,
> Phillip
>
>
>
> On 08/21/2017 10:21 PM, Greg Dropkin wrote:
>> 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
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
>