Robert, I have this quick hack to obtain approximate "Shewhart" prediction intervals for variance component models fit with lme (to nitpick slightly, "confidence" intervals have the interpretation of containing parameters, while "prediction" and "tolerance" intervals have the interpretation of containing future observations or statistics).
Back of the envelope documentation: the only argument that probably needs explaining is the "reps" argument to shewhart(). If your model is, e.g. fixed = Y ~ 1, random = ~ 1 | Batch then specify reps = c(1, 1) if you want to predict a single future observation from a single future batch, reps = c(1, 2) if you want to predict the mean of two future observations from a single future batch, reps = c(2, 2) if you want to predict the mean of 4 observations spread evenly over 2 future batches, ... Leave mult.check = 1, unless you want to do a Bonferroni correction. HTH, Jim Rogers valStats2 <- function (x, fixed, random, ...) { mod <- lme(fixed = fixed, data = x, random = random, ...) mn <- fixef(mod) vc <- VarCorr(mod) err <- "Expecting only random intercept terms and a single fixed intercept.\n" if (length(mn) > 1 || ncol(vc) > 2) stop(err) rn <- rownames(vc) skip <- regexpr("=", rn) > 0 if (!any(skip)) vnms <- attr(vc, "title") else vnms <- grep("=", rn, value = TRUE) vc <- vc[!skip, ] vnms <- trim(sub("=.*", "", vnms)) vnms <- c(vnms, "Residual") vnms <- paste("V", vnms, sep = ".") vars <- as.numeric(vc[, "Variance"]) stats <- c(mn, vars) names(stats) <- c("Intercept", vnms) stats } shewhart <- function (x, meancol = "Intercept", varcols = grep("^V\\.", names(x), value = TRUE), reps = c(1, 1), alpha = 0.02, mult.check = 1) { mn <- x[[meancol]] vr <- as.matrix(x[varcols]) totvar <- vr %*% (1/reps) totsd <- sqrt(totvar) LL.mean <- mn + qnorm(alpha/2/mult.check) * totsd UL.mean <- mn + qnorm(1 - alpha/2/mult.check) * totsd out <- data.frame(V.Total = totvar, LL.mean = LL.mean, UL.mean = UL.mean) out } ### Example, where x is your data.frame: foo <- valStats2(x, fixed = Y ~ 1, random = ~ 1|Batch) foo <- as.data.frame(t(as.matrix(foo))) data.frame(foo, shewhart(foo)) > -----Original Message----- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] Behalf Of Spencer Graves > Sent: Friday, September 03, 2004 11:58 AM > To: [EMAIL PROTECTED] > Cc: Robert Waters; [EMAIL PROTECTED] > Subject: Re: [R] confidence intervals > > > Hi, Robert: > > While it may be difficult to program this in general > (as suggested > by it's position on Doug's "To Do" list), all the pieces should be > available to support a special script for your specific application. > What fixed and random model(s) interest you most? > > hope this helps. spencer graves > > Douglas Bates wrote: > > > Robert Waters wrote: > > > >> Dear R users; > >> > >> Im working with lme and Id like to have an idea of how > >> can I get CI for the predictions made with the model. > >> Im not a stats guy but, if Im not wrong, the CIs > >> should be different if Im predicting a new data point > >> or a new group. Ive been searching through the web and > >> in help-lists with no luck. I know this topic had been > >> asked before but without replies. Can anyone give an > >> idea of where can I found information about this or > >> how can I get it from R? > >> > >> Thanks for any hint > > > > > > That's not currently implemented in lme. It's on the "To > Do" list but > > it is not very close to the top. > > LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}} ______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html