An interval estimate for the mean of the response for this set of predictors $x^\star = (1, x_i^\star, \dots, x_p^\star)$ is given by (in matrix-vector form)
$\hat y_i \pm t_{\alpha/2, \ n - (p + 1)} \sqrt{\hat \sigma^2 x_i^\star (X^TX)^{-1}{x_i^\star}^T}$
In [[R]], use the `predict` function with an `interval` parameter to get confidence intervals for the fitted values.
```R
predict(lm_data, interval="confidence", level=0.95)
```
# interpretation of the confidence interval
Imagine we could fix the predictors in the training data and resample the response many times. If we fit new models at the same predictor values and recomputed the confidence intervals, $(1-\alpha)\%$ of the models would include the true value of the parameter. In practice, this is not possible as it would be next to impossible to find more statistical units that match the original data in all values but the response value.