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