Computing by hand vs. R's magic wand

Question

I have fitted a polynomial to a data set that I have. Thus I have obtained coefficients $\beta_i$ for $i=0,1,2$ and have a relationship of the form is $$Y=\beta_0+\beta_{1}X+\beta_{2}X^{2}+\varepsilon$$ between the predictor $X$ and the response $Y$, where $\varepsilon$ is the error term.

My question is about computing confidence intervals.
I can compute for each of the coefficients $\beta_i$ the confidence interval using the confint()function from R. This gets me three (95%) confidence intervals $[a_0,b_0]$,$[a_1,b_1]$ and $[a_2,b_2]$ for $\beta_0$,$\beta_1$ and $\beta_2$.

Lets say I would like to predict the average value that I will obtain for $Y$ for a new sample point $x=0.25$, that is, I would like to know the confidence interval around this point.
Using the previously obtained confidence intervals for coefficients, I can compute that $x$ has to lie in (with probability 95%) in the interval $$[a_0+a_{1}0.25+a_{2}0.25^{2},b_0+b_{1}\tilde 0.25+b_{2}0.25^{2}].$$

Unfortunately, when I compare the values that I obtain in this way with the output of the R command predict(my_fitted_data, data.frame(x=0.25),interval = 'confidence'), that interval is way smaller than the one I obtaind by computing the interval, as described, by hand.
Did I do anything wrong? If so, what?

Thank you.

Answer 1

The mistake you made was to ignore the covariances between the estimators $\hat{\beta}_0$, $\hat{\beta}_1$ and $\hat{\beta}_2$.

You want a 95% confidence interval for $$ E(Y|X=0.25) = \beta_0 + 0.25\beta_1 + 0.25^2\beta_2 = a\beta, $$ with $a = [1\ \ 0.25\ \ 0.25^2]$ and $\beta = [\beta_0\ \ \beta_1\ \ \beta_2]^T$.

Under the usual assumptions (normally distributed errors, etc.), the confidence interval you want is $$ a\hat{\beta} \pm t_{n-p,\, 0.975} \sqrt{a\hat{\sigma}^2 (\mathbf{X}^T\mathbf{X})^{-1}a^T}, $$ with $n$ being the sample size, $p$ being the number of regression coefficients including the intercept (i.e. $p=3$), and $t_{n-p,\, 0.975}$ being the value such that $P(T_{n-p} \leq t_{n-p,\, 0.975}) = 0.975$ if $T_{n-p}$ is a Student t random variable with $n-p$ degrees of freedom. We use 0.975 because that gives 0.025 in the upper tail and 0.025 in the lower tail. The matrix $\hat{\sigma}^2 (\mathbf{X}^T\mathbf{X})^{-1}$ is the estimated covariance matrix of the estimators $\hat{\beta}$, which can be obtained in R using the vcov function ($\mathbf{X}$ is the design matrix).

I acknowledge NRH, who showed me this when I asked a similar question (the notation I've used here is slightly different).

An example:

set.seed(1)

d <- data.frame(x = rnorm(50), y = rnorm(50))

lm1 <- lm(y ~ x + I(x^2), data=d)

x.value <- 0.25
conf.level <- 0.95

predict(lm1, data.frame(x = x.value), interval="confidence", level = conf.level)

##         fit        lwr       upr
## 1 0.1254577 -0.2068385 0.4577539

a <- c(1, x.value, x.value^2)

# percentile of the t-distribution to use (0.975 if conf.level is 0.95)
perc <- 1 - (1 - conf.level)/2

# degrees of freedom = sample size - number of coefs
DF <- nrow(d) - length(coef(lm1))

# Confidence interval
(a %*% coef(lm1)) + c(-1, 1)*qt(perc, DF)*sqrt(a %*% vcov(lm1) %*% a)

## [1] -0.2068385  0.4577539