[R] How does predict() calculate prediction intervals?

[Rolf Turner]

Just look at the code of predict.lm().  It is reasonably perspicuous.
In particular look at res.var.

     cheers,

         Rolf Turner

On 01/31/2013 05:50 AM, Kurt Rinehart wrote:
> For a given linear regression, I wish to find the 2-tailed t-dist
> probability that Y-hat <= newly observed values. I generate prediction
> intervals in predict() for plotting, but when I calculate my t-dist
> probabilities, they don't agree. I have researched the issues with variance
> of individual predictions and been advised to use the variance formula
> below (in the code).
>
> I presume my variance function differs from that used in predict(). Can
> someone advise me as to why I cannot reproduce the results of predict()?
>
> As a test, I calculated the prediction intervals about Y-hat by hand and
> compared them to predict():
>
> X <- log10(c(780,54,520,703,1356,1900,2500,741,1500,600))
> Y <-
> log10(c(0.037,0.036,0.026,0.0315,0.028,0.0099,0.00955,0.0405,0.019,0.0245))
> dat <- data.frame(cbind(X, Y))
> mod <- lm(Y ~ X, data = dat)
> rm(X,Y)
>
> ## model predictions, y.hat, at the new X values
> pred <- data.frame(predict.lm(mod,newdata = obs, interval = "prediction",
> level = 0.95))
>
> # 3 subsequent observations
> obs <- log10(data.frame(cbind(c(40, 200, 450), c(0.3, 0.06, .034))))
> names(obs) <- c("X", "Y")
>
> ## Calculating t-dist, 2-tailed, 95% prediction intervals for new
> observations
> mse <- anova(mod)[2,3]
> new.x <- obs$X - mean(dat$X)
> sum.x2 <- sum((dat$X - mean(dat$X))^2)
> y.hat <- pred$fit
> var.y.hat <- mse*(1+new.x^2/sum.x2)
>
> upr <- y.hat + qt(c(0.975), df = 8) * sqrt(var.y.hat)
> lwr <- y.hat + qt(c(0.025), df = 8) * sqrt(var.y.hat)
> hand <- data.frame(cbind(y.hat, lwr, upr))
>
> #The limits are not the same
> pred
> hand
>
> --
> Thank you,
> Kurt
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
>