[R-sig-ME] problem bootstraping a mixed-model (lme)
David Villegas Ríos
chirleu at gmail.com
Wed Apr 8 08:55:11 CEST 2015
Hi Ken and Ben.
Unfortunately this seems to be beyond my statistical/coding skills.
Maybe the solution is to look for alternatives to get a CI for the
repeatability for my model, like likelihood profiling as suggested by the
original paper of Nakagawa and Schielzeth, 2010. Does it sound good? Any
better alternative?
I will give it a try...
Thanks for your comments.
David
2015-04-08 3:00 GMT+02:00 Ken Beath <ken.beath at mq.edu.au>:
> 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 m...{{dropped:10}}
More information about the R-sig-mixed-models
mailing list