3
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ The parameter estimates are correlated. vcov(my_fitted_data) shows you the variance-covariance matrix of the parameter estimates, cov2cor(vcov(my_fitted_data)) the correlation matrix. Your approach ignores this correlation. $\endgroup$ – Wolfgang Jan 29 '17 at 21:32
  • $\begingroup$ @Wolfgang I don't really understand what that means (how to interpret the output of these commands), unfortunately....Is the take-home message: Use the R function for the correct confidence interval? $\endgroup$ – l7ll7 Jan 29 '17 at 21:40
  • 1
    $\begingroup$ @user10324 I asked a similar question a while ago. The answer there should give you what you need. stats.stackexchange.com/questions/10439/… $\endgroup$ – mark999 Jan 29 '17 at 21:42
  • 1
    $\begingroup$ Adding to other comments, you are misunderstanding confidence intervals, moreover, you seem to be asking about prediction intervals. $\endgroup$ – Tim Jan 29 '17 at 21:44
  • 1
    $\begingroup$ @user10324 still CI's do not give you "95% probability that x has to lie in the interval". $\endgroup$ – Tim Jan 29 '17 at 21:47
8
$\begingroup$

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
$\endgroup$
  • $\begingroup$ +1 especially the R example clarified it for me. thanks a lot!! $\endgroup$ – l7ll7 Jan 29 '17 at 23:12
  • $\begingroup$ @user10324 You're welcome. $\endgroup$ – mark999 Jan 29 '17 at 23:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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