[R] using stimulate(model) for parametric bootstrapping in lmer repeatabilities

Jenni Sanderson Jennifer.Louise.Sanderson at googlemail.com
Wed Jun 8 09:56:04 CEST 2011


Hi all,

I am currently doing a consistency analysis using an lmer model and
trying to use parametric bootstrapping for the confidence intervals.

My model is like this:

model<-lmer(y~A+B+(1|C/D)+(1|E),binomial)

where E is the individual level for consistency analysis, A-D are
other fixed and random effects that I have to control for.

Following Nakagawa and Scheilzeth I can work out the repeatability
estimate using the following (as it is a binomial and the residual
variance is fixed at 1).

attr(lme4::VarCorr(model)$E, "stddev")^2 /
(1*(pi^2)/3 + attr(lme4::VarCorr(model)$E, "stddev")^2  +
attr(lme4::VarCorr(model)$C, "stddev")^2  +
attr(lme4::VarCorr(model)$D, "stddev")^2  )

My question is can I use stimulate(model) to generate values that I
can then use to do parametric bootstrap analysis and generate the
confidence intervals?

Something like this:

n<-length(A)
niter<-1000
y<-matrix(nrow=n,ncol=niter*2)
for (i in 1:niter) {
y[,I(i*2-1):I(i*2)]<-simulate(model)[,1]   }
rvalues<-numeric()
for (i in 1:niter) {
yboot<-cbind(y[,I(i*2-1)],y[,I(i*2)])
mboot<-lmer(y~A+B+(1|C/D)+(1|E),binomial)
rvalues[i]<- attr(lme4::VarCorr(mboot)$E, "stddev")^2/(1*(pi^2)/3 +
attr(lme4::VarCorr(mboot)$E, "stddev")^2  +
attr(lme4::VarCorr(mboot)$C, "stddev")^2  +
attr(lme4::VarCorr(mboot)$D, "stddev")^2  )}
confidence.intervals<-quantile(rvalues,c(0.05,0.95))

In the guide to lme4 it says that stimulate() "generate simulations
based on the estimated fitted models (conditional on the estimated
values of both the random and fixed effects)", which sounds like
exactly what I would need to generate values for parametric
bootstrapping but I can't find any examples of where anyone has done
this.

Any advice would be very much appreciated! Thank you very much.

Jenni

Jenni Sanderson
PhD student - Conflict and Cooperation in Vertebrate Societies
University of Exeter



More information about the R-help mailing list