Lesson 3.3 - Inference in SLR

Gauss-Markov Assumptions

In the last lesson, we introduced a set of assumptions that we will commonly make regarding the error term in our hypothetical model \(Y = \beta_0 + \beta_1 X + \varepsilon\). These assumptions are as follows:

1. The relationship between \(X\) and \(Y\) is determined by a linear model of the form \(Y = \beta_0 + \beta_1 X + \varepsilon\).

2. Seperate observations of \(\varepsilon\) are independent from one another.

3. The error term \(\varepsilon\) is normally distributed with a mean of 0 and standard deviation \(\sigma\). That is, \(\varepsilon \sim N(0,\sigma^2)\).

4. The error term \(\varepsilon\) is independent from \(X\). In particular, the variance \(\sigma^2\) does not depend on \(X\).

Note that these assumptions are not required in order to find or even apply the OLS regression line. They are, however, necessary if we wish to use the model for inferential purposes.

Residual Standard Error

If we are only interested in using a simple linear regression model for the purposes of making point predictions, then we need only to find estimates for the parameters \(\beta_0\) and \(\beta_1\). However, if we wish to quantify the uncertainty in these predictions, or if we would like to conduct any sort of statistical inference using our model, then we also need an estimate for \(\mathrm{Var}\left[\varepsilon \right] = \sigma^2\). Such an estimate is provided by the quantity \(s^2\) defined below:

\[s^2 = \frac{SSE}{n-2}\] If the Gauss-Markov assumptions hold, then \(s^2\) is an unbiased estimator of \(\sigma^2\).

The square root of \(s^2\) is referred to as the residual standard error. It is typically denoted by either \(s\) or \(RSE\).

Distribution of Parameter Estimates

Assume that we fix several \(X\) values, \(x_1, ..., x_n\) and that we use these \(X\) values to generate \(Y\) values \(y_1, ..., y_n\) according to a hypothetical model \(Y = \beta_0 + \beta_1 X + \varepsilon\). Suppose we then calculate parameter estimates \(\hat\beta_0\) and \(\hat\beta_1\).

If we repeat this process, we will get different values of \(y_1, ..., y_n\) each time, since the error terms \(\varepsilon_i\) are assumed to be random. Since our \(Y\) values will vary each time, so too will our parameter estimates \(\hat\beta_0\) and \(\hat\beta_1\). In this way, our parameter estimates are random variables whose distribution depends on that of \(\varepsilon\).

It can be shown that if the Gauss-Markov assumptions hold, then the distibutions of \(\hat\beta_0\) and \(\hat\beta_1\) are as follows:

• \(\hat\beta_0\) follows a normal distribution with:

  • \(\mu_{\hat\beta_0} = \mathrm{E}[\hat\beta_0] = \beta_0\)

  • \(\sigma_{\hat\beta_0}^2 = \mathrm{Var}[\hat\beta_0] = \sigma^2 \left[\frac{1}{n} + \frac{\bar x^2}{S_{XX}} \right]\)

• \(\hat\beta_1\) follows a normal distribution with:

  • \(\mu_{\hat\beta_1} = \mathrm{E}[\hat\beta_1] = \beta_1\)

  • \(\sigma_{\hat\beta_1}^2 = \mathrm{Var}[\hat\beta_1] = \frac{\sigma^2}{S_{XX}}\)

Unfortunately, we don’t typically know \(\mathrm{Var}[\varepsilon] = \sigma^2\). If we do not know \(\sigma^2\), then we will need to approximate it using \(s^2 = \frac{SSE}{n-2}\). Replacing \(\sigma^2\) with \(s^2\) in the variance formulas above, and then taking square roots results in the appoximations for \(\mathrm{SD}[\hat\beta_0]\) and \(\mathrm{SD}[\hat\beta_1]\) shown below. These approximations are called the standard errors of the parameter estimates.

• \(\mathrm{SD}[\hat\beta_0] \approx \mathrm{SE}[\hat\beta_0] = s \sqrt{\frac{1}{n} + \frac{\bar x^2}{S_{XX}} }\)

• \(\mathrm{SD}[\hat\beta_1] \approx \mathrm{SE}[\hat\beta_1] = \frac{s}{\sqrt{S_{XX}}}\)

Hypothesis Tests for Parameters

It is often useful to be able to test the hypothesis that a particular parameter is equal to a specific number. In other words, we want to be able to conduct hypothesis tests of the following form:

\(\hspace{30pt} H_0: \beta_i = k\)

\(\hspace{30pt} H_A: \beta_i \neq k\)

We can perform such a test using the following test statistic:

\(\large\hspace{30pt} t = \frac{\hat\beta_{i} - k}{\mathrm{SE}[\hat\beta_i]}\)

This test statistic will follow a t distribution with \(n-2\) degrees of freedom.

Although we can test either regression parameter for any value of \(k\), the test that we will be most interested in is the following:

\(\hspace{30pt} H_0: \beta_1 = 0\)

\(\hspace{30pt} H_A: \beta_1 \neq 0\)

If we are unable to reject the null hypothesis in this test, then we can not be confident that the parameter \(\beta_1\) is non-zero. If \(\beta_1 = 0\), then our hypothetical model reduces to \(Y = \beta_0 + e\), and there is no relationship between \(X\) and \(Y\).

Hypothesis Tests in R

Fortunately, R will test for the regression parameters being zero when you call the lm() function. The results of these tests are displayed in the summary output for the model. To demonstrate this, we will consider the following simulated data.

set.seed(45)

n = 100
x <- runif(n, 0, 10)
y <- 8 + 1.3 * x + rnorm(n, 0, 15)

m <- lm(y ~ x)

summary(m)

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-31.210 -11.303   2.015  10.053  54.098 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   9.4861     3.2981   2.876  0.00494 **
x             0.9987     0.5707   1.750  0.08326 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.61 on 98 degrees of freedom
Multiple R-squared:  0.0303,    Adjusted R-squared:  0.02041 
F-statistic: 3.062 on 1 and 98 DF,  p-value: 0.08326

The summary above tells use the the p-value for the following test is equal to \(p = 0.00494\).

\(\hspace{30pt} H_0: \beta_0 = 0\)

\(\hspace{30pt} H_A: \beta_0 \neq 0\)

The summary also tells use the the p-value for the test below is equal to \(p = 0.08326\).

\(\hspace{30pt} H_0: \beta_1 = 0\)

\(\hspace{30pt} H_A: \beta_1 \neq 0\)

If we had decided in advance to use a significance level of \(\alpha = 0.05\), then we would be able to reject the hypothesis that the intercept was non-zero, but we would not have enough evidence to reject the hypothesis that the slope was non-zero.

The figure below displays a plot of the data set, along with the population and fitted models.

plot(y ~ x, pch=21, bg="cyan", col="black", cex=1)
abline(m$coefficients, col="darkred", lwd=2)

abline(7, 1.3, col="darkblue", lwd=2, lty=2)

Confidence and Prediction Intervals

Assume that we have a fitted simple linear regression model \(\hat Y = \beta_0 + \beta_1 X\). Given a predictor value \(X = x_0\), we can use the the model to generate a point prediction \(y_0\) for the response variable \(Y\). It is often useful to be able to quantify the amount of uncertainty there is a prediction of this sort. We will provide two ways of describing this uncertainty: confidence intervals and prediction intervals.

Confidence Intervals

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

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

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 \(\varepsilon\).

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

\[\hat y_0 \pm t_{(n-2, \alpha / 2)} \cdot s \sqrt{\frac{1}{n} + \frac{(x_0 - \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_0\) is a range of \(Y\) values that we are \(100\beta\%\) confident will contain \(y_0 = \beta_0 + \beta_1 x_0 + \varepsilon_0\).

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 \(\varepsilon\).

Prediction intervals are ALWAYS wider than confidence intervals.

The \(100\beta\%\) confidence interval for \(y_0 = \beta_0 + \beta_1 x_0 + \varepsilon_0\) is given by:

\[\hat y_0 \pm t_{(n-2, \alpha / 2)} \cdot s \sqrt{1 + \frac{1}{n} + \frac{(x_0 - \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_0] = \beta_0 + \beta_1 x_0\) is given by:

\[\hat y_0 \pm t_{(n-2, \alpha / 2)} \cdot s \sqrt{\frac{1}{n} + \frac{(x_0 - \bar x)^2}{(n-1)s_X^2}} \] The \(100\beta\%\) confidence interval for \(y_0 = \beta_0 + \beta_1 x_0 + \varepsilon_0\) 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