problem bootstraping a mixed-model (lme)

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

--

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