[R] Problem trying to use boot and lme together

From: Michael Dewey <med_at_aghmed.fsnet.co.uk>
Date: Wed 22 Jun 2005 - 01:49:59 EST


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@aghmed.fsnet.co.uk
http://www.aghmed.fsnet.co.uk/home.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 02:21:19 2005

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