confidence intervals

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: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch]On Behalf Of Spencer Graves
Sent: Friday, September 03, 2004 11:58 AM
To: bates at wisc.edu
Cc: Robert Waters; r-help at stat.math.ethz.ch
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}}