[R] Problem trying to use boot and lme together

[Prof Brian Ripley]

We don't have the structure of your dataset.  But it seems pretty clear 
that your resampling is not preserving the random effects structure: you 
are always fitting to the same clusters labelled by gp and so missing the 
major source of variability.

You either need to resample clusters, or resample the cluster random 
effects, as well as resample the within-cluster effects.  I doubt if there 
is much theoretical support for such bootstrapping schemes, nor software 
support: I would prefer a so-called parametric bootstrapping approach, 
that is resample from the fitted model -- see simulate.lme.



On Tue, 21 Jun 2005, Michael Dewey wrote:

> The outcome variable in my dataset is cost. I prefer not to transform it as
> readers will really be interested in pounds, not log(pounds) or
> sqrt(pounds). I have fitted my model using lme and now wish to use boot on
> it. I have therefore plagiarised the example in the article in RNews 2/3
> December 2002
>
> after loading the dataset I source this file
>
> ====================================================
> require(boot)
> require(nlme)
> model.0 <- lme(tot ~ time + timepost + pct + totpat
>    + (time + timepost) : single + single
>    + (time + timepost) : train + train
>    + time:gpattend + timepost:gpattend + gpattend,
>    data = common,
>    random = ~time + timepost | gp
> )
> ints.9 <- intervals(model.0, which="fixed")$fixed[9,]
> #
> common$fit <- fitted(model.0)
> common$res <- resid(model.0, type = "response")
> cats.fit <- function(data) {
>    mod <- lme(tot ~ time + timepost + pct + totpat
>       + (time + timepost) : single + single
>       + (time + timepost) : train + train
>       + time:gpattend + timepost:gpattend + gpattend,
>    data = data,
>    random = ~ time + timepost | gp)
>    summ <- summary(mod)
>    c(fixef(summ), diag(summ$varFix))
> }
> model.fun <- function(d, i) {
>    d$tot <- d$fit+d$res[i]
>    cats.fit(d)
> }
> tot.boot <- boot(common, model.fun, R=999)
> ============================================
> So I fit the model and then generate fitted values and residuals which I
> use within the model.fun function to generate the bootstrap resample.
>
> If I do this the plot looks fine as does the jack.after.boot plot but the
> confidence intervals are about 10% of the width of the ones from the lme
> output. I wondered whether I was using the wrong fitted and residuals so I
> added level = 0 to the calls of fitted and resid respectively (level = 1 is
> the default) but this seems to lead to resamples to which lme cannot fit
> the model.
>
> Can anyone spot what I am doing wrong?
>
> I would appreciate a cc'ed response as my ISP has taken to bouncing the
> digest posts from R-help with probability approximately 0.3.
>
> FWIW I am using 2.1.0 under XP (SP2) with the versions of boot and nlme
> which shipped with the binary.


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595