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@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.
______________________________________________
R-help@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.
