After this lecture you should be able to
Recall our non-linear regression model: \[Y = f(X) + \epsilon\]
where \(f\) is an unknown function and \(\epsilon\) is a random error term with mean 0 and variance \(\sigma^2\). We now consider \(X\) a vector of \(p\) explanatory variables, and we can use a multiple linear regression model to approximate \(f\).
d <- Sleuth3::case0901 |>
mutate(Time = case_when(
Time == 1 ~ "Late",
Time == 2 ~ "Early"
))
g <- ggplot(d, aes(x = Intensity, y = Flowers, color = Time, shape = Time)) +
geom_point() +
labs(title = "Sleuth3::case0901")
g
The multiple (linear) regression model is
\[Y_i = \beta_0 + \beta_1 X_{i1} + \cdots + \beta_p X_{ip} + \epsilon_i\]
where \(\epsilon_i\) is a random error term with \(E[\epsilon_i] = 0\) and \(Var[\epsilon_i] = \sigma^2\).
Alternatively the model can be written
\[Y = X\beta + \epsilon\]
where
The error assumptions are now \(E[\epsilon] = 0\) and \(Var[\epsilon] = \sigma^2 I\), where \(I\) is the \(n \times n\) identity matrix.
Throughout we will assume the \(X\) matrix is of full rank, i.e. \(\mathrm{rank}(X) = p+1\), but later we will learn about methods such as ridge regression, lasso and principal component regression that can be used when this assumption is violated.
If you have categorical explanatory variables, you can include them in the model using dummy variables. Suppose you have a categorical variable \(Z\) with \(k\) levels. You can create \(k-1\) dummy variables to represent the levels of \(Z\).
For example, if \(Z\) has levels A, B, and C, let’s choose C as the reference level and create create two dummy variables for A and B:
\[D_A = \mathrm{I}(Z=A) = \left\{ \begin{array}{ll} 1 & \text{if $Z$ = A} \\ 0 & \text{otherwise} \end{array} \right.\] and
\[D_B = \mathrm{I}(Z=B) = \left\{ \begin{array}{ll} 1 & \text{if $Z$ = B} \\ 0 & \text{otherwise} \end{array} \right.\]
Then, you can include these dummy variables in the regression model as follows: \[Y_i = \beta_0 + \beta_1 D_{iA} + \beta_2 D_{iB} + \epsilon_i.\]
The interpretation of the parameters is as follows:
When using polynomial terms, e.g. quadratic, or interactions, the interpretation of the coefficients becomes more complex.
A common transformation is to take the logarithm of the response variable, the explanatory variable, or both. The interpretations of the parameters change accordingly:
Another tool that helps interpretation of models is the use of contrasts. Contrasts are linear combinations of the regression coefficients that allow us to test specific hypotheses about the relationships between the explanatory variables and the response variable.
The vector \(\beta = (\beta_0, \beta_1, \ldots, \beta_p)^\top\) is estimated using the method of least squares:
\[\begin{array}{rl}\hat{\beta} &= \arg\min_{\beta} \sum_{i=1}^n (y_i - \beta_0 - \beta_1 X_{i1} - \cdots - \beta_p X_{ip})^2\\ &= \arg\min_{\beta} (y-X\beta)^\top (y-X\beta) \\ &= \left(X^\top X\right)^{-1} X^\top y \end{array}\]
These are commonly referred to as the ordinary least squares (OLS) estimators.
The residual sum of squares (RSS) are
\[RSS = \left(y-X\hat\beta\right)^\top \left(y-X\hat\beta\right)\]
We can use this to estimate the error variance \(\sigma^2\) using
\[\hat{\sigma}^2 = \frac{RSS}{n-(p+1)}\]
where \(n-(p+1)\) is the degrees of freedom, since we have estimated \((p+1)\) coefficients (\(\beta_0\) through \(\beta_p\)).
m <- lm(Flowers ~ Intensity + Time, data = d) # Additive modelVisualization of residuals:
shared_jitter <- position_jitter(width = 0.5, height = 0, seed = 20260914)
ggplot(m, aes(x = .fitted, y = .fitted + .resid)) +
geom_point(position = shared_jitter) +
geom_segment(
aes(xend = .fitted, yend = .fitted),
position = shared_jitter,
color = "red") +
geom_abline(slope = 1, intercept = 0, color = "blue") +
labs(
title = "Fitted vs Observed",
x = "Predicted value",
y = "Observed value")Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
ℹ Please use `broom::augment(<lm>)` instead.
ℹ The deprecated feature was likely used in the ggplot2 package.
Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
Parameter estimates:
coef(m) (Intercept) Intensity TimeLate
83.46416667 -0.04047143 -12.15833333 summary(m)$sigma[1] 6.441073For inference, we need to assume normality of the errors, i.e. \[\epsilon_i \stackrel{ind}{\sim} N(0, \sigma^2).\] This results in the model \[Y_i \stackrel{ind}{\sim} N(\beta_0 + \beta_1 X_{i1}+\cdots+\beta_pX_{ip}, \sigma^2)\] or, equivalently,
\[Y \sim N_n(X\beta, \sigma^2 I)\] where \(N_n\) is the \(n\)-dimensional multivariate normal distribution.
With the normality assumption, we can write the likelihood function:
\[\begin{array}{rl} L(\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_1 X_{i1} - \cdots - \beta_pX_{ip}]^2\right) \\ &= (2\pi\sigma^2)^{-n/2} \exp\left(-\frac{1}{2\sigma^2} [y-X\beta]^\top [y-X\beta] \right) \end{array}\]
the estimates for \(\beta\) that maximize this likelihood are the same as those obtained from least squares. These are referred to as the maximum likelihood estimators (MLEs) of \(\beta\).
For testing hypotheses about the regression coefficients, we use t-tests. The test statistic for \(\beta_j\) is
\[t = \frac{\hat{\beta}_j - \beta_j^0}{\text{SE}(\hat{\beta}_j)}\]
which, under the null hypothesis \(H_0: \beta_j = \beta_j^0\), follows a t-distribution with \(n-(p+1)\) degrees of freedom.
summary(m)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 83.46416667 3.27377202 25.494801 2.784695e-17
Intensity -0.04047143 0.00513237 -7.885525 1.036787e-07
TimeLate -12.15833333 2.62955696 -4.623719 1.463776e-04A \((1-\alpha)100\%\) confidence interval for \(\beta_j\) is given by
\[\hat{\beta}_j \pm t_{n-(p+1), 1-\alpha/2} \cdot \text{SE}\left(\hat{\beta}_j\right)\]
where \(t_{n-(p+1), 1-\alpha/2}\) is the \((1-\alpha/2)\) quantile of the t-distribution with \(n-(p+1)\) degrees of freedom.
The standard error of \(\hat{\beta}_j\) is given by \[\text{SE}\left(\hat{\beta}_j\right) = \sqrt{\hat{\sigma}^2 \left(X^\top X\right)^{-1}_{jj}}\] where \(\left(X^\top X\right)^{-1}_{jj}\) is the \(j\)-th diagonal element of the inverse of the matrix \(X^\top X\).
confint(m, level = 0.95) 2.5 % 97.5 %
(Intercept) 76.65598505 90.27234828
Intensity -0.05114478 -0.02979808
TimeLate -17.62679640 -6.68987027Contrasts also utilize t-tools, e.g. t-tests and confidence intervals.
A \((1-\alpha)100\%\) confidence interval for the mean response at a specific value \(x_0\) is given by
\[\hat{f}(x_0) \pm t_{n-(p+1), 1-\alpha/2} \cdot \text{SE}\left(\hat{f}(x_0)\right)\]
where \(x_0\) is now a vector of feature values, \(\hat{f}(x_0) = x_0^\top\hat\beta = \hat{\beta}_0 + \hat{\beta}_1 x_{01} + \cdots + \hat\beta_p x_{0p}\), and \(\text{SE}(\hat{f}(x_0))\) is the standard error of the predicted mean response. This is the uncertainty in the line at a given point \(x_0\).
nd <- unique(d |> select(Intensity, Time))
cx <- bind_cols(nd,
predict(
m,
newdata = unique(d |> select(Intensity, Time)),
interval = "confidence",
level = 0.95)) |>
rename(
Flowers = fit
)
cx Intensity Time Flowers lwr upr
1 150 Late 65.23512 59.66986 70.80038
2 300 Late 59.16440 54.61256 63.71625
3 450 Late 53.09369 49.14491 57.04247
4 600 Late 47.02298 43.07420 50.97175
5 750 Late 40.95226 36.40042 45.50410
6 900 Late 34.88155 29.31629 40.44681
7 150 Early 77.39345 71.82819 82.95871
8 300 Early 71.32274 66.77090 75.87458
9 450 Early 65.25202 61.30325 69.20080
10 600 Early 59.18131 55.23253 63.13009
11 750 Early 53.11060 48.55875 57.66244
12 900 Early 47.03988 41.47462 52.60514g +
geom_ribbon(data = cx, aes(ymin = lwr, ymax = upr, fill = Time), color = NA, alpha = 0.2) +
geom_line(data = cx, aes(linetype = Time))
A \((1-\alpha)100\%\) prediction interval for a new observation at a specific value \(x_0\) is given by \[\hat{Y}_0 \pm t_{n-(p+1), 1-\alpha/2} \cdot \text{SE}_{\text{pred}}(\hat{Y}_0)\]
where \(\hat{Y}_0 = \hat{f}(x_0) = x_0^\top\hat\beta\) and \(\text{SE}_{\text{pred}}(\hat{Y}_0)\) is the standard error of the predicted value, which accounts for both the uncertainty in the mean response and the variability of individual observations. The key here is that the prediction interval is wider than the confidence interval for the mean response, reflecting the additional uncertainty in predicting a single new observation due to the residual error \(\epsilon\).
px <- bind_cols(nd,
predict(
m,
newdata = unique(d |> select(Intensity, Time)),
interval = "prediction", # THIS CHANGED
level = 0.95)) |>
rename(
Flowers = fit
)
px Intensity Time Flowers lwr upr
1 150 Late 65.23512 50.73006 79.74017
2 300 Late 59.16440 45.01719 73.31162
3 450 Late 53.09369 39.12883 67.05855
4 600 Late 47.02298 33.05811 60.98784
5 750 Late 40.95226 26.80504 55.09948
6 900 Late 34.88155 20.37649 49.38660
7 150 Early 77.39345 62.88840 91.89851
8 300 Early 71.32274 57.17552 85.46996
9 450 Early 65.25202 51.28716 79.21689
10 600 Early 59.18131 45.21645 73.14617
11 750 Early 53.11060 38.96338 67.25781
12 900 Early 47.03988 32.53483 61.54494g +
geom_ribbon(data = px, aes(ymin = lwr, ymax = upr, fill = Time), color = NA, alpha = 0.2) +
geom_line(data = px, aes(linetype = Time))
To evaluate the statistical evidence for a group of explanatory variables, e.g. a single categorical variable with \(k+1\) levels, we can use an F-test. The null hypothesis is that all coefficients for the group of variables are equal to zero, i.e \[H_0: \beta_{j_1} = \beta_{j_2} = \cdots = \beta_{j_k} = 0\] where \(j_1, j_2, \ldots, j_k\) are the indices of the coefficients for the group of variables. The alternative hypothesis is that at least one of the coefficients is not equal to zero. The F-statistic is calculated as \[F_{k, n-(p+1)} = \frac{(RSS_r - RSS_f)/k}{RSS_f/(n-[p+1])} = \frac{(RSS_r - RSS_f)/k}{\hat\sigma^2_{f}}\]
where
Under the null hypothesis, the F-statistic follows an F-distribution with \(k\) and \(n-(p+1)\) degrees of freedom. A \(p\)-value can be calculated from the F-distribution to assess the statistical significance of the group of variables.
# Default
anova(m)Analysis of Variance Table
Response: Flowers
Df Sum Sq Mean Sq F value Pr(>F)
Intensity 1 2579.75 2579.75 62.181 1.037e-07 ***
Time 1 886.95 886.95 21.379 0.0001464 ***
Residuals 21 871.24 41.49
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(m, test="F")Single term deletions
Model:
Flowers ~ Intensity + Time
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 871.2 92.205
Intensity 1 2579.75 3451.0 123.241 62.181 1.037e-07 ***
Time 1 886.95 1758.2 107.056 21.379 0.0001464 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tail(capture.output(summary(m)), n = 2)[1][1] "F-statistic: 41.78 on 2 and 21 DF, p-value: 4.786e-08"# Manual (for nested models)
m0 <- lm(Flowers ~ 1, data = d) # Intercept only model
anova(m0, m)Analysis of Variance Table
Model 1: Flowers ~ 1
Model 2: Flowers ~ Intensity + Time
Res.Df RSS Df Sum of Sq F Pr(>F)
1 23 4337.9
2 21 871.2 2 3466.7 41.78 4.786e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The coefficient of determination, \(R^2\), measures the proportion of the variance in the response variable that is predictable from the explanatory variable(s). It is defined as \[R^2 = 1 - \frac{RSS}{TSS}\]
where \(RSS\) is the residual sum of squares and \(TSS\) is the total sum of squares, i.e.
\[TSS = \sum_{i=1}^n (y_i - \bar{y})^2.\]
Thus, \(R^2\) can be interpreted as the proportion of the total variability in the response variable that is explained by the model.
summary(m)$r.squared[1] 0.7991589Default plot:
g + geom_smooth(method = "lm", formula = y ~ x) # lines are not parallel 
Model summary:
summary(m)
Call:
lm(formula = Flowers ~ Intensity + Time, data = d)
Residuals:
Min 1Q Median 3Q Max
-9.652 -4.139 -1.558 5.632 12.165
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 83.464167 3.273772 25.495 < 2e-16 ***
Intensity -0.040471 0.005132 -7.886 1.04e-07 ***
TimeLate -12.158333 2.629557 -4.624 0.000146 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.441 on 21 degrees of freedom
Multiple R-squared: 0.7992, Adjusted R-squared: 0.78
F-statistic: 41.78 on 2 and 21 DF, p-value: 4.786e-08