2
$\begingroup$

EDIT: It appears that for this particular problem I made a simple arithmetic error. I am now getting the same results. What prompted me to write this question was that this is the third time that I've had this problem. However, in the past, I was mistakenly using interval="predict" when I should have been using interval="confidence".

I am trying to compute a 95% confidence interval for a mean response on a small dataset, yet when I calculate this manually I get a very different interval. How is R calculating the interval when using predict.lm? Am I using the wrong function call? I verified that my manual calculation is correct with the solutions manual.

Computing a 95% confidence interval for a mean response at $x_0 = 170$ for the following dataset:

      x        y
1   170   162.50
2   140   144.00
3   180   147.50
4   160   163.50
5   170   192.00
6   150   171.75
7   170   162.00
8   110   104.83
9   120   105.67
10  130   117.58
11  120   140.25
12  140   150.17
13  160   165.17

I calculated the confidence interval using the following equation from my book: $$ \hat\mu_{y|x_0}\pm t_{\alpha/2, n-2}\sqrt{MS_{\rm res}\left(\frac 1 n + \frac{(x_0-\bar x)^2}{S_{xx}} \right)} $$ And I got an interval of $168.36\pm23.465$ (or alternatively, $(144.895, 191.825)$). When I did this manual calculation, the MSRes I calculated is the same as the output from calling anova(model). Additionally, I calculated Sxx in R as follows:

# Results are equivalent:
sum((data$x - mean(data$x))^2) # = 6230
(length(data$x)-1)*var(data$x) # = 6230

However, when I calculate this in R, I get the same fitted value but a different interval:

> predict.lm(data.lm, newdata=data.frame(x=170), interval="confidence", level=0.95)
       fit      lwr      upr
1 168.3611 153.9071 182.8152

Why is R giving me an incorrect interval? Or is my textbook simply wrong? This also applies to a 95% prediction interval on the same $x_0$, using the following equation: $$ \hat y_0 \pm t_{\alpha/2, n-2}\sqrt{MS_{\rm res}\left(1 + \frac 1 n + \frac{(x_0-\bar x)^2}{S_{xx}} \right)} $$

$\endgroup$
3
  • $\begingroup$ (For whomever: this seems like as much a statistical question as an R code question; I think it should be on topic here.) $\endgroup$ Nov 25, 2015 at 2:17
  • $\begingroup$ It seems that I made a silly error. Please see my edit at the top of my question. Should this question be closed at this point? $\endgroup$
    – ssahli
    Nov 25, 2015 at 3:47
  • $\begingroup$ I see no need for this question to be closed. IMO, clarifying these calculations has intrinsic value irrespective of the original motivation. $\endgroup$ Nov 25, 2015 at 3:49

1 Answer 1

4
$\begingroup$

It's hard to guess what you might have done wrong since Sxx and MSres are correct. Perhaps the t-value? Here is the calculation done by hand in R to confirm that the R answer agrees with the formula:

data <- structure(list(x = c(170, 140, 180, 160, 170, 150, 170, 110, 
120, 130, 120, 140, 160), y = c(162.5, 144, 147.5, 163.5, 192, 
171.75, 162, 104.83, 105.67, 117.58, 140.25, 150.17, 165.17)), .Names = c("x", 
"y"), row.names = c(NA, -13L), class = "data.frame")

data.lm <- lm(y~x, data=data)

# get prediction
pred <- predict.lm(data.lm, newdata=data.frame(x=170), interval="confidence", level=0.95)

# correct t-value
tval <- qt((1-0.95)/2, df=13-2)

# correct Sxx
Sxx <- sum((data$x - mean(data$x))^2)

# correct MSres - note division is by the number of df (but you have this right anyway)
MSres <- sum(data.lm$residuals^2)/11

# confidence interval calculated by hand
sqrt(MSres * (1/13 + (170 - mean(x))^2/Sxx)) * tval * c(1, -1) + pred[1]

# compare with R interval calculation
pred
$\endgroup$
1
  • $\begingroup$ Thanks for the input. I made a silly mistake in my calculation, and I ended with the same result as you. Also, it appears that my misconception is with what interval type to use for a given problem--the mistakes I have made in the past is because I was using interval="predict" when in fact I should have been using interval="confidence". I would make the mistake of assuming predict whenever the problem gave me a new value for x. Thank you for helping me see this! Our instructor couldn't help me since he is not familiar with R. $\endgroup$
    – ssahli
    Nov 25, 2015 at 3:39

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.