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

[Søren Højsgaard]

Dear Jenni,

In the newest version of the doBy package there is a function called PBrefdist (and PBrefdist.mer) for calculating the reference distribution of the likelihood ratio statistic for comparing nested models. Looking into this function may help you. Perhaps the functions PBmodcomp and BCmodcomp (and their .mer methods) for calculating p-values based on parametric bootstrap can also be of inspiration.

Regards
Søren 

________________________________________
Fra: r-help-bounces at r-project.org [r-help-bounces at r-project.org] På vegne af Jenni Sanderson [Jennifer.Louise.Sanderson at googlemail.com]
Sendt: 8. juni 2011 09:56
Til: r-help at r-project.org; r-sig-mixed-models at r-project.org
Emne: [R] using stimulate(model) for parametric bootstrapping in lmer   repeatabilities

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

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.