How does predict() calculate prediction intervals?

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