[R] confidence intervals

[Rogers, James A [PGRD Groton]]

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}}