How to compute a prediction interval from ordinary least squares regression output alone?

Question

This question has been haunting me for a long time, When I'm given an R output of linear regression, and asked to calculate 95% prediction interval, I feel there's something missing.

In this output, how am I supposed to calculate 95% prediction interval for X_b value of 10?

I can evaluate y_hat and MSE from the output but I fail to understand how to get the mean of X_b and Sxx from here.

Answer 1

$S_{xx},$ the sum of squares of the explanatory variable, is easy to obtain from the formula

$$\operatorname{se}(\hat\beta_1) = \sqrt{\frac{MS_{Res}}{S_{xx}}}$$

where the left hand side is the standard error of the slope, given as $1373$ in the question, and $MS_{Res}$ is the mean squared residual, whose square root (the "residual standard error" is given as $36600$ in the question.

The mean of the explanatory variable can almost be recovered from the formula for the estimated sampling variance of the intercept,

$$\widehat{\operatorname{Var}}(\hat\beta_0) = MS_{Res}\left(\frac{1}{n} + \frac{\bar x^2}{S_{xx}}\right).$$

In the question, the left hand side is the square of the standard error, $\widehat{\operatorname{Var}}(\hat\beta_0) = 8004^2$ and $n = 98 + 2$ is found by adding the number of estimated coefficients to the "degrees of freedom" reported for the $F$ ratio statistic. Solving this for $\bar x$ usually gives two possible values. Unless you have some sense of what the value should be (one solution is positive and the other is negative), you're stuck (because, as you clearly are aware, the prediction interval at any value $x_0$ depends on its distance from $\bar x$ and the only value where that distance does not depend on the solution is $x_0=0$).

---

As an example of the problem, here is R code to manufacture two different datasets with differing values of $\bar x$ and identical ordinary least squares output.

x <- seq(2, 10, length.out = 30)
y <- x + rnorm(length(x))

fit <- lm(y ~ x)

b <- coefficients(fit)
x.bar <- mean(x)
x <- x - 2 * x.bar
y <- y - 2 * b[2] * x.bar

all.equal(summary(fit), summary(lm(y ~ x)))

It alters the initial data by subtracting $2\bar x$ from all $x$ values, subtracting $2\hat\beta_1 \bar x$ from all $y$ values to keep the coefficient estimates the same, and comparing their summaries. Its output is

Component “cov.unscaled”: Mean relative difference: 0.4387476 

That is, the only difference between the two datasets lies in the estimated covariance between $\hat\beta_0$ and $\hat\beta_1$ -- but that is not part of your regression output. (The sign of this covariance will differ in the two datasets.) If it could be recovered from the output then some other numbers in the output would have to differ, too, but that's not the case.

Here is a plot of the original data (blue; $\bar x = 6$) and their transformed version (red; $\bar x = -6$). The line is the common least squares fit.