How to find se.fit in R by hand?

Question

I'm trying to find by hand the se.fit. Focusing on the first observation (weight=$1.9$), when we write the following code:

install.packages("DAAG")
library("DAAG")

roller.lm <- lm(depression~weight,data=roller)
roller.pred <- predict(roller.lm,se.fit=T)
summary(roller.lm)[4]
roller.pred$fit[1]
roller.pred$se.fit[1]

We have this output

  > summary(roller.lm)[4]

  $coefficients
                Estimate  Std. Error   t value    Pr(>|t|)
  (Intercept)  -2.087148  4.7542813  -0.4390038  0.672274166
     weight     2.666746  0.7002426   3.8083171  0.005175013

   > roller.pred$fit[1]
   1 
   2.979669

   > roller.pred$se.fit[1]
   [1] 3.614297

So I want to use these informations to find the se.fit. Using the se of the coefficients, we have:

$$\hat y_1=$$

$=(-2.087148\pm 4.7542813)+(2.666746\pm0.7002426)\times1.9$

$=(-2.087148+2.666746\times1.9)\pm (4.7542813+ 0.7002426\times 1.9)$

$=2.9796694\pm 5.4545239$

So why R is giving me $3.614297$ instead of my calculation $5.4545239$? How can I discover the se.fit using these data?

Answer 1

You can’t just add up the standard errors. You need to work with the variances and covariances of the coefficient estimates, then take the square root to get back to the standard error scale.

The formula for the variance of a weighted sum of correlated variables is the key, as the estimates of the intercept and slope are correlated random variables, not independent. Use vcov(roller.lm) to get the variance/covariance matrix (variances are the diagonal elements, covariances off-diagonal). Then apply the formula to get the estimated variance of the weighted sum (weight 1 for the intercept, weight for the slope the value of x for which you want the estimate). The square root gives the standard error.

This answer works through the underlying code in detail, showing how you can then use the standard error to get confidence limits and prediction intervals.