Linear Models

Learning Objectives

• Understand how to find the best fitting line between the explanatory variable and the dependent variable.

• Know how to make confidence intervals for each regression coefficient and regression line.

• Know how to perform a hypothesis testing for each regression coefficient and regression line.

• Understand how to analyze real data using the simple and multiple regression models.

• Know how to interpret computer outcomes from R.

   

Introduction to Linear Models

Often, we are interested in understanding the relationship between one variable, called the dependent variable or response, on a number of other variables, called the independent variables or predictors. The dependent variable is typically denoted by $Y$, and the independent variables are denoted as $X_1,X_2,\ldots,X_p$.

Example: Suppose we are interested in predicting a final score of stat 415 class based on a number of predictors.

Response: 415 final score ($Y$)

Predictors: hours of study per week ($X_1$), number of previous math/stat classes ($X_2$), number of cups of coffee before the exam ($X_3$), final percentage of stat/math 414 ($X_4$), etc.

Goal: predict a 415 final score based on given predictor values. For example, what would be the expected 415 final score when

• hours of study per week = 3, number of previous math/stat classes = 5, number of cups of coffee before the exam =2, final percentage of stat/math 414 = 89%, or
• hours of study per week = 5, number of previous math/stat classes = 2, number of cups of coffee before the exam =5, final percentage of stat/math 414 = 72%?

Suppose

where $\epsilon$ is a random error ($E[\epsilon] = 0$). The goal is to learn $f$ from data. Once we know what $f$ is, we can predict the expected value of $Y$ given different predictor values $X_1,\dots,X_p$. In other words,

For example, the expected final score when hours of study per week = 3, number of previous math/stat classes = 5, number of cups of coffee before the exam =2, final percentage of stat/math 414 = 89% is $f(3,5,2,.89)$.

A linear model assumes that $f$ is linear. In other words,

Then the goal is to learn $\beta_0,\cdots,\beta_p$ from data. When $p=1$ (one explanatory model), a linear model is called as a simple linear regression model. Otherwise, the model is called multiple linear regression. In this class, we will mostly focus on simple linear regression.

The method of least squares

Suppose in the previous example only hours of study per week ($X$) matter in predicting the final score ($Y$). That is $Y= f(X) + \epsilon$. Furthermore, the data suggest that the relationship between $X$ and $Y$ is linear. We let $f$ be linear and assume $f(x) = \beta_0+\beta_1 x$, i.e.,

Final Score ($y_i$) = $\beta_0$ + $\beta_1$ Hours of study per week ($x_i$) + error, $i=1,\dots,n$.

Now the question is to find the best $b_0$ and $b_1$ based on the data $(x_1,y_1)$, $(x_2,y_2)$, ..., $(x_n,y_n)$ where $y_i$ is a measurement of $Y$ at $X=x_i$.

Remark: depending on the nature of the data, we may treat an independent variable $X$ either to be random or deterministic. For example, we may consider $(X_1, Y_1), . . . ,(X_n, Y_n)$ be $n$ independent and identically distributed (i.i.d.) random pair of random variables (a.k.a random design), or we may regard $X_1,\dots,X_n$ be as fixed numbers (a.k.a fixed design). In the fixed design case, $Y_1,\dots,Y_n$ are independent but not identically distributed, since each $Y_i$ has different mean value which depends on $X_i$.

In regression analysis, we are interested in estimating the expectation of $Y_i$ when $X_i = x$. In the case of fixed design, estimating the regression line is equivalent to the expectation of $Y_i$. On the other hand, in random design, the regression line is the conditional expectation of $Y_i$ given $X_i=x_i$.

In this class, we will always assume deterministic predictors (fixed design), and to emphasize this fact, we will use lower-case letters $x_1, . . . , x_n$ to denote predictors.

Suppose each $Y_i$ has the following representation:

where $\epsilon_i \sim N(0, \sigma^2)$ and each $x_i$ is a fixed number. Note this implies that $Y_i \sim N(\beta_0 + \beta_1x_i, \sigma^2)$. Given $(x_1,y_1),\dots,(x_n,y_n)$, how can we find $f(x) = \beta_0 + \beta_1 x$?

The basic idea is to choose a function of the predictors, $f$, so that the sum of squared distances from the data to the function values is minimized. Since $f$ is linear, we find $\beta_0,\beta_1$ such that $\sum_{i=1}^n (y_i - f(x_i))^2 = \sum_{i=1}^n (y_i - (\beta_0 +\beta_1x_i))^2$ is minimized.

Let $H(\beta_0,\beta_1) = \sum_{i=1}^n (y_i - (\beta_0 +\beta_1x_i))^2$. Then the minimizer $(\hat{\beta}_0,\hat{\beta}_1)$ of $H(\beta_0,\beta_1)$ are called the least square estimates of $β_0$ and $β_1$. The corresponding statistics (i.e., replace observed $y_i$ with random $Y_i$) are called least square estimators (LSE).

Theorem (LSE of $\beta_0$ and $\beta_1$ in simple linear regression) The least square estimators of $(\beta_0,\beta_1)$ are

Remark: It is important to note that both $\hat{\beta}_0$ and $\hat{\beta}_1$ are linear functions of the $Y_i$.

With some algebra, we can show

where $w_i = \frac{(x_i-\bar{x})}{\sum_{i=1}^{n}(x_i-\bar{x})^2}$, $i=1,\dots,n$.

Proof. Following the standard procedures for minimizing a function, we find the partial derivatives

$$\begin{align*} \frac{\partial S}{\partial \beta_0}&= -2\sum_{i=1}^{n}(y_i-\beta_0-\beta_1 x_i)\\ \frac{\partial S}{\partial \beta_1}&= -2\sum_{i=1}^{n}x_i(y_i-\beta_0-\beta_1 x_i). \end{align*}$$

Setting both to zero we get

$$\begin{align*} \sum_{i=1}^{n}y_i &= n \hat{\beta}_0 + \hat{\beta}_1\sum_{i=1}^{n}x_i\\ \sum_{i=1}^{n} x_i y_i &= \hat{\beta}_0 \sum_{i=1}^{n} x_i + \hat{\beta}_1 \sum_{i=1}^{n} x_i^2. \end{align*}$$

From these, we solve for ${\beta}_0$ and ${\beta}_1$, and treat them as estimates (i.e., putting on hats).

$$\begin{align*} \hat{\beta}_1 &= \frac{n\sum_{i=1}^{n}x_i y_i - \left(\sum_{i=1}^{n}x_i\right)\left(\sum_{i=1}^{n}y_i\right)}{n\sum_{i=1}^{n}x_i^2 - \left(\sum_{i=1}^{n}x_i\right)^2}\\ \hat{\beta}_0 &= \bar{y}-\hat{\beta}_1\bar{x}. \end{align*}$$

 

 

 

 

 

Remark: the same least squares method can be applied when we have more than one predictor variable. Namely, when $f(x_1,\dots,x_p) = \beta_0 + \beta x_1+\dots+\beta_px_p$, the least square estimates $(\hat{\beta}_0,\hat{\beta}_1,\dots,\hat{\beta}_p)$ of $(\beta_0,\beta_1,\dots,\beta_p)$ are

Theorem (MLE of $\sigma^2$, $\beta_0$ and $\beta_1$ in simple linear regression) The maximum likelihood estimates of $\beta_0$ and $\beta_1$ are the same as least squares estimates. The maximum likelihood estimator of $\sigma^2$ is given by

Theorem (Distribution of LSE of $\beta_0$ and $\beta_1$) Under the simple linear model (1) with Normally distributed errors $\epsilon_i \sim N(0,\sigma^2)$,

$\hat{\beta_1} \sim N(\beta_1, \frac{\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2})$ ,

$\hat{\beta_0}\sim N(\beta_0, \sigma^2(\frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i-\bar{x})^2}))$

Regarding optimality of Least squared estimators

• Least squared estimators $\hat{\beta_0}$ and $\hat{\beta}_1$ are unbiased for $\beta_0$ and $\beta_1$. In fact, under the linear model (1) and normally distributed $\epsilon_i \sim N(0,\sigma^2)$, $\hat{\beta_0}$ and $\hat{\beta}_1$ are minimum variance unbiased estimators (MVUEs).
• Note that without the normality assumption of error terms, the least squared estimators $\hat{\beta_0}$ and $\hat{\beta}_1$ are still unbiased. In what sense the least squared estimators $\hat{\beta_0}$ and $\hat{\beta_1}$ are optimal beyond being unbiased? If the distribution of an error term is unknown but has the constant variance (i.e., $\rm{Var}(\epsilon_i) = \sigma^2$ for all $i=1,\dots,n$), the least squared estimators $\hat{\beta_0}$ and $\hat{\beta}_1$ can be proved to be the Best Linear Unbiased Estimators (BLUEs) for $\beta_0$ and $\beta_1$ (Gauss-Markov Theorem). That is, the least square estimators are the estimators with the smallest variance among all unbiased estimators that are linear functions of the $Y_i$’s.

Confidence intervals and hypothesis testing related to regression models

Recall the distribution of LSE $\hat{\beta}_0$ and $\hat{\beta}_1$ under the linear model $Y_i = \beta_0 + \beta_1 x_i + \epsilon_i$ with $\epsilon_i \sim N(0,\sigma^2)$:

$\hat{\beta_1} \sim N(\beta_1, \frac{\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2})$ ,

$\hat{\beta_0}\sim N(\beta_0, \sigma^2(\frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i-\bar{x})^2})).$

Thus, the standardized regression coefficients are $N(0, 1)$ random variables, i.e.,

Unfortunately, $\sigma^2$ is unknown. Therefore, we replace $\sigma^2$ with a sample variance estimator $\widehat{\sigma^2}$ defined as

We note that the distribution of the residual sum of squares divided by $\sigma^2$ is $\chi^2(n-2)$, i.e.,

In particular, $S^2$ is an unbiased estimator of $\sigma^2$.

Moreover, we can show that $\bar{Y}, \hat{\beta}_1$, and $S^2$ are independent.

Replacing $\sigma^2$ with $S^2$ in (2), we have,

since,

These results allow us to construct confidence intervals for $β_0$ and $β_1$, and to carry out hypothesis tests using t-tests.

Confidence intervals for the regression parameters

Using

as pivotal quantities,

• a (two-sided) $1-\alpha$ confidence interval for $\beta_0$ is

  $$\hat{\beta}_0\pm t_{\alpha/2}(n-2)\sqrt{s^2(\frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i-\bar{x})^2})}$$

• a (two-sided) $1-\alpha$ confidence interval for $\beta_1$ is

  $$\hat{\beta}_1\pm t_{\alpha/2}(n-2)\sqrt{\frac{s^2}{\sum_{i=1}^n(x_i-\bar{x})^2}}$$

A $1-\alpha$ confidence interval for the average response value at $x^\star$

Often we are interested in estimating the average response value at a particular $X=x^\star$, i.e., estimate $E[Y] = β_0 + β_1x^\star$. Let $w = E[Y]$. We want to construct a $1-\alpha$ confidence interval for $w$.

The estimator of the expected value of $Y$ at $X = x^\star$ is $\hat{w} = \hat{\beta}_0 + \hat{\beta}_1x^\star$.

The distribution of the estimator $\hat{w}$ of $w$ is

Standardization gives

and by replacing unknown $\sigma^2$ with $S^2$, we obtain the following pivotal quantity:

Therefore, we can construct a (two-sided) $1-\alpha$ confidence interval for $w = \beta_0 + \beta_1 x^\star$ as

A $1-\alpha$ prediction interval for a new response value at $x^\star$

We are sometimes interested in predicting new value of $Y$ given $X=x^\star$. Where can we expect to see the next data point sampled? Note that even though we know the true regression line, we cannot perfectly predict the new value of $Y$ because of the existence of a random error. In particular, when $X = x^\star$,

where $\epsilon_{n+1} \sim N(0,\sigma^2)$.

Suppose we know the true regression line $\beta_0 + \beta_1 x^\star$ and also $\sigma^2$. Since $Y_{n+1} \sim N(\beta_0 + \beta_1 x^\star, \sigma^2)$, we can predict that $100(1-\alpha)$% of the times, the next value will be within $\beta_0 + \beta_1 x^\star \pm z_{\alpha/2} \sigma$.

Of course $(\beta_0,\beta_1, \sigma^2)$ are unknown, and we have to estimate these unknown quantities. The variation in the prediction of a new response depends on two components:

1. the variation due to estimating the mean $\beta_0 + \beta_1 x^\star$ with $\hat{\beta_0} + \hat{\beta}_1 x^\star$.
2. the variation due to new random error $\epsilon_{n+1}$.

Adding the two variance components, we obtain

A $1-\alpha$ prediction interval for $Y_{n+1}$ is

If you sample responses at $X=x^\star$ many times, we would expect that next value to lie within this prediction interval in 95% of the samples.

Remark: comparing confidence intervals and prediction intervals

• $(1-\alpha)$ prediction intervals are always wider than $(1-\alpha)$ confidence intervals

• As we sample an increasing number of points $(x_i,Y_i)$, the width of a $(1-\alpha)$ confidence interval for $\beta_0 + \beta_1 x^\star$ will go to $0$

• As we sample an increasing number of points $(x_i,Y_i)$, the width of a $(1-\alpha)$ prediction interval for future draws $Y_{n+1}$ will not go to $0$

   

Hypothesis tests for the regression parameters

Suppose we are interested in testing whether the slope $\beta_1$ is the same as the hypothesized value $\beta_{10}$ or not (usually we are interested in testing whether the slope $\beta_1$ is zero or not, i.e., $\beta_{10}=0$)

1. $H_0: \beta_1 = \beta_{10}$, vs. $H_1: \beta_1 \ne \beta_{10}$.

2. Test statistic

   $$T = \frac{\hat{\beta_1}-\beta_{10}}{\sqrt{\frac{S^2}{\sum_{i=1}^n(x_i-\bar{x})^2}}}$$

3. Under $H_0$, $T \sim t(n-2)$.

4. Reject $H_0$ if the observed test statistic $t_{obs}$ is in the rejection region. For example, when $\alpha=.05$, reject $H_0$ if $t_{obs}\le -t_{.025}(n-2)$ or $t_{obs}\ge t_{.025}(n-2)$.

The same steps are applied for testing about $\beta_0$ ($H_0: \beta_0 = \beta_{00}$), with the test statistic

Simulated example

​x
# true parameters
beta0 <- 1
beta1 <- 10
sigma <- 40
​
# generate a random sample of size n
set.seed(1234)
n <- 100
x <- rnorm(n, 50, 5)
y <- rnorm(n, beta0 + beta1 * x, sigma)
plot(x, y, pch = 16, col = 4)

By hand

xxxxxxxxxx
> cor(x, y)
[1] 0.7659389
> mean(x)
[1] 49.21619
> mean(y)
[1] 494.8116
> var(x)
[1] 25.22075
> var(y)
[1] 4121.471

R output

xxxxxxxxxx
# Fitting a linear model
model <- lm(y ~ x)
​
# LM output
> summary(model)
​
Call:
lm(formula = y ~ x)
​
Residuals:
     Min       1Q   Median       3Q      Max 
-115.450  -24.560    0.094   23.458  119.509 
​
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.9201    41.0697   0.315    0.754    
x             9.7913     0.8302  11.794   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
​
Residual standard error: 41.48 on 98 degrees of freedom
Multiple R-squared:  0.5867,    Adjusted R-squared:  0.5824 
F-statistic: 139.1 on 1 and 98 DF,  p-value: < 2.2e-16
​
​
# Confidence intervals
> confint(model,level = .95)
                 2.5 %   97.5 %
(Intercept) -68.581325 94.42161
x             8.143801 11.43884
​
# 95% confidence interval when x = 50
> predict(model,newdata = data.frame(x=50),interval = "confidence",level = .95)
       fit      lwr      upr
1 502.4862 494.1531 510.8192
​
# 95% prediction interval when x = 50
> predict(model,newdata = data.frame(x=50),interval = "prediction",level = .95)
       fit      lwr     upr
1 502.4862 419.7414 585.231