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

[Ben Bolker]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 15-04-07 06:50 PM, Ken Beath wrote:
> 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.

  But the OP said that simulate() doesn't allow models with a
correlation structure (which doesn't really surprise me).
There's another problem, which is that simulate() for lme models does
*not*
work the same way that it does for most model type (the nlme package
is old enough that it had a simulate() method quite a while before
it was introduced as a more general method); it returns the
log-likelihoods
of null and full models rather than new (parametric-bootstrapped)
response vectors.

This might be hackable (especially to get simulate.lme() to return
vectors rather than log-likelihoods), but the internals of nlme are
not particularly simple & friendly ...

  cheers
    Ben Bolker

Reproducible example:

library("nlme")
fm1Dial <-
    lme(rate ~(pressure + I(pressure^2) + I(pressure^3) +
I(pressure^4))*QB,
        random= ~1|Subject, data= Dialyzer)
simulate(fm1Dial)
fm3Dial <- update(fm1Dial,
                 corr = corAR1(0.771, form = ~ 1 | Subject))
simulate(fm3Dial)
## Error in simulate.lme(fm3Dial) :
##   models with "corStruct" and/or "varFunc" objects not allowed

> 
> 
> 
> 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.
>>> 
>>> 
>> 
> 
> 

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (GNU/Linux)

iQEcBAEBAgAGBQJVJHdkAAoJEOCV5YRblxUHc7UIAMnTzI3edqGmcdKKOJsWqhB+
DlhDK3N+FnEtyncLJZqSfG473IJqp6UpujUDsbKFCnGdiUNHK8mmuCuDFvkFJCX/
XSquyTPZwjCYSmECmvsXep7C18PS/BNj+HuUQSqM7AT82T4Rr2tSJbw7nhhAyncI
xlQpmjCEQSAbmn67aGenhY1ac7B4faieJ8TCkkZ+j4zYEbtP5Ov7ofexGZ1acQYg
9F0H1E6JlVqDyy2iJXcM8AHbB5wU7xLg+fHATPNUE18fp+VTc8PvSBPBvbLIP0QS
EVDbmnzRyLe6zyjMiIgQx7XhayMyLKDCjHBMNfCpzQuWTxliIJTVPuiFLkDQTHo=
=Up8Y
-----END PGP SIGNATURE-----