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

Wed Apr 8 03:00:58 CEST 2015

Misinterpreted that. Anyway, the solution is then from first principles
simulate the required data sets.This may not be difficult, but it is not
something I have done for this type of model.

On 8 April 2015 at 10:33, Ben Bolker <bbolker at gmail.com> wrote:

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