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. enter image description here


1 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.

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.