The variance in \(Y\) is explained by \(X\)!
This intuition is going to be very important: there’s variance in \(y\), but incorporating \(x\) explains some of that variability. Instead of variance around the mean of \(y\), we’re looking at the variance around the line.
Why is \(\sigma^2 < \sigma^2_y\)? Imagine a horizontal line at \(\bar y\) (also show that \(\hat\beta_0 = \bar y\) is the least squares estimator of the model \(y_i = \beta_0 + \epsilon\)). The variance above and below the line is simply the variance of \(y\), i.e. \(\sigma^2_y\). This is the absolute worst case scenario. If that line is moved to fit the data as well as possible, then the variance above and below the line will be smaller!
\[\begin{align*} \hat\epsilon_i &= y_i - \hat y_i \\&= y_i - \bar y - (\hat y_i - \bar y)\\ \implies ... \implies \sum_{i=1}^n\hat\epsilon_i^2 &= \sum_{i=1}^n(y_i-\bar y)^2 - \sum_{i=1}^n(\hat y_i - \bar y)^2 \end{align*}\] where we’ve simply added and subtracted \(\bar y\). The final line skips a few steps (try them yourself, or see textbook page 28)!
Note that the “…” is actually fairly difficult to prove. You will need to plug in the least squares estimates \(\hat\beta_0\) and \(\hat\beta_1\) in order to get to the final answer.
The last line is often written as: \(SS_E = SS_T - SS_{Reg}\).
\[ SS_E = SS_T - SS_{Reg} \]
Def: The number of “pieces of information” from \(y_1, y_2, ..., y_n\) to construct a new number.
Estimating one parameter takes away one degree of freedom for the rest!
| Source of Variation | Degress of Freedom \(df\) | Sum of Squares (\(SS\)) | Mean Square (\(MS\)) |
|---|---|---|---|
| Regression | 1 | \(\sum_{i=1}^n(\hat y_i - \bar y)^2\) | \(MS_{Reg}\) |
| Error | \(n-2\) | \(\sum_{i=1}^n(y_i - \hat y_i)^2\) | \(MS_E=s^2\) |
| Total (corrected) | \(n-1\) | \(\sum_{i=1}^n(y_i - \bar y)^2\) | \(s_y^2\) |
Okay, so, just test for \(SS_{Reg} = 0\)?
But how??? We need some measure of how far from 0 is statistically significant!!!
Recall that \(SS_E = \sum_{i=1}^n(y_i - \hat y_i)^2\).
Statistics is the process of adding context to numbers. We expect a model in which there’s no relationship to have an \(SS_{Reg}\) of 0. Due to random sampling, we’re never actually going to get this.
x <- runif(50, -5, 5)
ssreg_vals <- double(1000)
for (i in 1:1000) {
# There is no relationship here!
y <- 4 + 0*x + rnorm(50, 0, 3)
mylm <- lm(y ~ x)
yhat <- predict(mylm)
ybar <- mean(y)
ssreg_vals[i] <- sum((yhat - ybar)^2)
}
hist(ssreg_vals)
sum(ssreg_vals == 0)[1] 0There’s no relationship, but none of the estimated values of \(SS_{Reg}\) are 0! Some of them were as high as 80, even though there’s truly no relationship!
So, what value of \(SS_{Reg}\) counts as “too large”? We actually don’t test this directly. Because of the way we chose to do statistics, it’s actually quite difficult to create a test at the boundary of possible values (in this case, testing if \(SS_{Reg} = 0\) when \(SS_{Reg}\) can only be greater than or equal to 0).
The whole point of an ANOVA table is to get down to:
\[ F_{df_{Reg}, df_E} = \dfrac{MS_{Reg}}{MS_E} = \dfrac{MS_{Reg}}{s^2} \]
A study was made on the effect of temperature on the yield of a chemical process. The data are shown to the right.
library(aprean3) # Data from textbook
head(dse03a) # Data Set Exercise Ch03A x y
1 -5 1
2 -4 5
3 -3 4
4 -2 7
5 -1 10
6 0 8nrow(dse03a)[1] 11Answers:
# 1.
coef(lm(y ~ x, data = dse03a))(Intercept) x
9.272727 1.436364 # 2.
anova(lm(y ~ x, data = dse03a))Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 226.945 226.94 96.18 4.207e-06 ***
Residuals 9 21.236 2.36
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 3.
confint(lm(y ~ x, data = dse03a), conf.level = 0.95) 2.5 % 97.5 %
(Intercept) 8.225007 10.320447
x 1.105045 1.767682As homework, calculate all of these using the formulas in R!
Consider the following four data sets
| data_set | xbar | sigma_x | ybar | sigma_y | corr |
|---|---|---|---|---|---|
| 1 | 9 | 3.316625 | 7.500909 | 2.031568 | 0.8164205 |
| 2 | 9 | 3.316625 | 7.500909 | 2.031657 | 0.8162365 |
| 3 | 9 | 3.316625 | 7.500000 | 2.030424 | 0.8162867 |
| 4 | 9 | 3.316625 | 7.500909 | 2.030578 | 0.8165214 |
Can you guess what the plots look like?
\[\begin{align*} R^2 &= \frac{SSReg}{SST}\\& = \frac{S_{XY}^2}{S_{XX}S_{YY}} \end{align*}\]
layout(mat = matrix(c(1,2,3), nrow = 1), widths = c(0.5,1,1))
set.seed(18)
x <- runif(25, 0, 10)
y <- rnorm(25, 2 + 5*x, 6)
plot(rep(1, 25), y, xlab = "y", ylab = "y has variance", xaxt = "n")
abline(h = mean(y))
axis(2, mean(y), bquote(bar(y)), las = 1)
plot(x, y, ylab = "There's variance around the line")
abline(lm(y~x))
abline(h = mean(y))
axis(2, mean(y), bquote(bar(y)), las = 1)
mids <- predict(lm(y~x))
for(i in seq_along(mids)){
lines(x = rep(x[i], 2), y = c(y[i], mids[i]), col = 1)
}
mids <- predict(lm(y~x))
plot(mids ~ x, type = "n", ylab = "The line varies around the mean of y!")
for(i in seq_along(mids)){
lines(x = rep(x[i], 2), y = c(mean(y), mids[i]))
}
axis(2, at = mean(y), labels = bquote(bar(y)), las = 1)
abline(h = mean(y))
We must check our assumptions!
The statistical tests try to give a p-value for the hypothesis that the residuals are normal, but looking at the residual plot will always be superior!
Why \(\hat{\underline{y}}\) instead of \(\underline y\)?
Why not \(\underline x\)?
library(ggplot2)
theme_set(theme_bw())
library(patchwork)
library(broom)
library(palmerpenguins)
penguins <- penguins[complete.cases(penguins),]
g1 <- ggplot(penguins) +
aes(x = flipper_length_mm, y = body_mass_g) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, formula = y~x) +
labs(x = "Flipper Length (mm)",
y = "Body Mass (g)",
title = "y versus x")
plm <- lm(body_mass_g ~ flipper_length_mm, data = penguins)
p2 <- augment(plm)
g2 <- ggplot(p2) +
aes(x = .fitted, y = .resid) +
geom_point() +
geom_hline(yintercept = 0, col = "grey") +
labs(x = "Fitted", y = "Residuals", title = "Residuals versus Fitted")
g1 / g2
I just want to add a little more context to the “unit change” idea. In a linear model, the estimated height of the line is \(\hat y_i = \hat\beta_0 + \hat\beta_1x_i\). To go from Celcius to Fahrenheit, we use the equation \(F_i = 32 + \frac{9}{5}C_i\), where \(C_i\) is the Celcius value that we have and \(F_i\) is the Fahrenheit value that we get. A line is literally a change of units!
Mathematics is the process of making assumptions and seeing if we can break them.
x <- runif(200, 0, 10)
y0 <- 2 - 3*x
y <- y0 + rnorm(length(x), 0, 2 * x)
g1 <- ggplot() +
aes(x = x, y = y) +
geom_point() +
geom_smooth(se = FALSE, method = "lm", formula = y~x) +
labs(x = NULL, y = NULL, title = "y versus x")
xydf <- augment(lm(y ~ x))
g2 <- ggplot(xydf) +
aes(x = .fitted, y = .resid) +
geom_point() +
geom_hline(yintercept = 0, col = "grey") +
labs(x = "Fitted", y = "Residuals", title = "Residuals versus Fitted")
g1 + g2
The idea of heteroscedastic data is that \(V(\epsilon_i) = \sigma^2_i\), which depends on \(i\). The most obvious case is shown above, where the variance increases along the values of \(x\). This pattern is obviously present in the plot of \(y\) versus \(x\), but this can be very hard to see in higher dimensions!
## fig-height: 4
## fig-width: 8
g1 <- ggplot(mtcars) +
aes(x = disp, y = mpg) +
geom_point() +
geom_smooth(se = FALSE, method = "lm", formula = y~x) +
labs(x = "Engine Displacement", y = "Miles per Gallon", title = "y versus x")
mtdf <- augment(lm(mpg ~ disp, data = mtcars))
g2 <- ggplot(mtdf) +
aes(x = .fitted, y = .resid) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, col = "grey") +
labs(x = "Fitted", y = "Residuals", title = "Residuals versus Fitted - non-linear trend?")
g1 + g2
The plots above show a (mild) violation of the assumption that \(y_i = \beta_0 + \beta_1x_i + \epsilon_i\), i.e., that the true trend is linear. This can be seen in the plot of \(y\) against \(x\), but it’s much more obvious in the residual plot!
Note that, in this course, the phrase “The Residual Plot” always refers to a plot of residuals against fitted. There are other residual plots, but the residual plot is resids versus fitted.
Consider the data (2,3,3,4,5,5,6).
qnorm(0.25) = -0.67qnorm(0.75) = 0.67If perfectly normal, expect a straight line!
mydata <- c(2, 3, 3, 4, 5, 5, 6)
quants <- qnorm(c(
0.125, 0.25, 0.375, 0.5,
0.625, 0.75, 0.875))
plot(mydata ~ quants)
The example link is a fun thing to play with! It’s written in Python, so don’t worry about the code. Change the sample sizes to get an idea of just how hard it is to definitively show that data are normal.
Scale-Location
Cook’s Distance
Leverage
Find the values \(\hat\beta_0\), \(\hat\beta_1\), and \(\hat\sigma^2\) that maximize the likelihood of seeing our data.
Under the assumptions that \(X\) is fixed, \(Y = \beta_0 + \beta_1X + \epsilon_i\), and \(\epsilon_i\stackrel{iid}{\sim}N(0,\sigma^2)\), \[ Y \sim N(\beta_0 + \beta_1X, \sigma^2) \]
The probability of observering a data point is: \[ f_Y(y_i|x_i, \beta_0, \beta_1, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-1}{2\sigma^2}(y_i - \beta_0 - \beta_1x_i)^2\right) \]
The likelihood of the parameters, given the data, is: \[ L(\beta_0, \beta_1, \sigma^2|x_i, y_i) = \prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-1}{2\sigma^2}(y_i - \beta_0 - \beta_1x_i)^2\right) \]
Suppose we flipped 10 bottle caps and got 6 “crowns”. Assume the probability of “crown” (\(C\)) is unknown, labelled \(p\).
\[ L(\beta_0, \beta_, \sigma^2) = \prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-1}{2\sigma^2}(y_i - \beta_0 - \beta_1x_i)^2\right) \]
HWK: Show that the estimates for \(\hat\beta_0\) and \(\hat\beta_1\) are the same as the OLS estimates. The estimate for \(\hat\sigma^2\) should come out to: \[ \hat\sigma^2 = \frac{1}{n}\sum_{i=1}^n\left(y_i - \hat\beta_0 - \hat\beta_1x_i\right)^2 \]
The basic question of statistics: “How big is this number?”
Hypothesis test for \(\sigma^2 = \sigma_0^2\) versus \(\sigma^2 > \sigma_0^2\), where \(\sigma_0\) is the value from a previous study.
It can be shown that \(E(MS_{Reg}) = \sigma^2 + \beta_1S_{XX}\).
Consider the Null hypothesis \(\beta_1 =0\) (why is this a good null?).
From a previous class, we know that \[ \frac{(n-2)s^2}{\sigma^2}\sim\chi^2_{n-2} \]
From wikipedia, we know that the mean of a \(\chi^2_k\) distribution is \(k\). Therefore, \[ E\left(\frac{(n-2)s^2}{\sigma^2}\right) = n-2 \Leftrightarrow E(s^2) = \sigma^2 \] and thus \(s^2\) is unbiased.
This does not necessarily mean that \(s^2\) is the best estimator for \(\sigma^2\)!
Even though \(s^2\) is an unbiased estimator, \(s = \sqrt{s^2}\) is biased! Specifically, \(E(s) < \sigma\)
To see why, first note that \[ V(s) = E(s^2) - (E(s))^2 \Leftrightarrow E(s) = \sqrt{E(s^2) - V(s)} \] since \(V(s) > 0\), \(E(s^2) - V(s) < E(s^2)\), and therefore \[ E(s) < \sqrt{E(s^2)} = \sqrt{\sigma^2} = \sigma \]
ln(), then derivative, then set to 0 and solve.Recommended textbook exercises: A, C, E, F, H, I, J, O, P, W
library(aprean3) # install.packages("aprean3")
data(dse03a)
head(dse03a) x y
1 -5 1
2 -4 5
3 -3 4
4 -2 7
5 -1 10
6 0 8# part 1
lm(y ~ x, data = dse03a) |> coef()(Intercept) x
9.272727 1.436364 The prediction equation is \[ \hat y_i = \hat\beta_0 + \hat\beta_1x_i = 9.27 + 1.436 * x_i \]
# part 2
lm(y ~ x, data = dse03a) |> anova() # Reject the null, beta_1 is sig. diff. from 0Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 226.945 226.94 96.18 4.207e-06 ***
Residuals 9 21.236 2.36
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# part 3
e3lm <- lm(y ~ x, data = dse03a)
X <- cbind(1, dse03a$x)
XTX1 <- solve(t(X) %*% X)
betas <- XTX1 %*% t(X) %*% dse03a$y
df <- nrow(dse03a) - 2
preds <- X %*% betas
s <- sqrt(sum((dse03a$y - preds)^2) / df)
XTX <- t(X) %*% X
betas[2] + c(1, -1) * qt(0.025, df = df) * s / sqrt(XTX[2, 2])[1] 1.105045 1.767682confint(e3lm) # check that it's correct 2.5 % 97.5 %
(Intercept) 8.225007 10.320447
x 1.105045 1.767682waiting versus eruptions in the faithful data. Compare the the anova() output, shown below.data(faithful)
y <- faithful$waiting
y_bar <- mean(y)
n <- length(y)
faithful_lm <- lm(waiting ~ eruptions, data = faithful)
y_hat <- predict(faithful_lm)
anova(faithful_lm)Analysis of Variance Table
Response: waiting
Df Sum Sq Mean Sq F value Pr(>F)
eruptions 1 40644 40644 1162.1 < 2.2e-16 ***
Residuals 270 9443 35
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1n <- 100
x <- runif(n, 0, 10)
f_vals <- c()
msreg_vals <-c()
for (i in 1:1000) {
e <- rnorm(n, 0, 1)
y <- -3 + 0.005*x + e
mylm <- lm(y ~ x)
y_hat <- predict(mylm)
y_bar <- mean(y)
ms_reg <- sum((y_hat - y_bar)^2) / 1
msreg_vals[i] <- ms_reg
ms_e <- sum((y - y_hat)^2) / (n - 2)
f_vals <- c(f_vals, ms_reg / ms_e)
}
hist(f_vals, freq = FALSE, breaks = 40)
curve(df(x, 1, n - 2), col = 2, lwd = 3, add = TRUE)
Why is this only true when \(\beta_1=0\)? For intuition, this only should be true when the null is true: recall that p-values are always calculated assuming that the null is true, and the F-stat is no different.
Mathematically, this has to do with the fact that \(E(MS_{Reg}) = \sigma^2 + \beta_1S_{XX}\), not just \(\sigma^2\). Since this is not equal to \(\sigma^2\), this does not follow a Chi-Square distribution.
par(mfrow = c(1, 2))
hist(msreg_vals, breaks = 40, freq = FALSE)
curve(dchisq(x, 1), add = TRUE, col = 2, lwd = 3)
# qqplot, just to check
qqplot(x = qchisq(ppoints(length(msreg_vals)), 1), y = msreg_vals)
abline(a = 0, b = 1)
n <- 100
x <- runif(n, 0, 10)
e <- rnorm(n, 0, 3)
y <- 5 + 2*x + e
mylm <- lm(y ~ x)
plot(mylm$residuals ~ y)
Note that \(y_i = \hat y_i + \hat \epsilon_i\)
\[ cov(\hat\epsilon_i, y_i) = cov(\epsilon_i, \hat y_i + \hat \epsilon_i) = cov(\epsilon_i, \hat y_i) + cov(\hat\epsilon_i, \hat \epsilon_i) = 0 + \sigma^2 \]
Empiricially:
cov(mylm$residuals, y)[1] 7.987392var(mylm$residuals)[1] 7.987392This is why we use residuals versus fitted, not residuals versus observed!
Show that the estimates for \(\hat\beta_0\) and \(\hat\beta_1\) are the same as the OLS estimates.
Find the MLE estimate for \(\sigma^2\).
Derivation left as an exercise. Hint: You are allowed to take a derivative with respect to \(\sigma^2\), not just \(\sigma\).
n <- 100
true_sigma <- 3
n_minus <- seq(-3, 3, 1)
vals <- matrix(ncol = length(n_minus), nrow = 10000)
for (i in 1:10000) {
y <- rnorm(n, 0, true_sigma)
vals[i, ] <- sum((y - mean(y))^2) / (n - n_minus)
}
apply(vals, 2, mean) |> plot(x = n_minus, y = _)
abline(h = true_sigma^2, lwd = 3, col = 2)
From the plot, we get closest to the true value when we subtract 1 from \(n\). However, the MLE subtracts nothing from \(n\)! Why not? Let’s check the MSE.
apply(vals, 2, function(x) {
sum((x - true_sigma^2)^2)
}) |> plot(x = n_minus, y = _)
The MSE is lowest when we divide by… \(n+1\)???? Yes, this is just how it is. See here for a bit more discussion on this issue. Basically, the estimator with \(n-1\) is unbiased, and \(n\) is lower variance without too much bias, while \(n+1\) is much more bias than we’re comfortable with.
ggplot(mtcars) +
aes(x = disp, y = mpg) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, formula = y ~ 1, colour = "darkgrey", linewidth = 2) +
geom_smooth(method = "lm", se = FALSE, formula = y ~ x, colour = "green", linewidth = 2) +
geom_smooth(method = "lm", se = FALSE, formula = y ~ 0 + x, colour = "red", linewidth = 2)
# Model with an intercept (the usual model)
lm_with <- lm(mpg ~ disp, data = mtcars)
# Using "0 +" forces the intercept to be 0
lm_without <- lm(mpg ~ 0 + disp, data = mtcars)anova(lm_with)Analysis of Variance Table
Response: mpg
Df Sum Sq Mean Sq F value Pr(>F)
disp 1 808.89 808.89 76.513 9.38e-10 ***
Residuals 30 317.16 10.57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm_without)Analysis of Variance Table
Response: mpg
Df Sum Sq Mean Sq F value Pr(>F)
disp 1 7599.9 7599.9 36.57 1.073e-06 ***
Residuals 31 6442.4 207.8
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Some interpretations:
lm_with:
lm_without: