[R-sig-ME] problem bootstraping a mixed-model (lme)

[Ken Beath]

For parametric bootstrapping what is required is a model for the data, and
using that you generate the bootstrap samples, there is no resampling
involved.

For non-parametric bootstrapping with multilevel data, resampling
individual observations is not sufficient, what is required is to resample
whole clusters/subjects. This is possible in R by converting teh data to a
wide format, then in thefitting function taking the samples and expanding
them back to long and fitting the model. You will also need to create
unique id for each cluster before converting to long.

I recommend reading something like Manly's book on resampling and
bootstrapping, although I don't think it talks about multilevel models.


On 1 April 2015 at 23:30, David Villegas Ríos <chirleu at gmail.com> wrote:

> Dear all,
> I'm trying to boostrap repeatability estimated from a lme output.
> The model includes one fixed factor (month), one random factor (ID) and one
> correlation term to account for temporal autocorrelation of the replicates.
> I prefer parametric bootstraping since it is the recommended option
> according to Nakagawa and Schielzeth, 2010 (Biological reviews)
>
> These have been my attepmts so far:
>
> *Option 1: parametric bootstraping of the full model (what I really need)*
>
> bootcoef<-function(data, index){
>   dat<-data[index,]
>
>
> mod<-lme(dvm~factor(month),random=~1|ID,data=dat,correlation=corAR1(form=~month))
>
>
> return(as.numeric(VarCorr(mod))[1]/(as.numeric(VarCorr(mod))[1]+as.numeric(VarCorr(mod))[2]))
> # this is the repeatability estimate
> }
> output=boot(depm3,bootcoef,100,sim=parametric)
>
> *Error*: output$t yields 100 identical values.
>
> *Option 2: non-parametric bootstraping of the full model*
>
> bootcoef<-function(data, index){
>   dat<-data[index,]
>
>
> mod<-lme(dvm~factor(month),random=~1|ID,data=dat,correlation=corAR1(form=~month))
>
>
> return(as.numeric(VarCorr(mod))[1]/(as.numeric(VarCorr(mod))[1]+as.numeric(VarCorr(mod))[2]))
> # this is the repeatability estimate
> }
> output=boot(depm3,bootcoef,100)
>
> *Error*: Error in Initialize.corAR1(X[[2L]], ...) : covariate must have
> unique values within groups for "corAR1" objects
>
> *Option 3: parametric bootstraping of the model without the autocorrelation
> term*
>
> bootcoef<-function(data, index){
>   dat<-data[index,]
>   mod<-lme(dvm~factor(month),random=~1|ID,data=dat)
>
>
> return(as.numeric(VarCorr(mod))[1]/(as.numeric(VarCorr(mod))[1]+as.numeric(VarCorr(mod))[2]))
> # this is the repeatability estimate
> }
> output=boot(depm3,bootcoef,100,sim=parametric)
>
> *Erro*r: output$t yields 100 identical values which in addition I don't
> like because the autocorrelation term is not int he model
>
> *Option 4: non-parametric bootstraping of the model without the
> autocorrelation term*
>
> bootcoef<-function(data, index){
>   dat<-data[index,]
>   mod<-lme(dvm~factor(month),random=~1|ID,data=dat)
>
>
> return(as.numeric(VarCorr(mod))[1]/(as.numeric(VarCorr(mod))[1]+as.numeric(VarCorr(mod))[2]))
> # this is the repeatability estimate
> }
> output=boot(depm3,bootcoef,100)
>
> *Result*: I got 100 different values (this is ok), but I really need the
> autocorrelation term to be in.
>
> Is this something that you can comment about without reproducible data? Any
> suggestion would be greatly appreciated.
>
> David
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



-- 

*Ken Beath*
Lecturer
Statistics Department
MACQUARIE UNIVERSITY NSW 2109, Australia

Phone: +61 (0)2 9850 8516

Building E4A, room 526
http://stat.mq.edu.au/our_staff/staff_-_alphabetical/staff/beath,_ken/

CRICOS Provider No 00002J
This message is intended for the addressee named and may...{{dropped:9}}