Lesson 11 - Prediction and Confidence Intervals

Confidence Intervals

A confidence interval can be thought of as a vertical interval that we are reasonably confident the line representing our theoretical model passes through.

More precisely, a \(100\beta\%\) confidence interval given \(X = x\) is a range of \(Y\) values that we are \(100\beta\%\) confident will contain \(\mathrm{E}[Y | X=x] = \beta_0 + \beta_1 x\).

Such an interval accounts for the possible errors in the approximations of \(\hat\beta_0\) and \(\hat\beta_1\), but does not account for the uncertainty resulting from the noise term \(e\).

The \(100\beta\%\) confidence interval for \(\mathrm{E}[Y | X=x] = \beta_0 + \beta_1 x\) is given by:

\[\hat y \pm t_{(n-2, \alpha / 2)} \cdot s \sqrt{\frac{1}{n} + \frac{(x - \bar x)^2}{(n-1)s_X^2}} \] The \(\alpha\) in the formula above is given by \(\alpha = 1 - \beta\).

Prediction Intervals

A prediction interval can be thought of as a vertical interval that we are reasonably confident will contain the \(y\) coordinate of an observation with a given \(x\) coordinate.

More precisely, a \(100\beta\%\) prediction interval given \(X = x\) is a range of \(Y\) values that we are \(100\beta\%\) confident will contain \(y = \beta_0 + \beta_1 x + e\).

Such an interval accounts for the possible errors in the approximations of \(\hat\beta_0\) and \(\hat\beta_1\), as well as the uncertainty resulting from the noise term \(e\).

Prediction intervals are ALWAYS wider than confidence intervals.

The \(100\beta\%\) confidence interval for \(y = \beta_0 + \beta_1 x + e\) is given by:

\[\hat y \pm t_{(n-2, \alpha / 2)} \cdot s \sqrt{1 + \frac{1}{n} + \frac{(x - \bar x)^2}{(n-1)s_X^2}} \] The \(\alpha\) in the formula above is given by \(\alpha = 1 - \beta\).

Summary of Formulas

The \(100\beta\%\) confidence interval for \(\mathrm{E}[Y | X=x] = \beta_0 + \beta_1 x\) is given by:

\[\hat y \pm t_{(n-2, \alpha / 2)} \cdot s \sqrt{\frac{1}{n} + \frac{(x - \bar x)^2}{(n-1)s_X^2}} \] The \(100\beta\%\) confidence interval for \(y = \beta_0 + \beta_1 x + e\) is given by:

\[\hat y \pm t_{(n-2, \alpha / 2)} \cdot s \sqrt{1 + \frac{1}{n} + \frac{(x - \bar x)^2}{(n-1)s_X^2}} \] In each formula, we have that \(\alpha = 1 - \beta\).

Simulated Example

set.seed(1)
n = 10
x <- 8 + 1:n / 5
y <- 2*x + 3 + rnorm(n,0,0.8)
model <- lm(y~x)
          
plot(y ~ x, xlim=c(8,10.1), ylim=c(17,26), col="black", bg="red", pch=21, cex=1)
abline(model$coefficients, col="blue", lwd="2")

Let find the 95% confidence and prediction intervals for \(X=9\).

n = 10
x_ = 9
yhat = model$coefficients[1] + model$coefficients[2] * x_
s2 = sum((model$residuals)^2) / (n - 2)
s = sqrt(s2)
ci_sqrt = sqrt(1/n + (x_ - mean(x))^2 / ((n-1) * var(x)) )
pi_sqrt = sqrt(1 + 1/n + (x_ - mean(x))^2 / ((n-1) * var(x)) )
t = qt(0.975, n-2)
ci = c(yhat - t * s * ci_sqrt, yhat + t * s * ci_sqrt)
pi = c(yhat - t * s * pi_sqrt, yhat + t * s * pi_sqrt)
names(ci) = c('lower', 'upper')
names(pi) = c('lower', 'upper')
rbind(ci, pi)

      lower    upper
ci 20.60478 21.56296
pi 19.51630 22.65144

Let’s check our answer using the predict() function.

nd <- data.frame(x=9)
conf <- predict(model, interval="confidence", level=0.95, newdata=nd)
pred <- predict(model, interval="prediction", level=0.95, newdata=nd)
rbind(conf, pred)

       fit      lwr      upr
1 21.08387 20.60478 21.56296
1 21.08387 19.51630 22.65144

Example Using Pearson Data

We will now walk through an another example of using R to create prediction and confidence intervals. In this example, we will once again use the Pearson data set.

We start by loading the data and creating our model.

myData <- read.table("father_son.txt", sep="\t", header=TRUE)
mod <- lm(sheight ~ fheight, myData)
summary(mod)

Call:
lm(formula = sheight ~ fheight, data = myData)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8772 -1.5144 -0.0079  1.6285  8.9685 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 33.88660    1.83235   18.49   <2e-16 ***
fheight      0.51409    0.02705   19.01   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.437 on 1076 degrees of freedom
Multiple R-squared:  0.2513,    Adjusted R-squared:  0.2506 
F-statistic: 361.2 on 1 and 1076 DF,  p-value: < 2.2e-16

We generate a scatter plot of the data, with the regression line included.

plot(myData$sheight ~ myData$fheight, pch=21, bg="cyan", col="black", cex=1,
     main="Scatter Plot", xlab="Father Height", ylab="Son Height")
abline(mod$coefficients, col="darkred", lwd=2)

The following code uses our model to predict the height of an adult male whose father is 5 feet tall.

nd1 = data.frame(fheight=60)
predict(mod, newdata = nd)

Error in eval(expr, envir, enclos) : object 'fheight' not found

Let’s create a 95% confidence interval for the height of an adult male whose father is 5 feet tall. Recall that this interval has a 95% chance of including the true population mean height for such individuals. It does NOT claim to contain the heights of 95% of these individuals.

predict(mod, newdata=data.frame(fheight=60), interval="confidence", level=.95)

       fit      lwr      upr
1 64.73219 64.29899 65.16538

If we want an interval that we believe will contain the heights of 95% of all adult males whose fathers are 5 feet tall, we need to consider a prediction interval instead.

predict(mod, data.frame(fheight=60), interval="predict", level=.95)

       fit      lwr      upr
1 64.73219 59.93166 69.53271

Let’s generate a prediction interval for the height of an adult male whose father is 6 feet tall.

predict(mod, data.frame(fheight=72), interval="predict", level=.95)

      fit      lwr      upr
1 70.9013 66.11267 75.68993

We can also use predict() to create predictions for several values of the predictor at once.

predict(mod, newdata=data.frame(fheight=c(60, 62, 64, 66, 68, 70, 72)), interval="predict", level=.95)

       fit      lwr      upr
1 64.73219 59.93166 69.53271
2 65.76037 60.96770 70.55304
3 66.78856 62.00140 71.57572
4 67.81674 63.03275 72.60074
5 68.84493 64.06175 73.62812
6 69.87312 65.08839 74.65785
7 70.90130 66.11267 75.68993