R help - Jun 2005 - Problem trying to use boot and lme together

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.


Michael Dewey
med at aghmed.fsnet.co.uk
http://www.aghmed.fsnet.co.uk/home.html

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

The problem with simulate.lme is that it only returns logL for a given model
fitted to a simulated data set  - not the simulated data set itself (which one
might have expected a function with that name to do...). It would be nice with
that functionality...
S??ren
 

________________________________

Fra: r-help-bounces at stat.math.ethz.ch p?? vegne af Prof Brian Ripley
Sendt: ti 21-06-2005 18:53
Til: Michael Dewey
Cc: r-help-stat.math.ethz.ch
Emne: Re: [R] Problem trying to use boot and lme together



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

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html