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

[Ken Beath]

There is a simulate function for lme, look for it under the help. Then you
should be able just to modify the code on page 182 of Faraway.



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

> Thanks Ken.
> I tried to follow Faraday's book (Extending the linear model with R)
> guidelines for parametric bootstraping but simulate() does not allow models
> with a correlation structure, which is my case.
> Is there any way to circumvent this issue?
>
> David
>
> 2015-04-02 0:03 GMT+02:00 Ken Beath <ken.beath at mq.edu.au>:
>
>> 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 contain
>> confidential information.  If you are not the intended recipient, please
>> delete it and notify the sender.  Views expressed in this message are those
>> of the individual sender, and are not necessarily the views of the Faculty
>> of Science, Department of Statistics or Macquarie University.
>>
>>
>


-- 

*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}}