[R] Problem trying to use boot and lme together

Michael Dewey med at aghmed.fsnet.co.uk
Tue Jun 21 17:49:59 CEST 2005


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




More information about the R-help mailing list