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

Ken Beath ken.beath at mq.edu.au
Thu Apr 9 12:01:32 CEST 2015 

Profile likelihood isn't available with lme, only with lmer. Non-parametric
bootstrap should be possible. Here is some code for it, including 2 where
the model is misspecified. I think the code should work with unequal
numbers per subject but that would need to be checked, and the code should
be checked anyway. It will cope with errors, but replacing the results with
NA is not the best but only option I can think of.


library(nlme)
library(lme4)
library(boot)

fm.lme <- lme(Reaction ~ Days, data=sleepstudy, random=~1 | Subject)

summary(fm.lme)

sleepstudy.wide <-
reshape(sleepstudy,v.names="Reaction",idvar="Subject",timevar="Days",direction="wide")

oneboot <- function(d,i) {
    d <- d[i,]
    # create new Subject id so each is unique
    d$Subject <- 1:dim(d)[1]
    d <- reshape(d,direction="long")
    thecoef <- tryCatch({bootfm.lme <- lme(Reaction ~ Days, data=d,
random=~1 | Subject)
             fixef(bootfm.lme)},
      error = function(e) {print(e)
             return(rep(NA,length(fixef(fm.lme))))})
#    browser()
    return(thecoef)
}

sleep.boot <- boot(sleepstudy.wide, oneboot, R = 999)

# the standard errors will change as the model is misspecified as it is
missing
# the random effect for Days which gave some problems with lme and
bootstrapping
# due to lack of convergence in some samples
print(sleep.boot)
boot.ci(sleep.boot,index=2,type=c("norm","perc"))

fm2.lme <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)

summary(fm2.lme)

Orthodont.wide <-
reshape(Orthodont,v.names=c("distance","Sex"),idvar="Subject",timevar="age",direction="wide")

oneboot2 <- function(d,i) {
  d <- d[i,]
  # create new Subject id so each is unique
  d$Subject <- 1:dim(d)[1]
  d <- reshape(d,direction="long")
  thecoef <- tryCatch({bootfm2.lme <- lme(distance ~ age + Sex, data = d,
random = ~ 1|Subject)
                       fixef(bootfm2.lme)},
                      error = function(e) {print(e)
                      return(rep(NA,length(fixef(fm2.lme))))})
  #    browser()
  return(thecoef)
}

Orthodont.boot <- boot(Orthodont.wide, oneboot2, R = 999)

print(Orthodont.boot)
boot.ci(Orthodont.boot,index=2,type=c("norm","perc"))
boot.ci(Orthodont.boot,index=3,type=c("norm","perc"))

# fit as incorrect linear model
# then bootstrap to obtain correct standard errors
# if used for nonlinear models the coefficients will change as model will
# be population averaged
fm3.lm <- lm(distance ~ age + Sex, data = Orthodont)

summary(fm3.lm)

oneboot4 <- function(d,i) {
  d <- d[i,]
  # create new Subject id so each is unique
  d$Subject <- 1:dim(d)[1]
  d <- reshape(d,direction="long")
  thecoef <- tryCatch({bootfm3.lm <- lm(distance ~ age + Sex, data = d)
                       coef(bootfm3.lm)},
                      error = function(e) {print(e)
                       return(rep(NA,length(coef(fm3.lm))))})
  #    browser()
  return(thecoef)
}

Orthodont.boot3 <- boot(Orthodont.wide, oneboot4, R = 999)

print(Orthodont.boot3)
boot.ci(Orthodont.boot3,index=2,type=c("norm","perc"))
boot.ci(Orthodont.boot3,index=3,type=c("norm","perc"))


On 8 April 2015 at 16:55, David Villegas Ríos <chirleu at gmail.com> wrote:

> 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 may...{{dropped:9}}
>>
>> _______________________________________________
>> 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}}

More information about the R-sig-mixed-models mailing list