Re: [R] Problem trying to use boot and lme together

From: Søren Højsgaard <Soren.Hojsgaard_at_agrsci.dk>
Date: Wed 22 Jun 2005 - 06:28:42 EST


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

______________________________________________
R-help@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
Received on Wed Jun 22 06:37:34 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:56 EST