Functions from these packages will be used throughout this document:
[R code]
library(conflicted) # check for conflicting function definitions# printr is not loaded here: it inserts help-file output into markdown outputlibrary(rmarkdown) # Convert R Markdown documents into a variety of formats.library(pander) # format tables for markdownlibrary(ggplot2) # graphicslibrary(ggfortify) # help with graphicslibrary(dplyr) # manipulate datalibrary(tibble) # `tibble`s extend `data.frame`slibrary(magrittr) # `%>%` and other additional piping toolslibrary(haven) # import Stata fileslibrary(knitr) # format R output for markdownlibrary(tidyr) # Tools to help to create tidy datalibrary(plotly) # interactive graphicslibrary(dobson) # datasets from Dobson and Barnett 2018library(parameters) # format model output tables for markdownlibrary(haven) # import Stata fileslibrary(latex2exp) # use LaTeX in R code (for figures and tables)library(fs) # filesystem path manipulationslibrary(survival) # survival analysislibrary(survminer) # survival analysis graphicslibrary(KMsurv) # datasets from Klein and Moeschbergerlibrary(parameters) # format model output tables forlibrary(webshot2) # convert interactive content to static for pdflibrary(forcats) # functions for categorical variables ("factors")library(stringr) # functions for dealing with stringslibrary(lubridate) # functions for dealing with dates and times# Summarizes key information about statistical objects in tidy tibbleslibrary(broom)# Provides suite of functions to work with regression model# `broom::tidy()` tibbleslibrary(broom.helpers)
Here are some R settings I use in this document:
[R code]
rm(list =ls()) # delete any data that's already loaded into Rconflicts_prefer(dplyr::filter)ggplot2::theme_set( ggplot2::theme_bw() +# ggplot2::labs(col = "") + ggplot2::theme(legend.position ="bottom",text = ggplot2::element_text(size =12, family ="serif") ))knitr::opts_chunk$set(message =FALSE)options(digits =6)panderOptions("big.mark", ",")pander::panderOptions("table.emphasize.rownames", FALSE)pander::panderOptions("table.split.table", Inf)# use the `filter()` function from dplyr() by defaultconflicts_prefer(dplyr::filter)legend_text_size <-9run_graphs <-TRUE
Richard McElreath’s Geocentric Models (Statistical Rethinking 2026, Lecture A03) accompanies Chapter 4 of McElreath [12], which introduces Gaussian linear regression. The lecture fits these models using Bayesian inference rather than the least-squares and maximum likelihood estimation used in this chapter [13].
UC Davis STA 108 (“Applied Statistical Methods: Regression Analysis”) is a prerequisite for this course, so everyone here should have some understanding of linear regression already.
We will review linear regression to:
make sure everyone is caught up
provide an epidemiological perspective on model interpretation.
Here \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\) is the observed covariate vector for unit \(i\). Here \(\ \sim_{\perp\!\!\!\perp}\ \) means independently distributed. Conditional on \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_1,\ldots,\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_n\), the responses are independent across units, and each \(Y_i|\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\) follows the same conditional distributional form with common variance\(\sigma^2\) while allowing \(\mu_i\) to vary by observation.
2.2 Motivating example: birthweights and gestational age
Suppose we want to learn about the distributions of birthweights (outcome\(Y\)) for (human) babies born at different gestational ages (covariate\(A\)) and with different chromosomal sexes (covariate\(S\)) (Dobson and Barnett [1] Example 2.2.2).
pred_female <-coef(bw_lm1)["(Intercept)"] +coef(bw_lm1)["age"] *36## or using built-in prediction:pred_female_alt <-predict(bw_lm1, newdata =tibble(sex ="female", age =36))
Here’s another way we could rewrite this model (by collecting terms involving \(S\)):
\[
E[Y|M,A] = \beta_0 + \beta_M M+ (\beta_A + \beta_{AM} M) A
\]
Remark 1 (Interactions and the number of subgroups). If every subgroup has its own line, then each line is estimated from less data: the more subgroups we compare, the fewer observations each one has. When there are too many subgroups, we restrict how they may differ. The parallel-lines model without an interaction (Equation 2 with \(\beta_{AM} = 0\)) is one such restriction: the subgroups share one slope, and only their intercepts differ. Adding the interaction term relaxes the restriction, so that the slopes can also differ.
This relaxation is symmetric. If the slope of birthweight with respect to age differs by sex, then the difference in mean birthweight between the sexes differs by age, which is visible in Figure 3 as vertical gaps that are narrower at some ages than at others.
The same idea applies to two numeric covariates. The slope of the outcome with respect to one numeric covariate \(A\) may depend on the value of another numeric covariate \(Z\). Then each value of \(Z\) defines its own subgroup and its own slope, so there are as many subgroups as there are distinct values of \(Z\), and most of them contain no observations at all. To make progress, we have to assume a restrictive shape for how the slope changes with \(Z\); a product term \(\beta_{AZ}\, a \cdot z\) assumes that the slope changes linearly with \(Z\).
Exercise 7 According to this model, what’s the mean birthweight for a female born at 36 weeks?
Exercise 10 What is the interpretation of \(\beta_{M}\) in Model 2?
Solution.
Mean birthweight among males with gestational age 0 weeks: \[
\begin{aligned}
\mu(1,0) &= {\mathbb{E}\mathopen{}\left[Y|M = 1,A = 0\right]\mathclose{}}\\
&= \beta_0 + {\color{red}{\beta_M}} \cdot 1 + \beta_A \cdot 0 + \beta_{AM}\cdot 1 \cdot 0\\
&= \beta_0 + {\color{red}{\beta_M}}
\end{aligned}
\] Mean birthweight among females with gestational age 0 weeks: \[
\begin{aligned}
\mu(0,0) &= {\mathbb{E}\mathopen{}\left[Y|M = 0,A = 0\right]\mathclose{}}\\
&= \beta_0 + {\color{red}{\beta_M}} \cdot 0 + \beta_A \cdot 0 + \beta_{AM}\cdot 0 \cdot 0\\
&= \beta_0
\end{aligned}
\]
\[
\begin{aligned}
\beta_{M} &= \mu(1,0) - \mu(0,0) \\
&= {\mathbb{E}\mathopen{}\left[Y|M = 1,A = 0\right]\mathclose{}} - {\mathbb{E}\mathopen{}\left[Y|M = 0,A = 0\right]\mathclose{}}
\end{aligned}
\]\(\beta_M\) is the difference in mean birthweight between males with gestational age 0 weeks and females with gestational age 0 weeks.
Exercise 11 What is the interpretation of \(\beta_{AM}\) in Model 2?
Definition 2 (General procedure for interpreting a regression coefficient) To derive the interpretation of coefficient \(\beta_j\) in a model with linear predictor\(\eta(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) = \eta(x_1, \ldots, x_p)\):
Write down the linear predictor\(\eta(x_1, \ldots, x_p)\).
Differentiate or difference\(\eta\) with respect to \(x_j\):
If \(x_j\) is continuous: compute the partial derivative\(\frac{\partial}{\partial x_j}\eta(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})\).
If \(x_j\) is discrete: compute the one-unit difference \(\left.{\eta(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})}\right|_{x_j = a+1} - \left.{\eta(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})}\right|_{x_j = a}\).
For an interaction coefficient\(\beta_{jk}\) (multiplying the product \(x_j \cdot x_k\)), apply this step twice in sequence: first differentiate or difference \(\eta\) with respect to \(x_j\) to obtain the slope of \(x_j\) as a function of \(x_k\); then differentiate or difference that slope with respect to \(x_k\) (choosing continuous/discrete as appropriate for \(x_k\)) to obtain the change in slopes, which equals \(\beta_{jk}\).
If any covariates remain, set them to their reference levels. If covariates other than the one(s) used in step 2 still appear in the resulting expression, substitute each covariate’s reference level for it. This substitution eliminates the corresponding interaction terms and isolates \(\beta_j\) alone on the right-hand side.
Identify \(\beta_j\): the expression from the previous steps equals \(\beta_j\); write \(\beta_j = \cdots\) to state the interpretation explicitly.
Special case — intercept (\(\beta_0\)): If the coefficient of interest is the intercept \(\beta_0\), skip step 2. Apply step 3 directly by setting all covariates to their reference levels: \(\beta_0 = \eta(\text{reference pattern})\).
Example 1 (Applying the procedure to Model 1) For the birthweight parallel-lines model \[{\mathbb{E}\mathopen{}\left[Y|M=m, A=a\right]\mathclose{}} = \beta_0 + \beta_M m + \beta_A a,\] where \(M \in \{0,1\}\) (sex; \(1\) = male, \(0\) = female) and \(A\) is gestational age (weeks):
\(\beta_A\) (gestational age, continuous):\(\frac{\partial}{\partial a}\eta(m,a) = \beta_A\), so \(\beta_A = \frac{\partial}{\partial a}{\mathbb{E}\mathopen{}\left[Y|M=m, A=a\right]\mathclose{}}\). (Step 3 not needed: no remaining covariates appear after differentiating.)
so \(\beta_M = {\mathbb{E}\mathopen{}\left[Y|M=1, A=a\right]\mathclose{}} - {\mathbb{E}\mathopen{}\left[Y|M=0, A=a\right]\mathclose{}}\). (Step 3 not needed: \(a\) cancels in the differencing, so the result holds for all values of \(A\).)
\(\beta_0\) (intercept, special case): Skip step 2; set \(m=0\) and \(a=0\) (the reference levels) in step 3: \(\beta_0 = \eta(0,0) = {\mathbb{E}\mathopen{}\left[Y|M=0, A=0\right]\mathclose{}}\).
2.6 Stratified regression
\[
{\mathbb{E}\mathopen{}\left[Y|A=a, S=s\right]\mathclose{}} =
\beta_M m + \beta_{AM} (a \cdot m) +
\beta_F f + \beta_{AF} (a \cdot f)
\tag{3}\]
Compare this stratified model (Equation 3) with our interaction model, Equation 2:
\[
{\mathbb{E}\mathopen{}\left[Y|A=a, S=s\right]\mathclose{}} =
\beta_0 + \beta_A a + \beta_M m + \beta_{AM} (a \cdot m)
\]
[R code]
bw_lm2 <-lm(weight ~ sex + age + sex:age, data = bw)bw_lm2 |>parameters() |>print_md(include_reference = include_reference_lines,select ="{estimate}" )
Let the coefficients of this model be \(\gamma\)s instead of \(\beta\)s.
What are the interpretations of the \(\gamma\)s? How do they relate to the \(\beta\)s in Model 2? Which have the same interpretation? Which are different, and how do they differ? What is the pattern?
Solution. Interpretation of \(\gamma_0\):
From Model 4, \(\gamma_0\) is the mean birthweight among females (\(M = 0\)) with \(A^* = 0\) (i.e., \(A = 32\) weeks):
This centered intercept parameter differs from \(\beta_0 = {\mathbb{E}\mathopen{}\left[Y|M=0, A=0\right]\mathclose{}}\), which is the mean birthweight among females at \(A = 0\) weeks.
Interpretation of \(\gamma_M\):
From Model 4, \(\gamma_M\) is the sex difference in mean birthweight at \(A^* = 0\) (i.e., \(A = 32\) weeks):
Since shifting \(A\) by a constant does not change the slope, \(\gamma_{A^*} = \beta_A\): these two coefficients have the same value and interpretation.
Interpretation of \(\gamma_{A^*M}\):
From Model 4, \(\gamma_{A^*M}\) is the difference in slope with respect to \(A^*\) between males and females:
Since shifting \(A\) by a constant does not change slopes, \(\gamma_{A^*M} = \beta_{AM}\): these two coefficients have the same value and interpretation.
The pattern:
Slope coefficients (\(\gamma_{A^*}\) and \(\gamma_{A^*M}\)) are unchanged by rescaling: they have the same values and interpretations as the corresponding \(\beta\)s.
Coefficients change only for variables that have interactions with the rescaled variable \(A\). This set of modified coefficients includes the intercept (which can be viewed as the main effect of a variable that interacts with \(A\) via \(\beta_A\)), and the main effect of \(M\) (which interacts with \(A\) via \(\beta_{AM}\)). Shifting \(A\) by 32 weeks changes the reference point from \(A = 0\) to \(A = 32\), so these coefficients now represent quantities evaluated at \(A = 32\) weeks rather than at \(A = 0\) weeks.
Exercise 13 Using R, fit the rescaled interaction model with \(A^* = A - 36\) weeks in place of \(A\) in Model 2. Compare the coefficient estimates with those from the original model. Which coefficients change, and which remain the same?
Solution.
[R code]
bw <- bw |>mutate(`age - mean`= age -mean(age),`age - 36wks`= age -36 )lm1_c <-lm(weight ~ sex +`age - 36wks`, data = bw)lm2_c <-lm(weight ~ sex +`age - 36wks`+ sex:`age - 36wks`, data = bw)parameters(lm2_c, ci_method ="wald") |>print_md()
Centering gestational age does not change predictions
Centering gestational age changes the coefficient parameterization, but it does not change fitted values, confidence bands, or prediction bands.
The next output reports maximum absolute differences between the uncentered and centered models. All values should be near zero (up to floating-point rounding), which confirms that centering changes parameterization only.
Figure 6: Two-panel comparison of the fitted birthweight model bw_lm2 with and without centering gestational age, showing identical fitted values, confidence bands, and prediction bands in both panels.
2.9 Categorical covariates with more than two levels
Example: birthweight
In the birthweight example, the variable sex had only two observed values:
[R code]
unique(bw$sex)#> [1] female male #> Levels: female male
If there are more than two observed values, we can’t just use a single variable with 0s and 1s.
If we want to model Sepal.Length by species, we could create a variable \(X\) that represents “setosa” as \(X=1\), “virginica” as \(X=2\), and “versicolor” as \(X=3\).
[R code]
data(iris) # this step is not always necessary, but ensures you're starting# from the original version of a dataset stored in a loaded packageiris <- iris |>tibble() |>mutate(X =case_when( Species =="setosa"~1, Species =="virginica"~2, Species =="versicolor"~3 ) )iris |>distinct(Species, X)
Table 10: iris data with numeric coding of species
Then we could fit a model like:
[R code]
iris_lm1 <-lm(Sepal.Length ~ X, data = iris)iris_lm1 |>parameters() |>print_md()
Table 11: Model of iris data with numeric coding of Species
Example 2 The hers data [3] is available in the rmb package:
[R code]
hers <- rmb::hers |> haven::zap_labels()
[R code]
hers |>head()
Table 16: HERS data
The physact variable records each participant’s self-reported physical activity, relative to other women her age, on an ordered five-level scale (1 = much less active, …, 5 = much more active). We make it an ordered factor so that lm() uses polynomial contrasts, which respect the ordering, instead of the default treatment contrasts:
Table 17: Mean LDL by physical-activity level, using polynomial contrasts
Parameter
Coefficient
SE
95% CI
t(2747)
p
(Intercept)
144.71
0.85
(143.05, 146.38)
170.10
< .001
physact ord (linear)
-0.31
2.29
(-4.80, 4.18)
-0.14
0.892
physact ord (quadratic)
-1.86
2.05
(-5.88, 2.15)
-0.91
0.363
physact ord (cubic)
0.47
1.74
(-2.94, 3.89)
0.27
0.785
physact ord (4th degree)
-0.58
1.42
(-3.36, 2.21)
-0.40
0.686
With polynomial contrasts, each coefficient summarizes a trend component across the ordered levels: the linear term measures the strength and direction of a straight-line (monotone) trend in mean LDL across the levels, and the quadratic and higher terms capture curvature. The contrasts are orthonormal, so the linear coefficient is scaled by the contrast normalization; it is not the raw change in mean LDL per level. Here the linear-trend coefficient is small and not statistically significant, so these data show little evidence of a monotone association between physical-activity level and mean LDL cholesterol.
To find the maximum likelihood estimator (MLE), we set the score vector equal to the zero vector. The score equation for \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) is:
Assume the design matrix\(\mathbf{X}\) has full column rank\(p+1\). This rank condition implies that \(\sum_{i=1}^n\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i {\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i}^{\top} = {\mathbf{X}}^{\top}\mathbf{X}\) is an invertible\((p+1) \times (p+1)\) matrix. To solve for \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\), use the identity \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i (\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}) = (\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i {\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i}^{\top})\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\):
Differentiating each term in the expanded sum-of-squared-deviations quadratic form with respect to \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\):
Multiplying the sum-of-squared-deviations gradient by the scalar factor \(-\frac{1}{2\sigma^2}\) yields the full log-likelihood score vector \(\ell'_{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\):
The Hessian matrix \(-\frac{1}{\sigma^2}{\mathbf{X}}^{\top}\mathbf{X}\) is negative definite when \(\mathbf{X}\) has full column rank \(p+1\), confirming that \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}= ({\mathbf{X}}^{\top}\mathbf{X})^{-1}{\mathbf{X}}^{\top}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}\) uniquely maximizes the log-likelihood and is the maximum likelihood estimator for \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\).
summary(bw_lm2)#> #> Call:#> lm(formula = weight ~ sex + age + sex:age, data = bw)#> #> Residuals:#> Min 1Q Median 3Q Max #> -246.7 -138.1 -39.1 176.6 274.3 #> #> Coefficients:#> Estimate Std. Error t value Pr(>|t|) #> (Intercept) -2141.7 1163.6 -1.84 0.08057 . #> sexmale 873.0 1611.3 0.54 0.59395 #> age 130.4 30.0 4.35 0.00031 ***#> sexmale:age -18.4 41.8 -0.44 0.66389 #> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#> #> Residual standard error: 181 on 20 degrees of freedom#> Multiple R-squared: 0.643, Adjusted R-squared: 0.59 #> F-statistic: 12 on 3 and 20 DF, p-value: 0.000101
3.8 Predicted and fitted values
Definition 3 (Predicted value) In a regression model\(\operatorname{p}(y|\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})\), the predicted value of \(y\) given \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) is the estimated mean of \(Y\) given \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}=\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\):
\[\hat y \stackrel{\text{def}}{=}{\mathop{\hat{\mathbb{E}}}\nolimits\mathopen{}\left[Y|\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}=\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right]\mathclose{}}\]
For linear models, the predicted value can be straightforwardly calculated by multiplying each predictor value \(x_j\) by its corresponding coefficient \(\beta_j\) and adding up the results:
Definition 4 (Fitted value) For a given dataset \((\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}, \mathbf{X})\) and corresponding fitted model \(\operatorname{p}_{\hat{\beta}}(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}|\mathbf{x})\), the fitted value of \(y_i\) is the predicted value (see Definition 3) of \(y\) when \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}=\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\) using the estimated parameters \(\hat{\beta}\):
Remark 2 (Fitted values versus predicted values). Fitted values and predicted values use the same formula, but they answer different questions:
A fitted value\(\hat y_i\) describes an observation that is already in the dataset. Its outcome \(y_i\) helped determine \(\hat{\beta}\), so the model has already “seen” \(y_i\).
A predicted value\(\hat y\) describes a new or hypothetical observation with covariate vector \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\), whose outcome was not used to estimate \(\hat{\beta}\) (and \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) need not match any \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\) in the data).
Because each observation pulls the fitted line toward itself, fitted values tend to be closer to their observed outcomes than predicted values are to the outcomes of new observations. So the residuals\(y_i - \hat y_i\) tend to understate how far off the model’s predictions will be for new data.
Example 3 (Fitted versus predicted value for one birthweight observation) Take the first observation in the birthweight data. Its fitted value comes from bw_lm1, which was estimated using that observation. To predict the same observation as if it were new, refit the model without it and predict at its covariate values:
where \(\ell(\hat\theta)\) is the log-likelihood evaluated at the maximum-likelihood estimates \(\hat\theta\), \(p\) is the number of estimated parameters (including \(\hat\sigma^2\) for Gaussian models), and \(n\) is the number of observations.
Definition 6 (Bayesian Information Criterion (BIC))\[\text{BIC} = -2 \ell(\hat\theta) + p \log(n)\]
where \(\ell(\hat\theta)\), \(p\), and \(n\) are defined as in Definition 5.
Conceptual basis
Each criterion has two components:
Fit term (\(-2\ell(\hat\theta)\)): measures lack of fit — lower is better. A model with more parameters will always achieve a higher (or equal) log-likelihood on the observed data.
Penalty term (\(2 p\) for AIC; \(p \log(n)\) for BIC): penalizes model complexity to guard against overfitting.
Together, they balance goodness of fit against parsimony.
AIC vs. BIC
Criterion
Penalty per parameter
Tends to select
AIC
\(2\)
larger models
BIC
\(\log(n)\)
smaller models
The BIC penalty exceeds the AIC penalty whenever \(\log(n) > 2\), i.e., when \(n > e^2 \approx 7.4\). In practice, BIC almost always penalizes additional parameters more heavily than AIC and therefore tends to select simpler (more parsimonious) models [3], [7].
AIC in R
[R code]
-2*logLik(bw_lm2) |>as.numeric() +2* (length(coef(bw_lm2)) +1) # sigma counts as a parameter here#> [1] 323.159AIC(bw_lm2)#> [1] 323.159
Lower values are better. There are no hypothesis tests or p-values associated with these criteria.
To compare models, calculate the criterion for each model and choose the model with the smallest value.
Differences of less than 2 units are generally considered negligible; differences greater than 10 are considered strong evidence in favor of the lower-criterion model.
BIC tends to favor simpler models than AIC, especially when the sample size is large.
(Residual) Deviance
Let \(q\) be the number of distinct covariate combinations in a data set.
For example, in the birthweight data, there are \(q = 12\) unique patterns (Table 18).
[R code]
bw_X_unique
Table 18: Unique covariate combinations in the birthweight data, with replicate counts
Definition 7 (Replicates) If a given covariate pattern has more than one observation in a dataset, those observations are called replicates.
Example 4 (Replicates in the birthweight data) In the birthweight dataset, there are 2 replicates of the combination “female, age 36” (Table 18).
Exercise 14 (Replicates in the birthweight data) Which covariate pattern(s) in the birthweight data has the most replicates?
Solution 1 (Replicates in the birthweight data). Two covariate patterns are tied for most replicates: males at age 40 weeks and females at age 40 weeks. 40 weeks is the usual length for human pregnancy (Polin et al. [19]), so this result makes sense.
[R code]
bw_X_unique |> dplyr::filter(n ==max(n))
Saturated models
The most complicated model we could fit would have one parameter (a mean) for each covariate pattern, plus a variance parameter. For the birthweight data, each covariate pattern is a combination \((s, a)\) of sex \(S\) and gestational age \(A\) that appears in the data, so this model has the structure:
Table 19: Saturated model for the birthweight data
Parameter
Coefficient
SE
95% CI
t(12)
p
sex (male) × age35
2925.00
187.92
(2515.55, 3334.45)
15.56
< .001
sex (female) × age36
2570.50
132.88
(2280.98, 2860.02)
19.34
< .001
sex (male) × age36
2625.00
187.92
(2215.55, 3034.45)
13.97
< .001
sex (female) × age37
2539.00
187.92
(2129.55, 2948.45)
13.51
< .001
sex (male) × age37
2737.50
132.88
(2447.98, 3027.02)
20.60
< .001
sex (female) × age38
2872.50
132.88
(2582.98, 3162.02)
21.62
< .001
sex (male) × age38
2982.00
108.50
(2745.60, 3218.40)
27.48
< .001
sex (female) × age39
2846.00
132.88
(2556.48, 3135.52)
21.42
< .001
sex (female) × age40
3152.25
93.96
(2947.52, 3356.98)
33.55
< .001
sex (male) × age40
3256.25
93.96
(3051.52, 3460.98)
34.66
< .001
sex (male) × age41
3292.00
187.92
(2882.55, 3701.45)
17.52
< .001
sex (female) × age42
3210.00
187.92
(2800.55, 3619.45)
17.08
< .001
We call this model the full, maximal, or saturated model for this dataset.
Definition 8 (Saturated model (full model, maximal model)) Consider a dataset with \(q\) distinct covariate patterns \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_1, \ldots, \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_q\), and an outcome distribution family (for example, Gaussian, Bernoulli, or Poisson). The saturated model for that dataset and family (also called the full model or maximal model) gives each distinct covariate pattern its own free mean parameter:
\[
{\mathbb{E}\mathopen{}\left[Y \mid \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_k\right]\mathclose{}} \stackrel{\text{def}}{=}\mu_k,
\quad k = 1, \ldots, q,
\]
plus any dispersion parameter of the family (such as \(\sigma^2\) for the Gaussian family), shared by all patterns [1, Ch. 5].
Example 5 (Saturated model for a small dataset) Suppose we observe \(n = 6\) outcomes at \(q = 3\) covariate patterns \(x \in \{0, 1, 2\}\), with two replicates (Definition 7) per pattern:
\(x\)
\(y\) values
pattern mean \(\bar{y}_x\)
0
1, 3
\((1 + 3)/2 = 2\)
1
5, 7
\((5 + 7)/2 = 6\)
2
6, 8
\((6 + 8)/2 = 7\)
The Gaussian saturated model for these data has three mean parameters, \(\mu_0\), \(\mu_1\), and \(\mu_2\), plus one variance parameter \(\sigma^2\). Its maximum likelihood estimates of the means are the pattern means: \(\hat{\mu}_0 = 2\), \(\hat{\mu}_1 = 6\), and \(\hat{\mu}_2 = 7\).
Example 6 (A model that is not saturated (counterexample)) For the same data as Example 5, the simple linear regression model \({\mathbb{E}\mathopen{}\left[Y \mid X = x\right]\mathclose{}} = \beta_0 + \beta_1 x\) has only \(p = 2\) mean parameters, fewer than the \(q = 3\) covariate patterns. Because the three pattern means \(2\), \(6\), and \(7\) do not lie on a straight line (the differences between consecutive means are \(6 - 2 = 4\) and \(7 - 6 = 1\), which are unequal), no choice of \(\beta_0\) and \(\beta_1\) reproduces all three pattern means, so this model is not saturated.
Figure 9: Model 2 and saturated model for birthweight data, with confidence and prediction intervals
We can calculate the log-likelihood of this model as usual:
[R code]
logLik(lm_max)#> 'log Lik.' -151.402 (df=13)
We can compare this model to our other models using chi-square tests, as usual:
[R code]
library(lmtest)lrtest(lm_max, bw_lm2)
The likelihood ratio statistic for this test is \[\lambda = 2 * (\ell_{\text{full}} - \ell) = 10.355374\] where:
\(\ell_{\text{full}}\) is the log-likelihood of the full model: -151.401601
\(\ell\) is the log-likelihood of our comparison model (two slopes, two intercepts): -156.579288
This statistic is called the deviance or residual deviance for our two-slopes and two-intercepts model; it tells us how much the likelihood of that model deviates from the likelihood of the maximal model.
The corresponding p-value tells us whether there we have enough evidence to detect that our two-slopes, two-intercepts model is a worse fit for the data than the maximal model; in other words, it tells us if there’s evidence that we missed any important patterns. (Remember, a nonsignificant p-value could mean that we didn’t miss anything and a more complicated model is unnecessary, or it could mean we just don’t have enough data to tell the difference between these models.)
Definition 9 (Deviance and residual deviance) Let \(\ell_{\text{sat}}\) be the maximized log-likelihood of the saturated model (Definition 8), and let \(\ell(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}})\) be the maximized log-likelihood of a model \(M\) with the same outcome distribution family. The deviance of model \(M\), also called its residual deviance, is:
\[
D \stackrel{\text{def}}{=}2\mathopen{}\left(\ell_{\text{sat}} - \ell(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}})\right)\mathclose{}
\]
When the family has a dispersion parameter (such as \(\sigma^2\) in a Gaussian model), both log-likelihoods are evaluated at the same value of that dispersion parameter.
Definition 10 (Null model) The null model is the model with a single mean parameter shared by all observations (\({\mathbb{E}\mathopen{}\left[Y \mid \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right]\mathclose{}} = \beta_0\), an intercept only).
Definition 11 (Null deviance) The null deviance\(D_0\) is the deviance (Definition 9) of the null model (Definition 10):
where \(\ell_0\) is the maximized log-likelihood of the null model [1, Ch. 5].
Example 7 (Residual and null deviance for a small dataset) Use the data from Example 5, and suppose that the variance \(\sigma^2 = 1\) is known. With \(\sigma^2\) fixed, a Gaussian log-likelihood depends on the data only through the residual sum of squares (RSS): \(\ell= -\frac{n}{2}{\operatorname{log}\mathopen{}\left\{2\pi\sigma^2\right\}\mathclose{}} - \frac{\text{RSS}}{2\sigma^2}\). For a model \(M\) with residual sum of squares \(\text{RSS}_M\):
Saturated model: the fitted values are the pattern means \(2\), \(6\), and \(7\), so \(\text{RSS}_{\text{sat}} = (1-2)^2 + (3-2)^2 + (5-6)^2 + (7-6)^2 + (6-7)^2 + (8-7)^2 = 6\).
Simple linear regression model: the least-squares estimates use \({\bar{x}}= 1\) and \(\bar{y} = 30/6 = 5\):
The fitted values are \(2.5\), \(5\), and \(7.5\) at \(x = 0, 1, 2\), so \(\text{RSS}_M = (1-2.5)^2 + (3-2.5)^2 + (5-5)^2 + (7-5)^2 + (6-7.5)^2 + (8-7.5)^2 = 9\), and the residual deviance is \(D = (9 - 6)/1 = 3\).
Null model: the fitted value is \(\bar{y} = 5\) for every observation, so \(\text{RSS}_0 = (1-5)^2 + (3-5)^2 + (5-5)^2 + (7-5)^2 + (6-5)^2 + (8-5)^2 = 34\), and the null deviance is \(D_0 = (34 - 6)/1 = 28\).
The linear model’s residual deviance (\(3\)) is much smaller than the null deviance (\(28\)): the linear model captures most of the variation in pattern means that the null model misses.
Null Deviance
Similarly, the least complicated model we could fit would have only one mean parameter, an intercept:
And we can compare it to more complicated models using a likelihood ratio test:
[R code]
lrtest(bw_lm2, lm0)
The likelihood ratio statistic for the test comparing the null model to the maximal model is \[\lambda = 2 * (\ell_{\text{full}} - \ell_{0}) = 35.106732\] where:
\(\ell_{\text{0}}\) is the log-likelihood of the null model: -168.954967
\(\ell_{\text{full}}\) is the log-likelihood of the maximal model: -151.401601
In R, this test is:
[R code]
lrtest(lm_max, lm0)
This log-likelihood ratio statistic is called the null deviance (Definition 11). It tells us whether we have enough data to detect a difference between the null and full models.
Figure 11: Four models for birthweight data, with 95% confidence and prediction intervals. The spectrum from null to saturated includes many other possible models, such as quadratic polynomial models (with or without interactions) and generalized additive models (GAMs).
The residual deviance is a general concept that applies to all GLMs, including Gaussian linear regression. For any GLM, the residual deviance (Definition 9) is:
where \(\ell_{\text{sat}}\) is the log-likelihood of the saturated model (Definition 8) and \(\ell(\hat\beta)\) is the log-likelihood of the fitted model. However, this formula simplifies differently depending on the distribution family, and the resulting test statistics have different null distributions.
The saturated model has one free mean parameter per distinct covariate pattern. When there are no replicates (each covariate pattern appears exactly once), it sets \(\hat\mu_i = y_i\) for every observation, so its residual sum of squares is zero:
Because \(\sigma^2\) is unknown and must be estimated, comparing deviances between two nested Gaussian models uses an F-test, which is exact under Gaussian assumptions. Substituting \(D = \text{RSS}/\hat\sigma^2\), the \(\hat\sigma^2\) factors cancel:
For Poisson and Binomial models, the dispersion parameter \(\phi = 1\) is fixed rather than estimated. Consequently, the difference in deviances between two nested models follows an approximate \(\chi^2\) distribution:
\[D_1 - D_2 \;\dot \sim\; \chi^2_{p_2 - p_1}\]
This asymptotic result replaces the exact F-test used for Gaussian models.
Saturated vs. fully parametrized models with replicates
When some covariate patterns appear more than once (i.e., when there are replicates), it is important to distinguish between two special models:
The saturated model has one free mean parameter per distinct covariate pattern (\(q\) parameters, where \(q \leq n\) is the number of unique patterns).
The fully parametrized model has one free mean parameter per observation (\(n\) parameters).
When there are no replicates (\(q = n\)), these two coincide. When there are replicates (\(q < n\)), the saturated model constrains all observations sharing a covariate pattern to have the same mean, but places no constraint on means across different patterns. See Kleinbaum and Klein [8] for further discussion of this distinction.
Deviance is always measured relative to the saturated model, not the fully parametrized model.
Gaussian deviance with replicates
When covariate pattern \(k\) has \(n_k\) replicates with sample mean \(\bar{y}_k\), the saturated model fits \(\hat\mu_k = \bar{y}_k\) for each pattern \(k\), so its residual sum of squares equals the within-group (pure error) SS:
This within-group residual sum of squares is nonzero whenever any covariate pattern has replicates with different response values. The Gaussian deviance relative to the saturated model is therefore:
R’s deviance() applied to an lm object returns the total fitted-model RSS, not the deviance relative to the saturated model. To compute the deviance relative to the saturated model, subtract the saturated model’s RSS: deviance(lm_fit) - deviance(lm_saturated).
For the birthweight data, we can verify this directly using bw_lm2 (the interaction model) and lm_max (the saturated model):
deviance(bw_lm2) # total RSS of fitted model#> [1] 652425deviance(lm_max) # within-group (pure error) SS#> [1] 423783deviance(bw_lm2) -deviance(lm_max) # lack-of-fit SS (deviance vs. saturated)#> [1] 228642
GLM deviance with replicates
For a Binomial GLM fit to data grouped by covariate pattern (with \(y_k\) events and \(n_k\) observations for pattern \(k\)), the saturated model sets \(\hat\pi_k^{\text{sat}} = y_k / n_k\) for each pattern \(k\). R’s deviance() correctly computes \(2(\ell_{\text{sat}} - \ell(\hat\beta))\) using this grouping:
When binomial data is ungrouped (individual Bernoulli observations, \(n_i = 1\)), R uses the fully parametrized model as its reference — assigning \(\hat\pi_i = y_i \in \{0, 1\}\) to each individual observation. Since each predicted probability is exactly 0 or 1, every Bernoulli likelihood contribution equals 1, and hence \(\ell_{\text{fp}} = 0\). Thus R’s deviance() returns \(-2\ell(\hat\beta)\).
We can verify this directly: observations are ungrouped Bernoulli with repeated covariate patterns (\(x \in \{0, 1\}\) appears three times each).
This reference model is the fully parametrized model (one parameter per observation), not the saturated model (one parameter per distinct covariate pattern). They coincide only when there are no repeated covariate patterns (\(q = n\)). When patterns do repeat (\(q < n\)), the saturated model sets \(\hat\pi_k = y_k/n_k\) per pattern, giving \(\ell_{\text{sat}} < 0\).
deviance() for ungrouped data cannot be used as a goodness-of-fit test against the \(\chi^2\) distribution when \(q < n\). The correct GOF statistic is \(2(\ell_{\text{sat}} - \ell(\hat\beta))\), but R’s deviance() for ungrouped data returns \(-2\ell(\hat\beta)\) (using \(\ell_{\text{fp}} = 0\)). These two quantities differ by \(-2\ell_{\text{sat}} > 0\) whenever \(q < n\). Even if each covariate pattern has many replicates (large \(n_k\)), R’s ungrouped deviance is the wrong statistic to compare against \(\chi^2(q - p)\). To obtain a valid GOF test when patterns repeat, fit the model using grouped data (one row per pattern with \(y_k\) and \(n_k\)), so that R uses the saturated model as its reference.
The model assumes that the distribution conditional on a given \(X\) value is Gaussian.
Correct Functional Form of Conditional Mean Structure (Linear Component)
The model assumes that the conditional means have the structure: \[{\mathbb{E}\mathopen{}\left[Y|\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right]\mathclose{}} = {\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}}^{\top}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\]
Homoskedasticity
The model assumes that variance\(\sigma^2\) is constant (with respect to \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\)).
Figure 12: Birthweight model with interaction term
Fitted model for hers data
Consider the hers data from Vittinghoff et al. [3].
The “heart and estrogen/progestin study” (HERS) was a clinical trial of hormone therapy for prevention of recurrent heart attacks and death among 2,763 post-menopausal women with existing coronary heart disease (CHD) [20].
The trial was conducted at 20 US clinical centers. Participants were randomized to receive either conjugated equine estrogens (0.625 mg/day) plus medroxyprogesterone acetate (2.5 mg/day) or a matching placebo [20]. Women were followed for a sample mean of 4.1 years [20].
The primary outcome was nonfatal myocardial infarction or CHD death [20].
Definition 12 (Deviation of an observation from its subpopulation mean) The deviation of an observation from its subpopulation mean in a probabilistic model \(p(Y \mid X)\), is the difference between an observed value and its model-implied mean given covariates:
Definition 13 (Signal (statistical sense)) In statistical modeling, the signal is the deterministic part of the model. For mean-based models, the signal is the model-implied mean function, for example \({\mathbb{E}\mathopen{}\left[Y \mid X\right]\mathclose{}}\) (or \({\mathbb{E}\mathopen{}\left[Y\right]\mathclose{}}\) when there are no covariates).
Theorem 1 (Deviations in Gaussian models) If \(Y\) has a Gaussian distribution, then \(\mathbb{d}(Y)\) also has a Gaussian distribution, and vice versa.
Proof. By Definition 12, \(\mathbb{d}(Y \mid X = x) = Y - {\mathbb{E}\mathopen{}\left[Y \mid X = x\right]\mathclose{}}\). Let \(\mu(x) = {\mathbb{E}\mathopen{}\left[Y \mid X = x\right]\mathclose{}}\). Since \(\mathbb{d}(Y \mid X = x) = Y - \mu(x)\) is an affine function of \(Y\), and Gaussian distributions are closed under affine transformations, if \(Y \mid X = x \sim \operatorname{N}(\mu(x), \sigma^2)\), then:
\[
\mathbb{d}(Y \mid X = x) = Y - \mu(x) \sim \operatorname{N}(\mu(x) - \mu(x), \sigma^2) = \operatorname{N}(0, \sigma^2)
\]
Conversely, since \(Y = {\mathbb{E}\mathopen{}\left[Y \mid X = x\right]\mathclose{}} + \mathbb{d}(Y \mid X = x) = \mu(x) + \mathbb{d}(Y \mid X = x)\), if \(\mathbb{d}(Y \mid X = x) \sim \operatorname{N}(0, \sigma^2)\), then \(Y \mid X = x = \mu(x) + \mathbb{d}(Y \mid X = x) \sim \operatorname{N}(\mu(x), \sigma^2)\).
Definition 14 (Residuals) A residual is the difference between an observed value and its corresponding fitted value, \(\hat y\): \[
\begin{aligned}
\hat{\varepsilon}&\stackrel{\text{def}}{=}y - \hat y
\end{aligned}
\]
For indexed observations, the residual is defined as: \[
\hat{\varepsilon}_i \stackrel{\text{def}}{=}y_i - \hat y_i
\]
For a particular fitted model, each residual \(\hat{\varepsilon}_i\) is tied to its corresponding fitted value \(\hat y_i\). The fitted value \(\hat y_i\) is often a sample mean or fitted conditional mean, but not always.
Theorem 2 (Residuals are orthogonal to the covariates) Fit a linear regression model by maximum likelihood, with design matrix \(\mathbf{X}\) of full column rank, coefficient estimate \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\), fitted values \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}} \stackrel{\text{def}}{=}\mathbf{X}\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\), and residual vector \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\hat{\varepsilon}}}\vphantom{\hat{\varepsilon}}}\stackrel{\text{def}}{=}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}- \hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}}\). Then:
In particular, if the model has an intercept, the residuals add up to zero, \(\sum_{i=1}^n\hat{\varepsilon}_i = 0\), and the residuals are orthogonal to the fitted values, \(\sum_{i=1}^n\hat y_i \hat{\varepsilon}_i = 0\).
Proof. In Section 3, setting the score equal to zero showed that \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) solves \({\mathbf{X}}^{\top}\mathbf{X}\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}= {\mathbf{X}}^{\top}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}\). Then:
If the model has an intercept, the first column of \(\mathbf{X}\) is all ones, so the first entry of Equation 12 is \(\sum_{i=1}^n1 \cdot \hat{\varepsilon}_i = 0\). For the fitted values:
\[
\begin{aligned}
\sum_{i=1}^n\hat y_i \hat{\varepsilon}_i
&= {\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}}}^{\top} \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\hat{\varepsilon}}}\vphantom{\hat{\varepsilon}}}
&& \text{(sum of products as an inner product)}
\\
&= {\mathopen{}\left(\mathbf{X}\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}}^{\top} \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\hat{\varepsilon}}}\vphantom{\hat{\varepsilon}}}
&& \text{(definition of } \hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}} \text{)}
\\
&= {\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}}^{\top} {\mathbf{X}}^{\top} \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\hat{\varepsilon}}}\vphantom{\hat{\varepsilon}}}
&& \text{(transpose of a product)}
\\
&= {\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}}^{\top} \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{0}}\vphantom{0}}_{(p+1) \times 1}
&& \text{(the orthogonality just shown)}
\\
&= 0
&& \text{(multiply)}
\end{aligned}
\]
Example 9 (Residual sums in the birthweight model) For the birthweight model with an interaction term (bw_lm2), we can compute \({\mathbf{X}}^{\top}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\hat{\varepsilon}}}\vphantom{\hat{\varepsilon}}}\), \(\sum_{i=1}^n\hat{\varepsilon}_i\), and \(\sum_{i=1}^n\hat y_i \hat{\varepsilon}_i\) directly:
All three are zero up to floating-point rounding error, as Theorem 2 guarantees. Because the model includes an indicator for male babies, the entry of \({\mathbf{X}}^{\top}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\hat{\varepsilon}}}\vphantom{\hat{\varepsilon}}}\) for that indicator is zero too: the residuals for the male babies add up to zero, and so, subtracting from the overall sum, do the residuals for the female babies.
Example 10 (Without an intercept, the residuals need not add up to zero) If we drop the intercept and regress birthweight on gestational age alone, \(\mathbf{X}\) has a single column, the ages, so Theorem 2 only guarantees \(\sum_{i=1}^n\text{age}_i \hat{\varepsilon}_i = 0\):
Different sources are not fully consistent about these terms. For terminology in this course, we use:
deviation for deviation of an observation from its population mean (typically conditional on covariates). This quantity is often called error in other sources, even though it doesn’t directly involve an estimate; cf., https://en.wikipedia.org/wiki/Errors_and_residuals;
Theorem 3 (Mean and variance of residuals) For an ordinary least squares linear model with fitted values \(\hat y_i = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) (and fitted-value vector \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}}\)), if the conditional mean is correctly specified so that:
and if the errors are homoskedastic so that \({\operatorname{Var}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\mid \mathbf{X}\right)\mathclose{}} = \sigma^2\mathbb{I}_n\) where \(\mathbb{I}_n\) is the \(n \times n\)identity matrix, then, treating \(\mathbf{X}\) as fixed in this derivation, the residual moments are:
The conditional mean result (Equation 13) uses unbiasedness of fitted values. The conditional variance result (Equation 14) uses the homoskedasticity condition with the OLS hat-matrix representation.
When leverage is roughly uniform, \(h_{ii} \approx k/n\), where \(k = \mathrm{rank}(\mathbf{X})\) (so \(k = p + 1\) for a full-rank model with an intercept and \(p\)predictors), so \({\operatorname{Var}\mathopen{}\left(\hat{\varepsilon}_i\right)\mathclose{}} \approx \sigma^2\) as \(n\) grows.
Theorem 4 (Conditional distribution of residuals in Gaussian models) In a Gaussian linear regression model, the residuals\(\hat{\varepsilon}_i\), conditional on the design matrix \(\mathbf{X}\), are normally distributed:
where \(h_{ii}\) is the \(i\)-th diagonal element of the hat matrix (Definition 15). When leverage is roughly uniform, \(\sigma^2(1-h_{ii}) \approx \sigma^2\) for large \(n\); replacing the unknown \(\sigma^2\) with its estimate \(\hat{\sigma}^2\) gives:
Proof. Using the hat matrix (Definition 15), the residual vector is \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\hat{\varepsilon}}}\vphantom{\hat{\varepsilon}}} = (\mathbb{I}_n-H)\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\), so the \(i\)-th residual is a linear combination of the observations:
where \(\delta_{ij}\) is the Kronecker delta (\(\delta_{ij} = 1\) if \(i = j\), else \(\delta_{ij} = 0\)).
Given \(\mathbf{X}\), the random variables\(Y_1, \ldots, Y_n\) are conditionally independent Gaussians by the Gaussian model assumption. Therefore, conditional on \(\mathbf{X}\), the residual \(\hat{\varepsilon}_i\) is a linear combination of independent Gaussians, so \(\hat{\varepsilon}_i\) is Gaussian.
In particular, for \(i \neq j\), these covariances are generally nonzero.
Marginal distributions of residuals
To look for problems with our model, we can check whether the residuals\(\hat{\varepsilon}_i\) and standardized residuals\({\hat{\varepsilon}'{}}_i\) look like they have the distributions that they are supposed to have, according to the model.
Figure 21: Marginal distribution of standardized residuals
QQ plot of standardized residuals
[R code]
library(ggfortify)# needed to make ggplot2::autoplot() work for `lm` objectsqqplot_lm2_auto <- bw_lm2 |>autoplot(which =2, # options are 1:6; can do multiple at oncencol =1 ) +theme_classic()print(qqplot_lm2_auto)
Figure 22: Typical QQ-plot patterns for six types of residual distribution. Adapted from Dunn and Smyth [2, Fig. 3.6], with thanks to the authors; reproduced here using simulated data (\(n = 150\)). Each column shows one distribution scenario: a QQ plot (top) paired with a histogram of the simulated data (bottom).
Figure 23: Three equivalent ways to produce a QQ plot of the standardized residuals for the birthweight model (Equation 2). All three plots show the same data and reference line.
Levene / Brown–Forsythe test (a median-centered Levene variant, where standard Levene centers on group sample means, and Brown–Forsythe centers on group medians for more robustness; e.g., via car::leveneTest(..., center = median) or equivalent code).
Levene / Brown–Forsythe test (homoskedasticity across groups)
Levene’s test transforms residuals to within-group absolute deviations: \[
z_{gi} = |e_{gi} - c_g|,
\] where \(c_g\) is the group center.
Classical Levene uses the group sample mean for \(c_g\). Brown–Forsythe uses the group median, which is more robust.
Then run a one-way ANOVA on \(z_{gi}\) by group: \[
F = \frac{\text{MS}_{\text{between}}}{\text{MS}_{\text{within}}}
\sim F_{G-1, N-G}
\quad\text{under }H_0.
\]
Small p-values suggest unequal residual variance.
For simple linear regression, Kutner et al. [4, pp. 116–117] describes the Brown–Forsythe test by splitting observations into two \(X\)-level groups (low versus high), computing absolute deviations from each group median, and applying a two-sample pooled-variance t test: let \[
z_{ij} = |e_{ij} - \tilde e_i|,
\] where \(j\) indexes observations within group \(i\), and \(\tilde e_i\) is the median residual in group \(i\). Then: \[
t_{\text{BF}} =
\frac{\bar z_{1} - \bar z_{2}}
{s_p \sqrt{1/n_{1} + 1/n_{2}}},
\quad
t_{\text{BF}} \approx t_{n_{1}+n_{2}-2}
\text{ under }H_0.
\] Here, \(\bar z_{1}\) and \(\bar z_{2}\) are the sample means of the \(z_{ij}\) values in groups \(i=1\) and \(i=2\), \(s_p\) is their pooled standard deviation, and \(n_{1}, n_{2}\) are the two group sample sizes. Large \(|t_{\text{BF}}|\) indicates nonconstant residual variance.
Breusch–Pagan test (variance changing with covariates)
The Fligner–Killeen and Levene tests compare residual variances across the levels of a grouping variable. The Breusch–Pagan test instead asks whether the error variance changes with the covariates themselves, which can be continuous [21].
Definition 17 (Auxiliary regression) Let \(\hat{\varepsilon}_1, \ldots, \hat{\varepsilon}_n\) be the residuals from a fitted linear regression model, and let \(z_i\) be a vector of \(q\) covariates for observation \(i\) (by default, the model’s own predictors, without the intercept).
The auxiliary regression is the ordinary least-squares regression of the squared residuals \(\hat{\varepsilon}_i^2\) on \(z_i\), with an intercept; let \(R^2_{\text{aux}}\) be its coefficient of determination.
Definition 18 (Studentized Breusch–Pagan statistic) With \(R^2_{\text{aux}}\) the coefficient of determination of the auxiliary regression (Definition 17), the studentized Breusch–Pagan statistic[22] is:
Definition 19 (Original Breusch–Pagan statistic) In contrast to Equation 17, the original Breusch–Pagan statistic[21] uses \(\hat{\sigma}^2 \stackrel{\text{def}}{=}\frac{1}{n}\sum_{i=1}^n\hat{\varepsilon}_i^2\) and regresses \(g_i \stackrel{\text{def}}{=}\hat{\varepsilon}_i^2 / \hat{\sigma}^2\) on \(z_i\); the statistic is half of that regression’s explained sum of squares.
Theorem 5 (Null distribution of the Breusch–Pagan statistic) Under the null hypothesis that the error variance is constant, \({\operatorname{Var}\mathopen{}\left(Y_i \mid z_i\right)\mathclose{}} = \sigma^2\) for every \(i\), both versions of the Breusch–Pagan statistic are approximately \(\chi^2_q\) in large samples. The original version also assumes Gaussian errors [21]; the studentized version does not [22].
Large values of the statistic, and small p-values, suggest that the error variance changes with \(z_i\).
Example 11 (Breusch–Pagan test for the birthweight interaction model) For the birthweight interaction model (Equation 2), \(z_i\) holds sex, gestational age, and their interaction, so \(q = 3\). We fit both auxiliary regressions by hand and compare the results with lmtest::bptest():
\[
\begin{aligned}
\text{BP}
&\stackrel{\text{def}}{=}n R^2_{\text{aux}}
&& \text{(definition of the studentized statistic)}
\\
&= 24 \times 0.10435
&& \text{(substitute } n \text{ and } R^2_{\text{aux}}\text{)}
\\
&\approx 2.504
&& \text{(arithmetic)}
\end{aligned}
\tag{18}\]
and the original statistic is:
\[
\begin{aligned}
\text{BP}_{\text{orig}}
&\stackrel{\text{def}}{=}\tfrac{1}{2} \text{ESS}_g
&& \text{(definition of the original statistic)}
\\
&= \tfrac{1}{2} \times 1.7742
&& \text{(substitute the explained sum of squares)}
\\
&\approx 0.887
&& \text{(arithmetic)}
\end{aligned}
\tag{19}\]
Comparing each value with a \(\chi^2_3\) distribution (Theorem 5) gives \(p \approx 0.47\) for the studentized statistic and \(p \approx 0.83\) for the original statistic, which match bptest(). Neither test gives evidence that the error variance changes with sex or gestational age. The two statistics differ noticeably. The original statistic divides the explained sum of squares by 2, which is the variance of \(\varepsilon_i^2 / \sigma^2\) when the errors are Gaussian, and so the variance that \(g_i\) approaches in large samples under Gaussian errors; the studentized statistic equals the same explained sum of squares divided by the mean squared deviation of the \(g_i\) values (their variance with divisor \(n\)) instead, which here is about 0.71. A non-significant result from only 24 babies is weak evidence that the variance is constant, so we would still look at the residual plots.
Shapiro–Wilk test (normality of standardized residuals)
For ordered standardized residuals\(r_{(1)} \le \cdots \le r_{(n)}\), the Shapiro–Wilk statistic is: \[
W =
\frac{\mathopen{}\left(\sum_{i=1}^n a_i r_{(i)}\right)\mathclose{}^2}
{\sum_{i=1}^n (r_i - \bar r)^2},
\] where \(a_i\) are constants from normal-order-statistic moments. The numerator uses ordered residuals \(r_{(i)}\), while the denominator uses the original (unordered) residuals.
If residuals are Gaussian, \(W\) tends to be close to 1. Small \(W\) (and small p-value) indicates departure from normality.
Interpretation rule: for all three tests in this table, a small p-value is evidence against the corresponding model assumption.
Compared with visual diagnostics:
Fligner–Killeen / Levene summarizes the same heteroskedasticity signal that we inspect in residuals-vs-fitted (Figure 24) and scale-location (Figure 32) plots.
This mean and variance structure should be true for every value of \(x\).
If we didn’t correctly guess the functional form of the linear component of the mean, for covariate vector \(x = (x_1,\ldots,x_p)\), \[{\mathbb{E}\mathopen{}\left[Y \mid X=x\right]\mathclose{}} = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p\]
Then the residuals might have nonzero mean.
Regardless of whether we guessed the mean function correctly, the variance of the residuals might differ between values of \(x\).
Residuals versus fitted values
[R code]
autoplot(bw_lm2, which =1, ncol =1) |>print()
Figure 24: birthweight model (Equation 2): residuals versus fitted values
Figure 31: Number of authors versus title length in PLOS Medicine
Scale-location plot
Definition 20 (Scale-location plot) A scale-location plot graphs \(\sqrt{\mathopen{}\left|{\hat{\varepsilon}'{}}_i\right|\mathclose{}}\), the square root of the absolute value of each standardized residual, against the corresponding fitted value\(\hat y_i\).
Example 12 (Scale-location plot for the birthweight model)Figure 32 shows the scale-location plot for the birthweight model with interactions (Equation 2).
[R code]
autoplot(bw_lm2, which =3, ncol =1) |>print()
Figure 32: Scale-location plot of birthweight data
The blue trend line is roughly flat, so this plot shows no clear sign that the residual variance changes with the fitted value.
Unlike the residuals-versus-fitted plot (Figure 24), the blue line in a scale-location plot doesn’t need to be near 0, but it should be flat. If the standardized residuals have a standard Gaussian distribution, the line should sit near \({\mathbb{E}\mathopen{}\left[\sqrt{\mathopen{}\left|Z\right|\mathclose{}}\right]\mathclose{}} \approx 0.82\), where \(Z \sim {\operatorname{N}\mathopen{}\left(0,1\right)\mathclose{}}\):
gamma(3/4) *2^(1/4) /sqrt(pi)#> [1] 0.822179
If the line trends upward or downward, the residual variance \(\sigma^2\) might not be constant. In a residuals-versus-fitted plot, the same problem shows up as a funnel shape; transforming the outcome \(Y\) (for example, with \(\log{Y}\) or \(\sqrt{Y}\)) or using a model that allows non-constant variance can help [23, Sec. 3.3.3].
Example 13 (Scale-location plot when the variance is not constant)Figure 33 shows a scale-location plot for simulated data whose residual standard deviation grows in proportion to the mean. The blue line slopes clearly upward, which is the pattern that signals non-constant variance.
[R code]
set.seed(1)sim_hetero <- tibble::tibble(x =runif(200, min =1, max =10),y =rnorm(200, mean =2* x, sd =0.5* x))lm(y ~ x, data = sim_hetero) |>autoplot(which =3, ncol =1) |>print()
Figure 33: Scale-location plot for simulated data with non-constant variance
Residuals versus leverage
[R code]
autoplot(bw_lm2, which =5, ncol =1) |>print()
Figure 34: birthweight model with interactions (Equation 2): residuals versus leverage
Cook’s distance
Tip
This section is adapted from Dobson and Barnett [1, Sec. 6.3].
Residuals and leverage each capture one aspect of an observation’s influence on the model fit. Cook’s distance[24] combines them into a single diagnostic.
Definition 21 (Cook’s distance) For ordinary least squares with \(p\) parameters, the Cook’s distance for observation \(i\) is
where \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}_{(-i)}\) is the least-squares estimate obtained after deleting observation \(i\).
Because refitting the model \(n\) times is expensive, Equation 20 has an equivalent closed-form expression that uses only quantities from the full-data fit:
resid_vs_fit <- bw |>ggplot(aes(x = predlm2, y = residlm2, col = sex, shape = sex) ) +geom_point() +theme_classic() +geom_hline(yintercept =0)print(resid_vs_fit)
Standardized residuals vs fitted
[R code]
bw |>ggplot(aes(x = predlm2, y = std_resid, col = sex, shape = sex) ) +geom_point() +theme_classic() +geom_hline(yintercept =0)
Standardized residuals vs gestational age
[R code]
bw |>ggplot(aes(x = age, y = std_resid, col = sex, shape = sex) ) +geom_point() +theme_classic() +geom_hline(yintercept =0)
sqrt(abs(std_resid)) vs fitted
Compare with autoplot(bw_lm2, 3) (note: autoplot uses leverage-adjusted rstandard())
[R code]
bw |>ggplot(aes(x = predlm2, y = sqrt_abs_std_resid, col = sex, shape = sex) ) +geom_point() +theme_classic() +geom_hline(yintercept =0)
Diagnostics for the independence assumption
The independence assumption means that the model errors are independent, or at least uncorrelated, across observations after conditioning on the predictors in the model. Because those errors are unobserved, we usually assess this assumption using residual-based plots and formal tests.
For data with a natural ordering (for example, time, spatial location, or clinic visit sequence), we usually assess independence with both plots and formal tests [4, Ch. 12], [5, Ch. 6].
Definition 22 (Common diagnostics for independence)
Residuals versus observation order (or time) plots: patterns, runs, or drifts can indicate dependence [4, Ch. 12].
Correlograms: plot the sample autocorrelation function (ACF), and often the partial autocorrelation function (PACF), to look for serial structure [5, Ch. 6].
Mathematical details for Durbin-Watson and Breusch-Godfrey
Suppose the observations have a meaningful order, indexed by \(t = 1, \ldots, n\). Let \(e_t\) denote the OLS residual from a fitted regression model.
Definition 23 (Durbin-Watson diagnostic) The test statistic is \[
\begin{aligned}
d
&= \frac{\sum_{t = 2}^{n} (e_t - e_{t-1})^2}
{\sum_{t = 1}^{n} e_t^2}
\end{aligned}
\] which is approximately \(2(1 - \hat r_1)\), where \(\hat r_1\) is the sample lag-1 residual autocorrelation [4, Ch. 12], [5, Ch. 6]. This approximation is most useful in standard large-sample linear-model settings, and is less straightforward in models with lagged outcomes [4, Ch. 12].
The null and alternatives are typically framed through an AR(1) error model: \[
\varepsilon_t = \rho \varepsilon_{t-1} + u_t.
\] Here, \(\varepsilon_t\) is the autocorrelated regression error process, and \(u_t\) is a white-noise innovation term. The null is \(H_0: \rho = 0\), and alternatives can be one-sided or two-sided, depending on whether we suspect positive, negative, or any serial correlation[4, Ch. 12].
Definition 24 (Breusch-Godfrey diagnostic) For a test of order \(p\), we run an auxiliary regression of residuals on: the original regressors, and lagged residuals \(e_{t-1}, \ldots, e_{t-p}\). If \(R^2_{\text{aux}}\) is from that auxiliary model, the LM statistic is \[
\text{LM} = (n - p)R^2_{\text{aux}},
\] and under \(H_0: \rho_1 = \cdots = \rho_p = 0\), it is asymptotically \(\chi^2_p\)[5, Ch. 6], [26, Ch. 11].
Example 15 (Simulated example for independence diagnostics) The following simulation generates ordered data with positively autocorrelated errors. This serial correlation creates a setting where independence should fail. We generate the error term from an AR(1) process with coefficient 0.7, matching the notation in the mathematical details section.
set.seed(204)n_obs <-120x <-seq(-1, 1, length.out = n_obs)ar_errors <-as.numeric(arima.sim(model =list(ar =0.7), n = n_obs, sd =1))y <-2+1.5* x + ar_errorsserial_dep_exm <- tibble::tibble(time_index =seq_len(n_obs),x = x,y = y)serial_dep_exm_lm <-lm(y ~ x, data = serial_dep_exm)summary(serial_dep_exm_lm)#> #> Call:#> lm(formula = y ~ x, data = serial_dep_exm)#> #> Residuals:#> Min 1Q Median 3Q Max #> -2.6747 -0.9128 -0.0552 1.0087 2.5022 #> #> Coefficients:#> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 1.998 0.110 18.13 < 2e-16 ***#> x 1.773 0.189 9.37 6.4e-16 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#> #> Residual standard error: 1.21 on 118 degrees of freedom#> Multiple R-squared: 0.426, Adjusted R-squared: 0.422 #> F-statistic: 87.7 on 1 and 118 DF, p-value: 6.36e-16
lmtest::dwtest(serial_dep_exm_lm, alternative ="two.sided")#> #> Durbin-Watson test#> #> data: serial_dep_exm_lm#> DW = 0.7623, p-value = 4.16e-12#> alternative hypothesis: true autocorrelation is not 0lmtest::bgtest(serial_dep_exm_lm, order =4)#> #> Breusch-Godfrey test for serial correlation of order up to 4#> #> data: serial_dep_exm_lm#> LM test = 45.94, df = 4, p-value = 2.53e-09
In this example, the residual-order plot and correlogram suggest serial dependence, and both formal tests reject independence.
Example 16 (Real-data examples for independence diagnostics) We can also run the same diagnostics on real ordered data sets.
Nile annual river flow
nile_df <- tibble::tibble(year =as.numeric(time(datasets::Nile)),flow =as.numeric(datasets::Nile))nile_lm <-lm(flow ~ year, data = nile_df)summary(nile_lm)#> #> Call:#> lm(formula = flow ~ year, data = nile_df)#> #> Residuals:#> Min 1Q Median 3Q Max #> -483.7 -98.2 -23.2 111.4 368.7 #> #> Coefficients:#> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 6132.174 1001.758 6.12 1.9e-08 ***#> year -2.714 0.522 -5.20 1.1e-06 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#> #> Residual standard error: 151 on 98 degrees of freedom#> Multiple R-squared: 0.217, Adjusted R-squared: 0.209 #> F-statistic: 27.1 on 1 and 98 DF, p-value: 1.07e-06lmtest::dwtest(nile_lm, alternative ="two.sided")#> #> Durbin-Watson test#> #> data: nile_lm#> DW = 1.247, p-value = 9.37e-05#> alternative hypothesis: true autocorrelation is not 0lmtest::bgtest(nile_lm, order =4)#> #> Breusch-Godfrey test for serial correlation of order up to 4#> #> data: nile_lm#> LM test = 15.94, df = 4, p-value = 0.0031
Figure 38: Residuals versus year for Nile flow model
Mauna Loa atmospheric carbon dioxide (co2)
co2_df <- tibble::tibble(decimal_year =as.numeric(time(datasets::co2)),co2_ppm =as.numeric(datasets::co2))co2_lm <-lm(co2_ppm ~ decimal_year, data = co2_df)summary(co2_lm)#> #> Call:#> lm(formula = co2_ppm ~ decimal_year, data = co2_df)#> #> Residuals:#> Min 1Q Median 3Q Max #> -6.040 -1.948 -0.002 1.911 6.515 #> #> Coefficients:#> Estimate Std. Error t value Pr(>|t|) #> (Intercept) -2.25e+03 2.13e+01 -106 <2e-16 ***#> decimal_year 1.31e+00 1.07e-02 122 <2e-16 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#> #> Residual standard error: 2.62 on 466 degrees of freedom#> Multiple R-squared: 0.969, Adjusted R-squared: 0.969 #> F-statistic: 1.48e+04 on 1 and 466 DF, p-value: <2e-16lmtest::dwtest(co2_lm, alternative ="two.sided")#> #> Durbin-Watson test#> #> data: co2_lm#> DW = 0.2124, p-value <2e-16#> alternative hypothesis: true autocorrelation is not 0lmtest::bgtest(co2_lm, order =12)#> #> Breusch-Godfrey test for serial correlation of order up to 12#> #> data: co2_lm#> LM test = 453.9, df = 12, p-value <2e-16
For both real-data examples, the diagnostics indicate residual dependence, consistent with their time-ordered structure.
No single diagnostic is definitive. In practice, we combine visual and formal diagnostics, and interpret them in the context of study design [4, Ch. 12], [26, Ch. 11].
4.9 Model selection
(adapted from Dobson and Barnett [1] §6.3.3; for more information on prediction, see James et al.[27] and Harrell [28]).
DAGs for variable selection
For explanatory models, variable inclusion should not rely on automated algorithms alone. Directed acyclic graphs (DAGs) can help us encode substantive assumptions about which variables are confounders, mediators, or colliders.
In a Dobson-style workflow, we use the DAG first to decide a defensible adjustment set, then compare candidate regression models within that set using predictive or likelihood-based criteria. This approach keeps model selection aligned with study design, rather than only with numerical fit.
In this workflow, the DAG encodes hypothesized time ordering and causal pathways. This structure helps us decide which variables belong in the candidate model set before we run stepwise, subset, or penalized selection methods. The key point is that a variable can improve apparent fit while still distorting the effect we care about if it sits on a causal pathway or induces collider bias [1, Sec. 6.3.3].
Example 17 (high temperatures and mortality) Suppose the research question is: “What is the effect of high temperatures on deaths?”
A plausible DAG includes: high temperatures \(\rightarrow\) deaths, high temperatures \(\rightarrow\) blackouts, high temperatures \(\rightarrow\) forest fires, blackouts \(\rightarrow\) deaths, forest fires \(\rightarrow\) air pollution, and both forest fires and air pollution \(\rightarrow\) deaths [1, Sec. 6.3.3].
[R code]
dag_heat_deaths <- ggdag::dagify( D ~ T + B + F + P, B ~ T, F ~ T, P ~ F,exposure ="T",outcome ="D",labels =c(T ="High temperatures",D ="Deaths",B ="Blackouts",F ="Forest fires",P ="Air pollution" ))dag_heat_deaths |> ggdag::ggdag(use_labels ="label") + ggdag::theme_dag()
Figure 39: DAG for the high temperatures and mortality example.
If our estimand is the total effect of high temperatures on deaths, then blackouts, forest fires, and air pollution are mediators on downstream pathways. Adjusting for them would remove part of the effect we are trying to estimate.
But if the estimand changes to the effect of blackouts on deaths, then high temperatures become a confounder because they are a common cause of both blackouts and deaths. So the same variable can be “adjust for” or “do not adjust for” depending on the causal question.
Example 18 (salmonellosis and weather) For daily salmonellosis counts, a plausible DAG has temperature and rainfall as direct causes of risk, with both variables also affecting humidity. If humidity has no direct arrow to salmonellosis, then humidity is not a causal driver; it mainly reflects shared variation from temperature and rainfall. Including humidity in the regression is therefore usually unnecessary for estimating the weather effects, and it may reduce interpretability because humidity is a collider of temperature and rainfall. Conditioning on humidity can induce misleading associations and bias weather-effect estimates when that conditioning opens a path between temperature and rainfall that is not otherwise blocked [1, Sec. 6.3.3].
[R code]
dag_salmonellosis_weather <- ggdag::dagify( S ~ T + R + C, H ~ T + R, T ~ Se, R ~ Se,labels =c(S ="Salmonellosis",T ="Temperature",R ="Rainfall",H ="Humidity",C ="Test change",Se ="Season" ))dag_salmonellosis_weather |> ggdag::ggdag(use_labels ="label") + ggdag::theme_dag()
Figure 40: DAG for the salmonellosis and weather example.
The same DAG can include a “change in test” indicator when surveillance switched to a more sensitive diagnostic test. That variable can improve predictive accuracy and face validity, but omitting it does not invalidate causal interpretation for temperature and rainfall effects if it is not on their causal paths [1, Sec. 6.3.3].
Likewise, season may cause both temperature and rainfall, but if season has no direct path to salmonellosis beyond those weather variables, including season can obscure interpretation of the proximal weather effects [1, Sec. 6.3.3].
Finally, DAG-informed variable definitions matter. In one intensive-care example, a “nutrition score” looked highly predictive of mortality until investigators clarified that people near death were often coded as zero because no assessment was completed. So the recorded variable combined nutrition with prognosis, and naive selection would overstate its scientific meaning [1, Sec. 6.3.3].
Mean squared error
We might want to minimize the population mean squared error (pMSE), \({\mathbb{E}\mathopen{}\left[(y-\hat y)^2\right]\mathclose{}}\), for new observations that weren’t in our data set when we fit the model.
Unfortunately, \[\frac{1}{n}\sum_{i=1}^n (y_i-\hat y_i)^2\] gives a biased estimate of \({\mathbb{E}\mathopen{}\left[(y-\hat y)^2\right]\mathclose{}}\) for new data. If we want an unbiased estimate, we will have to be clever.
This optimism in sample fit is one reason that \(R^2\) is not enough for model selection. \(R^2\) does not decrease as we add explanatory variables, even when those variables do not improve out-of-sample prediction. Selecting a model by maximizing \(R^2\) alone can lead to overfitting, since \(R^2\) does not penalize this non-decreasing property.
In many applications, we split the data into training, validation, and test sets. This validation approach extends the basic train/test split by separating model tuning from final model assessment. This framework follows [23, pp. 198–201].
Use the training set to estimate model parameters.
Use the validation set to compare candidate models, choose transformations, or tune hyperparameters.
Use the test set once at the end to estimate final out-of-sample performance.
Keeping the test set untouched during model building helps avoid optimistic bias from repeatedly trying many models. If data are limited, \(k\)-fold cross-validation can replace a single validation split for more stable tuning.
Example 19 (Numerical example) The following example (Example 19) uses R’s built-in mtcars dataset (\(n=32\) cars) to predict fuel efficiency (mpg) from vehicle weight (wt). It uses one random split to compare linear, quadratic, cubic, and quartic models (mpg ~ wt, mpg ~ wt + I(wt^2), mpg ~ wt + I(wt^2) + I(wt^3), and mpg ~ wt + I(wt^2) + I(wt^3) + I(wt^4)) on a validation set, then reports the chosen model’s test sRMSE on untouched test data [23, p. 213].
Figure 44: Quartic model fit superimposed on data partitions
[R code]
validation_results
Table 25: Training and validation sRMSE for candidate models
The chosen model is the linear model, because it has the lower validation sRMSE.
[R code]
performance_comparison
Table 26: Validation and test sRMSE for the chosen model
4.11 Cross-validation
Rather than one arbitrary split, \(k\)-fold cross-validation repeatedly partitions the data into training and test folds. For each split we compute prediction error, then summarize errors across folds and replications. The preferred model has lower prediction error, with simpler models favored when errors are similar.
When the number of candidate explanatory variables is small, we can compare all possible subsets (\(2^p\) models) using the same cross-validation metric.
[R code]
data("carbohydrate", package ="dobson")library(cvTools)full_model <-lm(carbohydrate ~ ., data = carbohydrate)cv_full <- full_model |>cvFit(data = carbohydrate, K =5, R =10,y = carbohydrate$carbohydrate )reduced_model <- full_model |>update(formula =~ . - age)cv_reduced <- reduced_model |>cvFit(data = carbohydrate, K =5, R =10,y = carbohydrate$carbohydrate )
When the number of candidate predictors is modest, we can use best subset selection. For each model size \(k = 0, 1, \ldots, p\), we fit all \(\binom{p}{k}\) models and keep the best model of size \(k\) (e.g., by lowest RSS in the training data).
Best subset selection yields at most \(p+1\) candidate models to compare across model sizes, using criteria such as cross-validated prediction error, \(C_p\), AIC/BIC, or adjusted \(R^2\). In that sense, best subset selection is more exhaustive than one-path methods like forward or backward stepwise selection.
Stepwise methods are another common approach. Forward selection adds variables one at a time. Backward selection starts with the full model and removes variables sequentially. Both approaches can select different models from the same data.
olsrr::ols_step_both_aic(full_model)#> #> #> Stepwise Summary #> -------------------------------------------------------------------------#> Step Variable AIC SBC SBIC R2 Adj. R2 #> -------------------------------------------------------------------------#> 0 Base Model 140.773 142.764 83.068 0.00000 0.00000 #> 1 protein (+) 137.950 140.937 80.438 0.21427 0.17061 #> 2 weight (+) 132.981 136.964 77.191 0.44544 0.38020 #> -------------------------------------------------------------------------#> #> Final Model Output #> ------------------#> #> Model Summary #> ---------------------------------------------------------------#> R 0.667 RMSE 5.505 #> R-Squared 0.445 MSE 30.301 #> Adj. R-Squared 0.380 Coef. Var 15.879 #> Pred R-Squared 0.236 AIC 132.981 #> MAE 4.593 SBC 136.964 #> ---------------------------------------------------------------#> RMSE: Root Mean Square Error #> MSE: Mean Square Error #> MAE: Mean Absolute Error #> AIC: Akaike Information Criteria #> SBC: Schwarz Bayesian Criteria #> #> ANOVA #> -------------------------------------------------------------------#> Sum of #> Squares DF Mean Square F Sig. #> -------------------------------------------------------------------#> Regression 486.778 2 243.389 6.827 0.0067 #> Residual 606.022 17 35.648 #> Total 1092.800 19 #> -------------------------------------------------------------------#> #> Parameter Estimates #> ----------------------------------------------------------------------------------------#> model Beta Std. Error Std. Beta t Sig lower upper #> ----------------------------------------------------------------------------------------#> (Intercept) 33.130 12.572 2.635 0.017 6.607 59.654 #> protein 1.824 0.623 0.534 2.927 0.009 0.509 3.139 #> weight -0.222 0.083 -0.486 -2.662 0.016 -0.397 -0.046 #> ----------------------------------------------------------------------------------------
4.14 Lasso
Lasso is a penalized regression method that shrinks coefficient estimates toward zero. It adds an \(L_1\) penalty term to the objective function, creating a trade-off between model fit and parsimony. The intercept is not penalized.
As the tuning parameter \(\lambda\) increases, more coefficients are shrunken strongly, and some become exactly zero. So lasso both regularizes the model and performs variable selection. In practice, \(\lambda\) is usually chosen by cross-validation to minimize prediction error.
coef(cvfit, s ="lambda.1se")#> 4 x 1 sparse Matrix of class "dgCMatrix"#> lambda.1se#> (Intercept) 34.5502388#> age . #> weight -0.0510494#> protein 0.5472281
The reduced model (Equation 22) has \(p_1\) parameters. The full model (Equation 23) adds \(q\) extra predictors and has \(p_2 = p_1 + q\) parameters. The reduced model is a special case of the full model with the constraint \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\gamma}}\vphantom{\gamma}} = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{0}}\vphantom{0}}\).
\(H_0: \beta_{AM} = 0\) vs. \(H_A: \beta_{AM} \neq 0\).
F-statistic
Definition 26 (Partial F-statistic) Let \(\text{RSS}_0\) and \(\text{RSS}_1\) denote the residual sums of squares from the reduced and full models, respectively. The partial F-statistic is:
\(\text{RSS}_0 = \sum_{i=1}^n(y_i - \hat{y}_i^{(0)})^2\) is the residual SS under \(H_0\)
\(\text{RSS}_1 = \sum_{i=1}^n(y_i - \hat{y}_i^{(1)})^2\) is the residual SS under \(H_A\)
\(q = p_2 - p_1\) is the number of constraints (extra parameters in the full model)
\(n - p_2\) is the residual degrees of freedom of the full model
Theorem 6 (Null distribution of the partial F-statistic) Under \(H_0\) and the Gaussian linear regression assumptions,
\[
F \ \sim \ F_{q,\; n - p_2}
\tag{25}\]
This exact distribution is non-asymptotic: it holds for any sample size \(n\) when the errors \(\epsilon_i \ \sim_{\operatorname{iid}}\ N(0, \sigma^2)\).
Proof. Under \(H_0\), the extra predictors \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}}\) contribute nothing, so the numerator \(\text{RSS}_0 - \text{RSS}_1 \ \sim \ \sigma^2 \chi^2_q\) and the denominator \(\text{RSS}_1 \ \sim \ \sigma^2 \chi^2_{n-p_2}\) are independent chi-squared random variables (this independence property follows from the Gauss-Markov theorem and the properties of projections in the column spaces of the design matrices). Therefore,
From Section 4.1.8, the Gaussian deviance of model \(k\) is \(D_k = \text{RSS}_k / \hat\sigma^2\), where \(\hat\sigma^2\) is an estimate of \(\sigma^2\). Because \(\hat\sigma^2\) appears in both the numerator and denominator of Equation 24, it cancels:
For Gaussian linear regression, the MLE of \(\sigma^2\) under model \(k\) is \(\hat\sigma^2_k = \text{RSS}_k / n\), so the log-likelihood at the MLE is:
So \(\lambda \approx q \cdot F\) for large \(n\), and both statistics have the same asymptotic null distribution \(\chi^2_q\) (since \(q \cdot F_{q, n-p_2} \dot \sim \chi^2_q\) as \(n \to \infty\)).
Why use the F-test instead of the LRT?
For Gaussian linear regression:
The F-test is exact — it has the correct \(F_{q,n-p_2}\) null distribution for any sample size \(n\), provided the errors are Gaussian. It accounts for the estimation of \(\sigma^2\) via the residual degrees of freedom.
The LRT is approximate — it relies on the asymptotic \(\chi^2_q\) distribution, which requires large \(n\) and treats \(\sigma^2\) as known (fixed at its MLE).
Both tests give similar p-values for large \(n\). For small to moderate \(n\), the F-test is preferred because it is exact.
For non-Gaussian GLMs (Poisson, Binomial), the F-test is not applicable; the LRT (or Wald test) is the standard approach.
For Gaussian linear regression, the covariance formula holds exactly rather than approximately, and \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) is unbiased:
Theorem 7 (Exact mean and covariance of \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\)) Under the Gaussian linear regression model (Definition 1), let \(\mathbf{X}\) be the \(n \times (p+1)\) design matrix whose rows are \({\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i}^{\top}\), and suppose \(\mathbf{X}\) has full column rank. Then, conditional on \(\mathbf{X}\), the maximum likelihood estimator \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}= {\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}{\mathbf{X}}^{\top}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\) is unbiased for \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) and has covariance matrix \(\sigma^2{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}\):
These results are exact for every sample size, not large-sample approximations.
Proof. Write \(\mathbf{A} \stackrel{\text{def}}{=}{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}{\mathbf{X}}^{\top}\), a \((p+1) \times n\) matrix that is fixed once we condition on \(\mathbf{X}\), so that \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}= \mathbf{A}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\). Under Definition 1, \({\mathbb{E}\mathopen{}\left[\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\mid \mathbf{X}\right]\mathclose{}} = \mathbf{X}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\), and because the \(Y_i\) are independent with common variance \(\sigma^2\), \({\operatorname{Cov}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\mid \mathbf{X}\right)\mathclose{}} = \sigma^2\mathbf{I}_n\).
For the mean, using linearity of expectation for a fixed matrix times a random vector:
\[
\begin{aligned}
{\mathbb{E}\mathopen{}\left[\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\mid \mathbf{X}\right]\mathclose{}}
&= \underbrace{\mathbf{A}}_{(p+1) \times n} \underbrace{{\mathbb{E}\mathopen{}\left[\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\mid \mathbf{X}\right]\mathclose{}}}_{n \times 1}
&& \text{(linearity of expectation)}
\\
&= {\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}{\mathbf{X}}^{\top} \mathbf{X}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}
&& \text{(substitute } \mathbf{A} \text{ and } {\mathbb{E}\mathopen{}\left[\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\mid \mathbf{X}\right]\mathclose{}} \text{)}
\\
&= \underbrace{\mathbf{I}_{p+1}}_{(p+1) \times (p+1)} \underbrace{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}_{(p+1) \times 1}
&& \text{(a matrix inverse times the matrix is the identity)}
\\
&= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}
&& \text{(identity matrix)}
\end{aligned}
\]
For the covariance, using \({\operatorname{Cov}\mathopen{}\left(\mathbf{A}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\right)\mathclose{}} = \mathbf{A}{\operatorname{Cov}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\right)\mathclose{}}{\mathbf{A}}^{\top}\) for a fixed matrix \(\mathbf{A}\), and \({\mathbf{A}}^{\top} = \mathbf{X}{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}\) (because \({\mathbf{X}}^{\top}\mathbf{X}\) is symmetric, so its inverse is too):
\[
\begin{aligned}
{\operatorname{Cov}\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\mid \mathbf{X}\right)\mathclose{}}
&= \underbrace{\mathbf{A}}_{(p+1) \times n}
\underbrace{{\operatorname{Cov}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\mid \mathbf{X}\right)\mathclose{}}}_{n \times n}
\underbrace{{\mathbf{A}}^{\top}}_{n \times (p+1)}
&& \text{(covariance of a linear transformation)}
\\
&= \mathbf{A} \mathopen{}\left(\sigma^2\mathbf{I}_n\right)\mathclose{} {\mathbf{A}}^{\top}
&& \text{(substitute } {\operatorname{Cov}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\mid \mathbf{X}\right)\mathclose{}} \text{)}
\\
&= \sigma^2\mathbf{A} \mathbf{I}_n {\mathbf{A}}^{\top}
&& \text{(move the scalar } \sigma^2\text{ to the front)}
\\
&= \sigma^2\mathbf{A} {\mathbf{A}}^{\top}
&& \text{(multiplying by } \mathbf{I}_n \text{ changes nothing)}
\\
&= \sigma^2{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}{\mathbf{X}}^{\top} \mathbf{X}{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}
&& \text{(substitute } \mathbf{A} \text{ and } {\mathbf{A}}^{\top} \text{)}
\\
&= \sigma^2\mathbf{I}_{p+1} {\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}
&& \text{(a matrix inverse times the matrix is the identity)}
\\
&= \sigma^2{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}
&& \text{(identity matrix)}
\end{aligned}
\]
Example 21 (Covariance matrix of \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) for the birthweight model) For the birthweight model with an interaction term (bw_lm2), Theorem 7 says that \({\operatorname{Cov}\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\mid \mathbf{X}\right)\mathclose{}} = \sigma^2{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}\). Replacing the unknown \(\sigma^2\) by its estimate \(\hat{\sigma}^2\) gives the estimated covariance matrix that vcov() reports:
The largest difference between the two matrices is zero, up to floating-point rounding error.
To see that \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) is unbiased with this covariance matrix, we can treat the fitted coefficients and \(\hat{\sigma}\) as the true \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) and \(\sigma\), simulate many new outcome vectors \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{Y}}\vphantom{Y}}\) with the same design matrix \(\mathbf{X}\), and refit the model to each one:
The sample mean of the 10,000 simulated estimates is close to \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) for every coefficient, and the standard deviation of the simulated estimates is close to the standard error from \(\sigma^2{\mathopen{}\left({\mathbf{X}}^{\top}\mathbf{X}\right)\mathclose{}}^{-1}\), as Equation 31 predicts.
In the Gaussian linear regression case, we also have exact results. To test \(H_0: \beta_j = \beta_{j,0}\) (typically \(\beta_{j,0} = 0\)):
The covariance matrix\(\widehat{\operatorname{Cov}}(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}})\) is a \(p \times p\)symmetric matrix, where \(p\) is the number of regression coefficients (including the intercept, if present). Its rows and columns correspond to those \(p\) coefficient estimates. When an intercept is included, the coefficients are typically written \(\hat\beta_0, \hat\beta_1, \ldots\), with \(\hat\beta_0\) denoting the intercept. The matrix entries themselves are still indexed by position, so matrix row/column index \(1\) corresponds to the intercept term when one is included.
The diagonal entries are the variances of the individual coefficient estimates: the entry in matrix position \((j+1, j+1)\) is \({\operatorname{Var}\mathopen{}\left(\hat\beta_j\right)\mathclose{}}\), for \(j = 0, 1, \ldots, p - 1\).
The off-diagonal entries are the covariances between pairs of coefficient estimates: the entry in matrix position \((i+1, j+1)\) is \({\operatorname{Cov}\mathopen{}\left(\hat\beta_i, \hat\beta_j\right)\mathclose{}}\), for \(i \neq j\).
The matrix is symmetric: \({\operatorname{Cov}\mathopen{}\left(\hat\beta_i, \hat\beta_j\right)\mathclose{}} = {\operatorname{Cov}\mathopen{}\left(\hat\beta_j, \hat\beta_i\right)\mathclose{}}\).
Example 22 For model 2, which has four coefficients \((\hat\beta_0, \hat\beta_M, \hat\beta_A, \hat\beta_{AM})\), the covariance matrix has the following layout:
Table 31: Covariance matrix of \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) for birthweight model 2 (with interaction term)
Statistical significance versus practical importance
A small \(p\)-value says that an estimated slope is distinguishable from zero. It does not say that the slope is large enough to matter. With enough data, the standard error of \(\hat\beta_1\) shrinks toward zero, so even a tiny true slope becomes statistically significant.
Example 23 (A statistically significant but negligible slope) We simulate \(n = 10^6\) observations with
so the true slope is \(\beta_1 = 0.01\): a one-unit increase in \(X\), which is one standard deviation of \(X\), changes the mean of \(Y\) by 0.01 standard deviations of \(Y\)’s noise.
The estimated slope is 0.011 (standard error 0.001), with a 95% confidence interval of (0.00906, 0.013) and \(p =\) <2e-16. The test of \(H_0: \beta_1 = 0\) is rejected, yet the model explains 0.012% of the variance of \(Y\). The confidence interval tells the same story: it excludes zero, and every value in it is tiny.
Remark 3 (Significance is not importance). Statistical significance depends on the sample size as well as on the size of the effect. When we report a slope, we should report the estimate and its confidence interval, and judge whether the values in the interval are large enough to matter in the scientific setting, not only whether the interval excludes zero. With \(n = 10^4\) in place of \(n = 10^6\), the standard error of the same slope would be about 0.01, so the same true slope would usually not be significant at all.
5.3 Inference for predicted means
Exercise 18 Given a maximum likelihood estimate\(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) and a corresponding estimated covariance matrix\(\hat{\Sigma}\stackrel{\text{def}}{=}\widehat{\operatorname{Cov}}(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}})\), calculate a 95% confidence interval for the predicted mean \(\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) = {\mathbb{E}\mathopen{}\left[Y|\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}=\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right]\mathclose{}}\).
Theorem 8 (Confidence interval for a predicted mean) In a linear regression model with \(n\) observations and \(p\) coefficients, let \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) be the OLS estimate of \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\), and let \(\hat\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}} \cdot \hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) be the estimated mean at covariate pattern \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\). A 95% confidence interval for \(\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) = {\mathbb{E}\mathopen{}\left[Y|\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{X}}\vphantom{X}}=\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right]\mathclose{}}\) is:
Example 24 (Confidence interval for the mean birthweight of 40-week males) For the birthweight model with an interaction (bw_lm2), the code below applies Theorem 8 at sex male and gestational age 40 weeks, and compares the result with the interval from predict(interval = "confidence"):
By the definition \(\hat\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) and the variance of a linear combination:
where \(\Sigma \stackrel{\text{def}}{=}\operatorname{Cov}(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}})\).
\[
\begin{aligned}
{\color{blue}{{\mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})\right)\mathclose{}}}}
&= \sqrt{{\color{orange}{{\mathop{\widehat{\operatorname{Var}}}\nolimits\mathopen{}\left(\hat\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})\right)\mathclose{}}}}}
&& \text{(take the square root)}
\\
&= {\color{blue}{\sqrt{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}}^{\top}\hat{\Sigma}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}}}}
&& \text{(substitute the previous display)}
\end{aligned}
\]
Theorem 9 (Estimated variance and standard error of predicted mean)\[
{\color{orange}{{\mathop{\widehat{\operatorname{Var}}}\nolimits\mathopen{}\left(\hat\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})\right)\mathclose{}}}} = {\color{orange}{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}}^{\top}\hat{\Sigma}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}}}
\tag{34}\]
In R, predict() with se.fit = TRUE computes the estimated mean \(\hat\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) and its estimated standard error for each covariate pattern:
Why are confidence bands narrower near the center of the data?
The confidence bands in Figure 47 are visibly narrower near the center of the gestational age range. To understand why, consider a simple linear regression with one predictor \(a\) (gestational age) and an intercept. The covariate vector for a new observation at age \(a\) is \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}= {(1, a)}^{\top}\), and the general variance formula for the predicted mean specializes to:
where \(\bar{a} = \frac{1}{n}\sum_{i=1}^na_i\) is the sample mean gestational age and \(S_{AA} = \sum_{i=1}^n(a_i - \bar{a})^2\) is the total variation in gestational age.
To derive Equation 36, first note that for \(n\) observations with design matrix rows \((1, a_i)\):
Since \((a - \bar{a})^2 \geq 0\) and equals zero when \(a = \bar{a}\), the standard error is minimized at the sample mean gestational age \(\bar{a}\) and increases as \(a\) moves away from \(\bar{a}\) in either direction. Consequently, confidence intervals are narrowest at \(\bar{a}\) and widen toward the extremes of the gestational age range.
Intuitively, the fitted line is “anchored” at the center of the data: in simple linear regression with an intercept, the OLS fitted line passes exactly through the sample mean point \((\bar{a}, \bar{y})\), so the estimated mean at \(\bar{a}\) is relatively stable across different samples. Moving away from \(\bar{a}\), small changes in the estimated slope cause the fitted line to “pivot” around that anchor, producing larger changes in the predicted mean the further \(a\) is from \(\bar{a}\).
Additionally, there is typically more nearby data near the center of the covariate range, so we have more information about the true mean response there. Near the edges of the covariate range, there is less nearby data, leaving us less confident in our estimate of the mean response at those values.
5.4 Inference for differences in means
Exercise 20 Given a maximum likelihood estimate\(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) and a corresponding estimated covariance matrix\(\hat{\Sigma}\stackrel{\text{def}}{=}\widehat{\operatorname{Cov}}(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}})\), calculate a 95% confidence interval for the difference in means comparing covariate patterns \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) and \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\):
Theorem 10 (Confidence interval for a difference in means) In a linear regression model with \(n\) observations and \(p\) coefficients, let \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) and \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\) be two covariate patterns, let \(\Delta \mu \stackrel{\text{def}}{=}\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) - \mu(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}})\) be the difference in means between them, and let \(\widehat{\Delta \mu}\) be its estimate from the OLS estimate \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\). A 95% confidence interval for \(\Delta \mu\) is:
Example 25 (Confidence interval for the male-female birthweight difference at 40 weeks) For the birthweight model with an interaction (bw_lm2), the code below applies Theorem 10 to the difference between males and females at 40 weeks, using the estimated covariance matrix of the coefficients. As a check, it refits the model with age centered at 40 weeks, where the sex coefficient is exactly this difference:
The estimated difference is 136.3 grams, with 95% confidence interval (-63.6, 336.2). The refitted model gives (-63.6, 336.2), which agrees.
Exercise 21 Express the difference in means \(\Delta \mu \stackrel{\text{def}}{=}\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) - \mu(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}})\) as a linear function of \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\).
where \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\stackrel{\text{def}}{=}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}- \mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\).
Theorem 11 (Difference in means as a linear combination of coefficients) Let \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) and \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\) be two covariate patterns, and let \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\stackrel{\text{def}}{=}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}- \mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\) be their difference. In a linear regression model with coefficient vector \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\), the difference in means \(\Delta \mu \stackrel{\text{def}}{=}\mu(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}) - \mu(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}})\) is:
Example 26 (Difference in mean birthweight as a linear combination) For the birthweight model with an interaction (bw_lm2), let \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) be the covariate pattern for males and \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\) the pattern for females, both at 40 weeks. The code below computes \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\), and compares \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}} \cdot \hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\) with the difference in the two predicted means:
where \(\Sigma \stackrel{\text{def}}{=}\operatorname{Cov}(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}})\).
Theorem 12 (Estimated variance and standard error of difference in means)\[
{\color{orange}{{\mathop{\widehat{\operatorname{Var}}}\nolimits\mathopen{}\left(\widehat{\Delta \mu}\right)\mathclose{}}}} = {\color{orange}{{\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}}^{\top}\hat{\Sigma}(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})}}
\tag{41}\]
In R, we can compute the CI for the difference in means by computing the predicted means for both covariate patterns and applying the variance formula directly:
Table 36: 95% CI for difference in birthweight means (male vs. female)
diff_mean
se
ci_lower
ci_upper
136.305
95.846
-63.626
336.236
5.5 Prediction
We can also construct prediction intervals for the value of a new observation \(Y^*\), given a covariate pattern \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\).
A new observation \(Y^*\) at \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\) follows the same linear model as the training data:
where \(\varepsilon^* \sim {\operatorname{N}\mathopen{}\left(0, \sigma^2\right)\mathclose{}}\) is independent of the training data (and therefore independent of \(\hat{\beta}\)).
\({\operatorname{Var}\mathopen{}\left(\varepsilon^*\right)\mathclose{}} = \sigma^2\): the variance of the future observation’s own noise (the residual variance), which is irreducible given \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\).
\({\operatorname{Var}\mathopen{}\left(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}} \cdot \hat{\beta}\right)\mathclose{}} = {(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}})}^{\top}{\operatorname{Var}\mathopen{}\left(\hat{\beta}\right)\mathclose{}}\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}\): the variance due to estimating the mean \(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) from the data.
The first component is not present in the confidence interval for \(\hat\mu(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}})\), which only accounts for estimation uncertainty. This extra variance component is why prediction intervals are always wider than the corresponding confidence intervals: a prediction interval must additionally account for the noise \(\varepsilon^*\) in a new observation.
The standard error of the prediction error is the square root of its variance:
Since \(\varepsilon^* \sim {\operatorname{N}\mathopen{}\left(0, \sigma^2\right)\mathclose{}}\) and \(\hat{\beta}\sim {\operatorname{N}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}, \sigma^2(\mathbf{X}'\mathbf{X})^{-1}\right)\mathclose{}}\) are both normally distributed, and since \(\varepsilon^* \perp\!\!\!\perp\hat{\beta}\) (by independence of the new error term and parameter estimates), their difference \(Y^* - \hat{Y}^* = \varepsilon^* - {(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}})}^{\top}(\hat{\beta}- \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}})\) is also normally distributed, with mean 0 and variance given by Equation 43. Standardizing by dividing by \(\sigma\sqrt{1 + {(\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}})}^{\top}(\mathbf{X}'\mathbf{X})^{-1}\mathord{{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^*}}}\) yields a standard normalrandom variable. Replacing \(\sigma\) with \(\hat{\sigma}\), which is estimated with \(n-p\) degrees of freedom and is independent of \(Y^* - \hat{Y}^*\), the standardized prediction error follows a \(t\) distribution:
x <-tibble(age =40, sex ="male")bw_lm2 |>predict(newdata = x, interval ="predict")#> fit lwr upr#> 1 3210.64 2805.71 3615.57
If you don’t specify newdata, you get a warning:
[R code]
bw_lm2 |>predict(interval ="predict") |>head()#> Warning in predict.lm(bw_lm2, interval = "predict"): predictions on current data refer to _future_ responses#> fit lwr upr#> 1 2552.73 2124.50 2980.97#> 2 2552.73 2124.50 2980.97#> 3 2683.13 2275.99 3090.27#> 4 2813.53 2418.60 3208.47#> 5 2813.53 2418.60 3208.47#> 6 2943.93 2551.48 3336.38
The warning from the last command is: “predictions on current data refer to future responses” (since you already know what happened to the current data, and thus don’t need to predict it).
Figure 48: Confidence and prediction intervals for birthweight model 2
6 Multiple linear regression: the HERS study
The birthweight example carried a single continuous predictor (gestational age) alongside a binary group (sex), so every fitted model was a line — or a pair of lines — in a plane.
To see how the same ideas extend to several continuous predictors at once, we now work through a larger example: data from the Heart and Estrogen/Progestin Replacement Study (HERS).
hers_scatter_data <- hers |> dplyr::filter(!is.na(age), !is.na(BMI), !is.na(LDL), !is.na(statins))hers_scatter_no <- hers_scatter_data |> dplyr::filter(statins =="no")hers_scatter_yes <- hers_scatter_data |> dplyr::filter(statins =="yes")hers_scatter_axes <-list(xaxis =list(title ="Age (yr)"),yaxis =list(title ="BMI (kg/m²)"),zaxis =list(title ="LDL (mg/dL)"))plotly::plot_ly() |> plotly::add_trace(data = hers_scatter_no,x =~age, y =~BMI, z =~LDL,type ="scatter3d", mode ="markers",name ="No statins", scene ="scene",marker =list(size =3, opacity =0.5, color ="steelblue") ) |> plotly::add_trace(data = hers_scatter_yes,x =~age, y =~BMI, z =~LDL,type ="scatter3d", mode ="markers",name ="Yes statins", scene ="scene2",marker =list(size =3, opacity =0.5, color ="darkorange") ) |> plotly::layout(annotations =list(list(text ="No statins", x =0.22, y =0.95,xref ="paper", yref ="paper", showarrow =FALSE ),list(text ="Yes statins", x =0.78, y =0.95,xref ="paper", yref ="paper", showarrow =FALSE ) ),scene =c( hers_scatter_axes,list(domain =list(x =c(0, 0.48), y =c(0, 1))) ),scene2 =c( hers_scatter_axes,list(domain =list(x =c(0.52, 1), y =c(0, 1))) ) )
Figure 49: hers data (Vittinghoff et al. [3]): age (years) and BMI (kg/m²) vs. baseline LDL (mg/dL), split by statin use (left: no statins; right: statins).
Figure 52: Fitted regression surface from Model 46 (patients taking statins)
The coefficient estimates in Table 38 quantify the shape of these surfaces. Both tabs display the same fitted flat plane: the parallel-planes model constrains the age and BMI slopes to be the same regardless of statin use. The statins-yes tab reveals that patients taking statins tend to have lower LDL values on average — a shift the parallel-planes model cannot capture, because it does not include statin use as a predictor.
6.3 Interaction model
\[
\begin{aligned}
Y|A,B &\ \sim_{\operatorname{ciid}}\ N(\mu(A,B), \sigma^2)\\
\mu(a,b) &= \beta_0 + \beta_A a + \beta_B b + \beta_{AB}(a \cdot b)
\end{aligned}
\tag{47}\]
The slope of mean LDL with respect to BMI now depends on age:
Figure 54: Fitted regression surface from Model 47 (patients taking statins)
Unlike the flat parallel-planes surface, the interaction model produces a twisted regression surface. The BMI slope changes with age: at any given age, the BMI slope equals the base BMI coefficient plus the interaction coefficient times that age. Both tabs plot the same interaction-model surface; statin use is used only to partition the observed data points for visual clarity.
Equivalently, we can fit the same model using reference (indicator) coding, with an overall intercept, a main effect for statin use, a BMI slope, and a statin–BMI interaction:
\[
{\mathbb{E}\mathopen{}\left[LDL|B=b, V=v\right]\mathclose{}} =
\beta_0 + \beta_V\, v + \beta_B\, b + \beta_{VB}\, (v \cdot b)
\tag{49}\]
Table 44: Likelihood ratio test comparing the hers parallel-planes (hers_lm1) and interaction (hers_lm2) models
Here, the likelihood ratio test is not statistically significant (Table 44: \(\chi^2_1 \approx 1.4\), \(p \approx 0.24\)): the age-BMI interaction does not fit the data meaningfully better than we would expect from adding one parameter by chance. So we prefer the simpler parallel-planes model (hers_lm1, Equation 46) on parsimony grounds, consistent with the AIC and BIC comparison in Table 42.
7 End-of-chapter exercises
These paper-and-pencil exercises use actual datasets. No numerical computation is required.
Exercise 23 Main-effects model interpretation.
Adapted from the main-effects interpretation exercise in [3, Ch. 4].
Dataset summary: The ToothGrowth dataset contains one row per guinea pig. Variable definitions:
\(Y\): tooth length (len)
\(X\): dose (dose)
\(Z\): supplement indicator (supp) with \(Z=1\) for VC and \(Z=0\) for OJ
Consider the additive model
\[
{\mathbb{E}\mathopen{}\left[Y|X=x,Z=z\right]\mathclose{}}
= \beta_0 + \beta_X x + \beta_Z z.
\]
State the interpretation of \(\beta_0\), \(\beta_X\), and \(\beta_Z\). Which interpretations depend on the reference levelOJ?
Figure 59: ToothGrowth data with additive (parallel-lines) fit
Solution.
\(\beta_0\) is the mean tooth length for guinea pigs given OJ at dose 0.
\(\beta_X\) is the change in mean tooth length for each one-unit increase in dose. In this additive model, that dose effect is the same for OJ and VC.
\(\beta_Z\) is the mean tooth-length difference between VC and OJ at the same dose. The interpretations that depend on OJ as the reference level are \(\beta_0\) and \(\beta_Z\).
Exercise 24 Interaction interpretation.
Adapted from the interaction-slope interpretation example in [1, Sec. 6.7.1].
Use ToothGrowth as introduced in Exercise 23. Variable definitions:
\(Y\): tooth length (len)
\(X\): dose (dose)
\(Z\): supplement indicator (supp) with \(Z=1\) for VC and \(Z=0\) for OJ
Consider the interaction model
\[
{\mathbb{E}\mathopen{}\left[Y|X=x,Z=z\right]\mathclose{}}
= \beta_0 + \beta_X x + \beta_Z z + \beta_{XZ}(x\cdot z).
\]
Using the plot, write the slope with respect to \(X\) for OJ and for VC. Then interpret \(\beta_{XZ}\).
eoc_mtcars_fit <-lm(mpg ~ wt + hp, data = eoc_mtcars)autoplot(eoc_mtcars_fit, which =1, ncol =1)autoplot(eoc_mtcars_fit, which =2, ncol =1)
Figure 62: Diagnostic plots for linear model fit to mtcars data
(a) Residuals versus fitted
(b) Normal Q-Q
Solution.
The residual-vs-fitted plot shows nonlinearity and non-constant spread. The Q-Q plot shows mild tail departures from normality.
The correct conclusion is that linearity and constant-variance assumptions are the most questionable here. A practical next step is to fit a model with a nonlinear weight term (for example, a spline or polynomial in wt) and then re-check diagnostics.
Exercises
Exercise 27 (Likelihood for simple linear regression) (adapted from Kleinbaum et al. [7], Chapter 6, and Dobson and Barnett [1], Chapter 6)
(a) Write the likelihood \(\mathcal{L}(\beta_0, \beta_1, \sigma^2; \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}, \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})\) for observed data \((y_i, x_i)\), \(i = 1, \ldots, n\).
(b) Write the log-likelihood\(\ell(\beta_0, \beta_1, \sigma^2; \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}, \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}})\) and show it simplifies to
Since \(-\frac{n}{2}{\operatorname{log}\mathopen{}\left\{2\pi\sigma^2\right\}\mathclose{}}\) does not depend on \(\beta_0\) or \(\beta_1\), and \(\sigma^2 > 0\), maximizing \(\ell\) over \((\beta_0, \beta_1)\) is equivalent to maximizing \(-\frac{\text{RSS}}{2\sigma^2}\), which is equivalent to minimizing \(\text{RSS}\).
Exercise 28 (Score equations for simple linear regression) (adapted from Dobson and Barnett [1], Chapter 6)
Using the log-likelihood from Exercise 27, derive the score equations for \(\beta_0\) and \(\beta_1\). Specifically, compute \(\frac{\partial \ell}{\partial \beta_0} = 0\) and \(\frac{\partial \ell}{\partial \beta_1} = 0\), and show they are equivalent to the normal equations:
These score equations are the normal equations, whose solution gives the OLS/MLE estimators\(\hat\beta_0\) and \(\hat\beta_1\).
Exercise 29 (Closed-form OLS estimators (adapted from Vittinghoff et al. [3], Chapter 4)) Using the normal equations from Exercise 28, show that the closed-form solutions are:
where \(Y\) is diastolic blood pressure (mmHg), \(X_1\) is age (years), and \(X_2\) is an indicator for current smoking (\(X_2 = 1\) for smoker, \(X_2 = 0\) for non-smoker). Suppose the estimates are \(\hat\beta_0 = 55\), \(\hat\beta_1 = 0.4\), \(\hat\beta_2 = 6\).
(a) Interpret \(\hat\beta_0\).
(b) Interpret \(\hat\beta_1\).
(c) Interpret \(\hat\beta_2\).
(d) Predict mean blood pressure for a 50-year-old smoker.
(e) What assumption does this model make about the relationship between age and blood pressure for smokers versus non-smokers?
Solution. (a)
\(\hat\beta_0 = 55\) mmHg is the estimated mean diastolic blood pressure among non-smokers (\(X_2 = 0\)) aged zero (\(X_1 = 0\)). This prediction is an extrapolation beyond the observed range of ages and may not be a meaningful quantity on its own.
(b)
\(\hat\beta_1 = 0.4\) mmHg/year is the estimated mean increase in diastolic blood pressure per one-year increase in age, holding smoking status constant. Or equivalently: among people of the same smoking status, those who are one year older have, on average, 0.4 mmHg higher diastolic blood pressure.
(c)
\(\hat\beta_2 = 6\) mmHg is the estimated mean difference in diastolic blood pressure between current smokers and non-smokers of the same age. Smokers have, on average, 6 mmHg higher diastolic blood pressure than non-smokers of the same age.
The model assumes that the slope of mean blood pressure with respect to age is the same (\(\hat\beta_1 = 0.4\) mmHg/year) for both smokers and non-smokers (parallel slopes or no interaction). An interaction term \(\beta_3 X_1 X_2\) would be needed to allow different slopes.
Exercise 31 (Coefficient of determination) (adapted from Kleinbaum et al. [7], Chapter 7)
(a) Explain in words what \(\text{TSS}\), \(\text{RSS}\), and \(R^2\) measure.
(b) If \(n = 20\), \(\text{TSS} = 400\), and \(R^2 = 0.75\), compute \(\text{RSS}\).
(c) Is a large \(R^2\) sufficient to conclude the model is correct? Briefly explain.
Solution. (a)
\(\text{TSS}\) (total sum of squares): the total variability in \(Y\) around its sample mean \(\bar{y}\), ignoring the predictor \(X\).
\(\text{RSS}\) (residual sum of squares): the residual variability in \(Y\) remaining after accounting for the linear regression on \(X\). Smaller values indicate a better-fitting model.
\(R^2\): the proportion of the total variability in \(Y\) explained by the linear regression on \(X\). \(R^2 = 0\) means \(X\) explains none of the variation; \(R^2 = 1\) means \(X\) perfectly predicts \(Y\).
No. A high \(R^2\) indicates good fit in terms of variance explained, but does not guarantee the model is correctly specified. For example:
The true relationship may be nonlinear, and a high \(R^2\) can still be achieved with a linear model over a limited range.
Important confounders or interaction terms may be omitted.
The residuals may violate the model assumptions (non-normality, heteroscedasticity, dependence), which would invalidate inference even with high \(R^2\).
References
[1]
A. J. Dobson and A. G. Barnett, An introduction to generalized linear models, 4th ed. CRC press, 2018. Available: https://doi.org/10.1201/9781315182780
E. Vittinghoff, D. V. Glidden, S. C. Shiboski, and C. E. McCulloch, Regression methods in biostatistics: Linear, logistic, survival, and repeated measures models, 2nd ed. Springer, 2012. Available: https://doi.org/10.1007/978-1-4614-1353-0
[4]
M. H. Kutner, C. J. Nachtsheim, J. Neter, and W. Li, Applied linear statistical models. McGraw-Hill, 2005.
S. Weisberg, Applied linear regression, vol. 528. John Wiley & Sons, 2005.
[12]
R. McElreath, Statistical rethinking : A Bayesian course with examples in R and Stan, Second edition. in Chapman & hall/CRC texts in statistical science series. Boca Raton, FL: CRC Press, 2020.
J. Starmer, “StatQuest with JoshStarmer.” YouTube channel. Available: https://statquest.org/
[17]
H. Akaike, “A new look at the statistical model identification,”IEEE Transactions on Automatic Control, vol. 19, no. 6, pp. 716–723, 1974, Available: https://doi.org/10.1109/TAC.1974.1100705
[18]
G. Schwarz, “Estimating the dimension of a model,”The Annals of Statistics, vol. 6, no. 2, pp. 461–464, 1978, Available: https://doi.org/10.1214/aos/1176344136
[19]
R. A. Polin, W. W. Fox, and S. H. Abman, Fetal and neonatal physiology, 4th ed. Elsevier health sciences, 2011.
[20]
S. Hulley et al., “Randomized trial of estrogen plus progestin for secondary prevention of coronary heart disease in postmenopausal women,”JAMA : the journal of the American Medical Association, vol. 280, no. 7, pp. 605–613, 1998.
[21]
T. S. Breusch and A. R. Pagan, “A simple test for heteroscedasticity and random coefficient variation,”Econometrica, vol. 47, pp. 1287–1294, 1979.
[22]
R. Koenker, “A note on studentizing a test for heteroscedasticity,”Journal of Econometrics, vol. 17, pp. 107–112, 1981.
[23]
G. James, D. Witten, T. Hastie, and R. Tibshirani, An introduction to statistical learning: With applications in R, 2nd ed. Springer, 2021. doi: 10.1007/978-1-0716-1418-1.
[24]
R. D. Cook, “Detection of influential observation in linear regression,”Technometrics, vol. 19, no. 1, pp. 15–18, 1977, doi: 10.1080/00401706.1977.10489493.
[25]
K. A. Bollen and R. W. Jackman, “Regression diagnostics: An expository treatment of outliers and influential cases,” in Modern methods of data analysis, J. Fox and J. S. Long, Eds., Newbury Park, CA: Sage, 1990, pp. 257–291.
G. James, D. Witten, T. Hastie, R. Tibshirani, et al., An introduction to statistical learning, vol. 112. Springer, 2013. Available: https://www.statlearning.com/
[28]
F. E. Harrell, Regression modeling strategies: With applications to linear models, logistic regression, and survival analysis, 2nd ed. Springer, 2015. doi: 10.1007/978-3-319-19425-7.
[29]
G. Heinze, C. Wallisch, and D. Dunkler, “Variable selection – A review and recommendations for the practicing statistician,”Biometrical Journal, vol. 60, no. 3, pp. 431–449, 2018, doi: 10.1002/bimj.201700067.
[30]
R. V. Hogg, E. A. Tanis, and D. L. Zimmerman, Probability and statistical inference, Ninth edition. Boston: Pearson, 2015.