[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