Model Selection
EPI 553 — Principles of Statistical Inference II (Spring 2026)
Emmanuel Nana Arko
March 24, 2026
Introduction
In previous lectures, we learned how to fit multiple linear regression models, include dummy variables for categorical predictors, test for interactions, and assess confounding. But we have not yet addressed a fundamental question: how do we decide which variables belong in the model?
This question has different answers depending on the goal of the analysis:
| Goal | What matters | Variable selection driven by |
|---|---|---|
| Prediction | Model accuracy and reliability in new data | Statistical criteria (Adj. \(R^2\), AIC, BIC, cross-validation) |
| Association | Validity of the exposure coefficient | Subject-matter knowledge, confounding assessment, 10% rule |
In predictive modeling, we search for the subset of variables that best predicts \(Y\) without overfitting. In associative modeling, the exposure variable is always in the model, and we decide which covariates to include based on whether they are confounders.
This lecture covers both approaches, with emphasis on when each is appropriate and the pitfalls of automated selection.
Setup and Data
library(tidyverse)
library(haven)
library(janitor)
library(knitr)
library(kableExtra)
library(broom)
library(gtsummary)
library(car)
library(leaps)
library(MASS)
options(gtsummary.use_ftExtra = TRUE)
set_gtsummary_theme(theme_gtsummary_compact(set_theme = TRUE))The BRFSS 2020 Dataset
We continue with the BRFSS 2020 dataset, predicting physically unhealthy days from a pool of candidate predictors.
# Load the pre-saved analytic dataset (use the RDS file in your working directory)
brfss_ms <- readRDS("/Users/emmanuelsmac/Desktop/brfss_ms_2020.rds")# Dataset already cleaned and saved; display dimensions only
tibble(Metric = c("Observations", "Variables"),
Value = c(nrow(brfss_ms), ncol(brfss_ms))) |>
kable(caption = "Analytic Dataset Dimensions") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Metric | Value |
|---|---|
| Observations | 5000 |
| Variables | 10 |
We have 10 variables: 1 outcome and 9 candidate predictors. If we considered all possible subsets of the 9 predictors (ignoring interactions and transformations), there would be \(2^9 - 1 = 511\) possible models.
Part 1: Guided Practice — Model Selection
1. Building the Maximum Model
1.1 What Is the Maximum Model?
The maximum model is the model that includes all candidate predictor variables. It represents the upper bound of complexity. The “correct” model will have \(p \leq k\) predictors, where \(k\) is the number in the maximum model.
The candidate variables in the maximum model can include:
- Main effects (e.g., age, sex, BMI)
- Higher-order terms (e.g., age\(^2\), age\(^3\))
- Transformations (e.g., log(BMI))
- Interactions (e.g., sex \(\times\) age)
These candidates are chosen based on a literature search and the research question, not by throwing in every available variable.
# The maximum model with all candidate predictors
mod_max <- lm(physhlth_days ~ menthlth_days + sleep_hrs + age + sex +
education + exercise + gen_health + income_cat + bmi,
data = brfss_ms)
tidy(mod_max, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Maximum Model: All Candidate Predictors",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 2.6902 | 0.8556 | 3.1441 | 0.0017 | 1.0128 | 4.3676 |
| menthlth_days | 0.1472 | 0.0121 | 12.1488 | 0.0000 | 0.1235 | 0.1710 |
| sleep_hrs | -0.1930 | 0.0673 | -2.8679 | 0.0041 | -0.3249 | -0.0611 |
| age | 0.0180 | 0.0055 | 3.2969 | 0.0010 | 0.0073 | 0.0288 |
| sexFemale | -0.1889 | 0.1820 | -1.0376 | 0.2995 | -0.5458 | 0.1680 |
| educationHS graduate | 0.2508 | 0.4297 | 0.5836 | 0.5595 | -0.5917 | 1.0933 |
| educationSome college | 0.3463 | 0.4324 | 0.8009 | 0.4233 | -0.5014 | 1.1940 |
| educationCollege graduate | 0.3336 | 0.4357 | 0.7657 | 0.4439 | -0.5206 | 1.1878 |
| exerciseYes | -1.2866 | 0.2374 | -5.4199 | 0.0000 | -1.7520 | -0.8212 |
| gen_healthVery good | 0.4373 | 0.2453 | 1.7824 | 0.0747 | -0.0437 | 0.9183 |
| gen_healthGood | 1.5913 | 0.2651 | 6.0022 | 0.0000 | 1.0716 | 2.1111 |
| gen_healthFair | 7.0176 | 0.3682 | 19.0586 | 0.0000 | 6.2957 | 7.7394 |
| gen_healthPoor | 20.4374 | 0.5469 | 37.3722 | 0.0000 | 19.3653 | 21.5095 |
| income_cat | -0.1817 | 0.0503 | -3.6092 | 0.0003 | -0.2803 | -0.0830 |
| bmi | 0.0130 | 0.0145 | 0.8997 | 0.3683 | -0.0153 | 0.0414 |
glance(mod_max) |>
dplyr::select(r.squared, adj.r.squared, sigma, AIC, BIC, df.residual) |>
mutate(across(everything(), \(x) round(x, 3))) |>
kable(caption = "Maximum Model: Fit Statistics") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| r.squared | adj.r.squared | sigma | AIC | BIC | df.residual |
|---|---|---|---|---|---|
| 0.386 | 0.384 | 6.321 | 32645.79 | 32750.06 | 4985 |
Interpretation: The maximum model explains approximately 38.6% of the variance in physically unhealthy days (R² = 0.386, Adjusted R² = 0.384). The strongest predictors are general health status (with “Poor” health associated with about 20 more unhealthy days compared to “Excellent”) and mental health days (each additional mentally unhealthy day is associated with 0.15 more physically unhealthy days). Exercise is also strongly associated, with exercisers reporting about 1.3 fewer physically unhealthy days. Several variables, including sex (p = 0.30), education (p > 0.40 for all levels), and BMI (p = 0.37), are not statistically significant, suggesting they may be candidates for removal in a more parsimonious model. The AIC is 32,645.8 and BIC is 32,750.1; these serve as baselines for comparing simpler models.
1.2 Overfitting vs. Underfitting
The goal of model building is to find the right balance:
| Problem | What happens | Consequence |
|---|---|---|
| Overfitting | Including variables with \(\beta = 0\) | No bias, but increased collinearity, inflated SEs, poor out-of-sample prediction |
| Underfitting | Omitting variables with \(\beta \neq 0\) | Bias in the remaining coefficients (omitted variable bias) |
Key insight: Underfitting is worse than overfitting in terms of bias. An overfit model gives unbiased estimates (just less precise), while an underfit model gives biased estimates. However, for prediction, overfitting degrades out-of-sample performance.
The objective is a parsimonious model: the simplest model that captures the important relationships without unnecessary complexity.
1.3 How Many Predictors Can We Include?
The error degrees of freedom must be positive: \(n - k - 1 > 0\), meaning \(n > k + 1\).
Rules of thumb for the minimum sample size:
| Rule | Requirement | Our data (n = 5,000) |
|---|---|---|
| Minimum 10 error df | \(n \geq k + 11\) | Can include up to 4,989 predictors |
| 5 observations per predictor | \(n \geq 5k\) | Can include up to 1,000 predictors |
| 10 observations per predictor | \(n \geq 10k\) | Can include up to 500 predictors |
With \(n = 5,000\), we are well within all rules of thumb for our 9 candidate predictors (plus dummy variables).
Caution with categorical variables: A categorical predictor with \(k\) levels uses \(k - 1\) degrees of freedom, not just 1. Our education (4 levels) uses 3 df, gen_health (5 levels) uses 4 df, so the maximum model actually uses 14 predictor df.
2. Selection Criteria
Given a set of candidate models, we need a criterion to compare them. We cover five: \(R^2\), Adjusted \(R^2\), \(F_p\) (partial F-test), AIC, and BIC.
2.1 \(R^2\) (Coefficient of Determination)
\[R^2 = 1 - \frac{SSE}{SST} = \frac{SSR}{SST}\]
\(R^2\) measures the proportion of variance in \(Y\) explained by the model. However, \(R^2\) always increases (or stays the same) when you add a predictor, regardless of whether it is useful. This makes raw \(R^2\) useless for model comparison across models of different sizes.
# Demonstrate that R2 always increases
models <- list(
"Sleep only" = lm(physhlth_days ~ sleep_hrs, data = brfss_ms),
"+ age" = lm(physhlth_days ~ sleep_hrs + age, data = brfss_ms),
"+ sex" = lm(physhlth_days ~ sleep_hrs + age + sex, data = brfss_ms),
"+ education" = lm(physhlth_days ~ sleep_hrs + age + sex + education, data = brfss_ms),
"+ exercise" = lm(physhlth_days ~ sleep_hrs + age + sex + education + exercise, data = brfss_ms),
"+ gen_health" = lm(physhlth_days ~ sleep_hrs + age + sex + education + exercise + gen_health, data = brfss_ms),
"+ mental health" = lm(physhlth_days ~ sleep_hrs + age + sex + education + exercise + gen_health + menthlth_days, data = brfss_ms),
"+ income" = lm(physhlth_days ~ sleep_hrs + age + sex + education + exercise + gen_health + menthlth_days + income_cat, data = brfss_ms),
"+ BMI (full)" = lm(physhlth_days ~ sleep_hrs + age + sex + education + exercise + gen_health + menthlth_days + income_cat + bmi, data = brfss_ms)
)
r2_table <- map_dfr(names(models), \(name) {
g <- glance(models[[name]])
tibble(
Model = name,
p = length(coef(models[[name]])) - 1,
`R²` = round(g$r.squared, 4),
`Adj. R²` = round(g$adj.r.squared, 4),
AIC = round(g$AIC, 1),
BIC = round(g$BIC, 1)
)
})
r2_table |>
kable(caption = "Model Comparison: R² Always Increases as Predictors Are Added") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | p | R² | Adj. R² | AIC | BIC |
|---|---|---|---|---|---|
| Sleep only | 1 | 0.0115 | 0.0113 | 35001.0 | 35020.6 |
|
2 | 0.0280 | 0.0276 | 34918.7 | 34944.8 |
|
3 | 0.0280 | 0.0274 | 34920.7 | 34953.3 |
|
6 | 0.0440 | 0.0428 | 34843.7 | 34895.9 |
|
7 | 0.0849 | 0.0836 | 34626.8 | 34685.5 |
|
11 | 0.3650 | 0.3636 | 32807.7 | 32892.4 |
|
12 | 0.3843 | 0.3828 | 32655.4 | 32746.6 |
|
13 | 0.3859 | 0.3843 | 32644.6 | 32742.4 |
|
14 | 0.3860 | 0.3843 | 32645.8 | 32750.1 |
Interpretation: Notice that R² increases monotonically from 0.012 (sleep only) to 0.386 (full model) as each predictor is added. However, Adjusted R² tells a different story: it plateaus at 0.384 after adding income (the 8th predictor), and adding BMI does not improve it further (still 0.384). The largest single jump in both R² and Adjusted R² occurs when general health is added (from 0.084 to 0.365), indicating it is by far the most powerful predictor. AIC and BIC both decrease sharply at that same step. AIC reaches its minimum at the full model (32,645.8), while BIC, which penalizes complexity more heavily, favors a slightly smaller model. This table illustrates a key lesson: R² will always reward you for adding variables, even useless ones, making it unreliable for model comparison.
2.2 Adjusted \(R^2\)
Adjusted \(R^2\) penalizes for model complexity:
\[R^2_{adj} = 1 - \frac{(n - i)(1 - R^2)}{n - p}\]
where \(i = 1\) if the model includes an intercept, \(n\) is the sample size, and \(p\) is the number of predictors. Unlike \(R^2\), Adjusted \(R^2\) can decrease when an uninformative predictor is added, because the penalty for using an extra degree of freedom outweighs the tiny increase in \(R^2\).
Selection rule: Choose the model with the largest Adjusted \(R^2\).
2.3 \(F_p\) (Partial F-Test)
The partial F-test compares a reduced model (with \(p\) predictors) to the maximum model (with \(k\) predictors):
\[F_p = \frac{\{SSE(p) - SSE(k)\} / (k - p)}{SSE(k) / (n - k - 1)}\]
This tests \(H_0\): the \(k - p\) omitted variables all have \(\beta = 0\).
- If \(F_p\) is not significant, the reduced model is adequate (the extra variables are not needed)
- If \(F_p\) is significant, at least one of the omitted variables is important
Selection rule: Choose the smallest model for which \(F_p\) is not significant when compared to the maximum model.
# Compare a small model to the maximum model
mod_small <- lm(physhlth_days ~ menthlth_days + gen_health + exercise, data = brfss_ms)
anova(mod_small, mod_max) |>
tidy() |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(caption = "Partial F-test: Small Model vs. Maximum Model") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| term | df.residual | rss | df | sumsq | statistic | p.value |
|---|---|---|---|---|---|---|
| physhlth_days ~ menthlth_days + gen_health + exercise | 4993 | 200472.4 | NA | NA | NA | NA |
| physhlth_days ~ menthlth_days + sleep_hrs + age + sex + education + exercise + gen_health + income_cat + bmi | 4985 | 199201.8 | 8 | 1270.601 | 3.9746 | 1e-04 |
Interpretation: The partial F-test compares the small model (mental health days + general health + exercise) to the maximum model (all 9 predictors). The F-statistic is 3.97 with p < 0.001, meaning the null hypothesis that the additional 6 variables all have β = 0 is rejected. In other words, at least one of the omitted variables (sleep, age, sex, education, income, BMI) contributes significantly to the model beyond the three core predictors. This means the small model, despite capturing much of the explained variance, is missing important information. We should look for a model between the small and maximum that retains the significant predictors while dropping the uninformative ones.
2.4 AIC (Akaike Information Criterion)
\[AIC = 2k - 2\log(\hat{L})\]
where \(k\) is the number of estimated parameters and \(\hat{L}\) is the maximized likelihood. AIC measures the relative information lost by a model. It balances goodness of fit against complexity.
- AIC is not a test; it is a relative comparison tool
- The model with the smallest AIC is preferred
- AIC differences < 2 suggest models are essentially equivalent
Selection rule: Choose the model with the smallest AIC.
2.5 BIC (Bayesian Information Criterion)
\[BIC = k \log(n) - 2\log(\hat{L})\]
BIC is similar to AIC but penalizes complexity more heavily, especially with large sample sizes (\(\log(n)\) vs. 2). BIC tends to select simpler models than AIC.
Selection rule: Choose the model with the smallest BIC.
2.6 MSE(p) (Mean Squared Error)
\[MSE(p) = \frac{SSE_p}{n - p - 1}\]
MSE(p) is the residual variance for a model with \(p\) predictors. It balances fit (smaller SSE) against model size (fewer df in the denominator).
Selection rule: Choose the model with the smallest MSE(p).
2.7 Comparing the Criteria
| Criterion | Direction | Penalizes | Best for |
|---|---|---|---|
| R² | Maximize | No | Never use alone |
| Adjusted R² | Maximize | Yes (df penalty) | Comparing nested models |
| Fp (partial F) | Not significant → keep reduced | Yes (F distribution) | Comparing to maximum model |
| AIC | Minimize | Yes (2k) | General comparison |
| BIC | Minimize | Yes (k log n) | Favors simpler models |
| MSE(p) | Minimize | Yes (df in denominator) | Similar to Adj. R² |
criteria_long <- r2_table |>
dplyr::select(Model, p, AIC, BIC) |>
pivot_longer(cols = c(AIC, BIC), names_to = "Criterion", values_to = "Value") |>
mutate(Model = factor(Model, levels = r2_table$Model))
ggplot(criteria_long, aes(x = p, y = Value, color = Criterion)) +
geom_line(linewidth = 1.1) +
geom_point(size = 3) +
labs(
title = "AIC and BIC Across Sequentially Larger Models",
subtitle = "Lower is better; BIC penalizes complexity more heavily",
x = "Number of Predictor Degrees of Freedom (p)",
y = "Criterion Value"
) +
theme_minimal(base_size = 13) +
scale_color_brewer(palette = "Set1")
AIC and BIC Across Sequential Models
Interpretation: Both AIC and BIC decrease sharply as the first several predictors are added, with the steepest drop occurring when general health enters the model. AIC continues to decrease (or remains flat) through the full model, suggesting it favors retaining most predictors. BIC, by contrast, reaches its minimum earlier and then begins to increase, reflecting its heavier penalty for model complexity. The divergence between AIC and BIC is typical in large samples: AIC tends to select larger models, while BIC favors parsimony. In practice, when AIC and BIC disagree, the choice depends on the modeling goal: AIC is better for prediction (it minimizes information loss), while BIC is better for identifying the “true” model (it is consistent, meaning it selects the correct model as n grows).
3. Variable Selection Strategies
3.1 All Possible Regressions (Best Subsets)
The most thorough approach is to fit every possible subset of predictors and compare them. With \(k\) predictors, there are \(2^k - 1\) models.
This is computationally feasible for moderate \(k\) (up to about 20-30 predictors). In R,
the leaps package implements this efficiently:
# Prepare a model matrix (need numeric predictors for leaps)
# Use the formula interface approach
best_subsets <- regsubsets(
physhlth_days ~ menthlth_days + sleep_hrs + age + sex + education +
exercise + gen_health + income_cat + bmi,
data = brfss_ms,
nvmax = 15, # maximum number of variables to consider
method = "exhaustive"
)
best_summary <- summary(best_subsets)subset_metrics <- tibble(
p = 1:length(best_summary$adjr2),
`Adj. R²` = best_summary$adjr2,
BIC = best_summary$bic,
Cp = best_summary$cp
)
p1 <- ggplot(subset_metrics, aes(x = p, y = `Adj. R²`)) +
geom_line(linewidth = 1, color = "steelblue") +
geom_point(size = 3, color = "steelblue") +
geom_vline(xintercept = which.max(best_summary$adjr2),
linetype = "dashed", color = "tomato") +
labs(title = "Adjusted R² by Model Size", x = "Number of Variables", y = "Adjusted R²") +
theme_minimal(base_size = 12)
p2 <- ggplot(subset_metrics, aes(x = p, y = BIC)) +
geom_line(linewidth = 1, color = "steelblue") +
geom_point(size = 3, color = "steelblue") +
geom_vline(xintercept = which.min(best_summary$bic),
linetype = "dashed", color = "tomato") +
labs(title = "BIC by Model Size", x = "Number of Variables", y = "BIC") +
theme_minimal(base_size = 12)
gridExtra::grid.arrange(p1, p2, ncol = 2)
Best Subsets: Adjusted R² and BIC by Model Size
## Best model by Adj. R²: 10 variables## Best model by BIC: 8 variables# Show which variables are in the BIC-best model
best_bic_idx <- which.min(best_summary$bic)
best_vars <- names(which(best_summary$which[best_bic_idx, -1]))
cat("\nVariables in BIC-best model:\n")##
## Variables in BIC-best model:## menthlth_days
## sleep_hrs
## age
## exerciseYes
## gen_healthGood
## gen_healthFair
## gen_healthPoor
## income_catInterpretation: The best subsets analysis confirms what the sequential analysis suggested. Adjusted R² reaches its maximum at 10 variables and plateaus, while BIC selects a more parsimonious model with 8 variables. The BIC-best model retains mental health days, sleep hours, age, exercise, three levels of general health (Good, Fair, Poor), and income. Notably, it drops sex, education, Very Good health (combining it implicitly with Excellent as the reference pattern), and BMI. These are exactly the variables that had the largest p-values in the maximum model. The fact that both criteria converge on a similar core set of predictors (mental health, general health, exercise) gives us confidence that these are the genuinely important variables.
3.2 Backward Elimination
Backward elimination starts with the maximum model and removes variables one at a time:
- Fit the maximum model
- Identify the predictor with the highest p-value (smallest partial F-statistic)
- If its p-value exceeds \(\alpha\) (typically 0.05 or 0.10), remove it
- Refit the model and repeat steps 2-3
- Stop when all remaining variables have p-values \(\leq \alpha\)
## === BACKWARD ELIMINATION ===## Step 1: Maximum model## Variables: menthlth_days, sleep_hrs, age, sexFemale, educationHS graduate, educationSome college, educationCollege graduate, exerciseYes, gen_healthVery good, gen_healthGood, gen_healthFair, gen_healthPoor, income_cat, bmi# Show p-values for the maximum model
pvals <- tidy(mod_back) |>
filter(term != "(Intercept)") |>
arrange(desc(p.value)) |>
dplyr::select(term, estimate, p.value) |>
mutate(across(where(is.numeric), \(x) round(x, 4)))
pvals |>
head(5) |>
kable(caption = "Maximum Model: Variables Sorted by p-value (Highest First)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| term | estimate | p.value |
|---|---|---|
| educationHS graduate | 0.2508 | 0.5595 |
| educationCollege graduate | 0.3336 | 0.4439 |
| educationSome college | 0.3463 | 0.4233 |
| bmi | 0.0130 | 0.3683 |
| sexFemale | -0.1889 | 0.2995 |
In R, the step() function automates backward elimination
using AIC:
# Automated backward elimination using AIC
mod_backward <- step(mod_max, direction = "backward", trace = 1)## Start: AIC=18454.4
## physhlth_days ~ menthlth_days + sleep_hrs + age + sex + education +
## exercise + gen_health + income_cat + bmi
##
## Df Sum of Sq RSS AIC
## - education 3 29 199231 18449
## - bmi 1 32 199234 18453
## - sex 1 43 199245 18454
## <none> 199202 18454
## - sleep_hrs 1 329 199530 18461
## - age 1 434 199636 18463
## - income_cat 1 521 199722 18466
## - exercise 1 1174 200376 18482
## - menthlth_days 1 5898 205100 18598
## - gen_health 4 66437 265639 19886
##
## Step: AIC=18449.13
## physhlth_days ~ menthlth_days + sleep_hrs + age + sex + exercise +
## gen_health + income_cat + bmi
##
## Df Sum of Sq RSS AIC
## - bmi 1 32 199262 18448
## - sex 1 40 199270 18448
## <none> 199231 18449
## - sleep_hrs 1 327 199557 18455
## - age 1 439 199670 18458
## - income_cat 1 520 199751 18460
## - exercise 1 1151 200381 18476
## - menthlth_days 1 5929 205159 18594
## - gen_health 4 66459 265690 19880
##
## Step: AIC=18447.92
## physhlth_days ~ menthlth_days + sleep_hrs + age + sex + exercise +
## gen_health + income_cat
##
## Df Sum of Sq RSS AIC
## - sex 1 42 199305 18447
## <none> 199262 18448
## - sleep_hrs 1 334 199596 18454
## - age 1 427 199690 18457
## - income_cat 1 514 199776 18459
## - exercise 1 1222 200484 18476
## - menthlth_days 1 5921 205184 18592
## - gen_health 4 67347 266609 19896
##
## Step: AIC=18446.98
## physhlth_days ~ menthlth_days + sleep_hrs + age + exercise +
## gen_health + income_cat
##
## Df Sum of Sq RSS AIC
## <none> 199305 18447
## - sleep_hrs 1 337 199641 18453
## - age 1 409 199713 18455
## - income_cat 1 492 199797 18457
## - exercise 1 1214 200518 18475
## - menthlth_days 1 5882 205186 18590
## - gen_health 4 67980 267285 19906tidy(mod_backward, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Backward Elimination Result (AIC-based)",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 3.1864 | 0.6663 | 4.7819 | 0.0000 | 1.8800 | 4.4927 |
| menthlth_days | 0.1461 | 0.0120 | 12.1352 | 0.0000 | 0.1225 | 0.1697 |
| sleep_hrs | -0.1951 | 0.0672 | -2.9038 | 0.0037 | -0.3269 | -0.0634 |
| age | 0.0174 | 0.0054 | 3.1981 | 0.0014 | 0.0067 | 0.0281 |
| exerciseYes | -1.2877 | 0.2336 | -5.5127 | 0.0000 | -1.7457 | -0.8298 |
| gen_healthVery good | 0.4617 | 0.2441 | 1.8914 | 0.0586 | -0.0169 | 0.9403 |
| gen_healthGood | 1.6368 | 0.2600 | 6.2953 | 0.0000 | 1.1271 | 2.1465 |
| gen_healthFair | 7.0787 | 0.3616 | 19.5735 | 0.0000 | 6.3697 | 7.7876 |
| gen_healthPoor | 20.5084 | 0.5423 | 37.8149 | 0.0000 | 19.4452 | 21.5716 |
| income_cat | -0.1657 | 0.0472 | -3.5115 | 0.0004 | -0.2582 | -0.0732 |
Interpretation: AIC-based backward elimination removed sex, education, and BMI from the maximum model, arriving at a 9-parameter model (counting dummy variables). These are the same three variables that were non-significant in the maximum model. The retained predictors (mental health days, sleep, age, exercise, general health, and income) all have p-values below 0.05. The resulting model has Adjusted R² = 0.385, essentially identical to the maximum model (0.384), confirming that the dropped variables contributed negligible explanatory power.
3.3 Forward Selection
Forward selection starts with the intercept-only model and adds variables one at a time:
- Start with no predictors (intercept only)
- For each candidate variable, compute the partial F-statistic for adding it to the current model
- Add the variable with the smallest p-value (largest partial F)
- If its p-value is \(\leq \alpha\), keep it and repeat steps 2-3
- Stop when no remaining variable has a p-value \(\leq \alpha\)
# Automated forward selection using AIC
mod_null <- lm(physhlth_days ~ 1, data = brfss_ms)
mod_forward <- step(mod_null,
scope = list(lower = mod_null, upper = mod_max),
direction = "forward", trace = 1)## Start: AIC=20865.24
## physhlth_days ~ 1
##
## Df Sum of Sq RSS AIC
## + gen_health 4 115918 208518 18663
## + menthlth_days 1 29743 294693 20386
## + exercise 1 19397 305038 20559
## + income_cat 1 19104 305332 20564
## + education 3 5906 318530 20779
## + age 1 4173 320263 20802
## + bmi 1 4041 320395 20805
## + sleep_hrs 1 3717 320719 20810
## <none> 324435 20865
## + sex 1 7 324429 20867
##
## Step: AIC=18662.93
## physhlth_days ~ gen_health
##
## Df Sum of Sq RSS AIC
## + menthlth_days 1 6394.9 202123 18509
## + exercise 1 1652.4 206865 18625
## + income_cat 1 1306.9 207211 18634
## + sleep_hrs 1 756.1 207762 18647
## + bmi 1 91.2 208427 18663
## <none> 208518 18663
## + sex 1 38.5 208479 18664
## + age 1 32.2 208486 18664
## + education 3 145.0 208373 18666
##
## Step: AIC=18509.19
## physhlth_days ~ gen_health + menthlth_days
##
## Df Sum of Sq RSS AIC
## + exercise 1 1650.52 200472 18470
## + income_cat 1 817.89 201305 18491
## + age 1 464.73 201658 18500
## + sleep_hrs 1 257.79 201865 18505
## + bmi 1 90.51 202032 18509
## <none> 202123 18509
## + sex 1 3.00 202120 18511
## + education 3 111.58 202011 18512
##
## Step: AIC=18470.19
## physhlth_days ~ gen_health + menthlth_days + exercise
##
## Df Sum of Sq RSS AIC
## + income_cat 1 509.09 199963 18460
## + age 1 333.74 200139 18464
## + sleep_hrs 1 253.06 200219 18466
## <none> 200472 18470
## + bmi 1 21.21 200451 18472
## + sex 1 10.74 200462 18472
## + education 3 26.94 200445 18476
##
## Step: AIC=18459.48
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat
##
## Df Sum of Sq RSS AIC
## + age 1 321.97 199641 18453
## + sleep_hrs 1 250.25 199713 18455
## <none> 199963 18460
## + bmi 1 27.98 199935 18461
## + sex 1 27.17 199936 18461
## + education 3 26.66 199937 18465
##
## Step: AIC=18453.42
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age
##
## Df Sum of Sq RSS AIC
## + sleep_hrs 1 336.79 199305 18447
## <none> 199641 18453
## + sex 1 45.31 199596 18454
## + bmi 1 42.00 199599 18454
## + education 3 22.62 199619 18459
##
## Step: AIC=18446.98
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age + sleep_hrs
##
## Df Sum of Sq RSS AIC
## <none> 199305 18447
## + sex 1 42.328 199262 18448
## + bmi 1 34.434 199270 18448
## + education 3 24.800 199280 18452tidy(mod_forward, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Forward Selection Result (AIC-based)",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 3.1864 | 0.6663 | 4.7819 | 0.0000 | 1.8800 | 4.4927 |
| gen_healthVery good | 0.4617 | 0.2441 | 1.8914 | 0.0586 | -0.0169 | 0.9403 |
| gen_healthGood | 1.6368 | 0.2600 | 6.2953 | 0.0000 | 1.1271 | 2.1465 |
| gen_healthFair | 7.0787 | 0.3616 | 19.5735 | 0.0000 | 6.3697 | 7.7876 |
| gen_healthPoor | 20.5084 | 0.5423 | 37.8149 | 0.0000 | 19.4452 | 21.5716 |
| menthlth_days | 0.1461 | 0.0120 | 12.1352 | 0.0000 | 0.1225 | 0.1697 |
| exerciseYes | -1.2877 | 0.2336 | -5.5127 | 0.0000 | -1.7457 | -0.8298 |
| income_cat | -0.1657 | 0.0472 | -3.5115 | 0.0004 | -0.2582 | -0.0732 |
| age | 0.0174 | 0.0054 | 3.1981 | 0.0014 | 0.0067 | 0.0281 |
| sleep_hrs | -0.1951 | 0.0672 | -2.9038 | 0.0037 | -0.3269 | -0.0634 |
Interpretation: Forward selection arrived at the same final model as backward elimination, including the same 9 predictor terms. The order of entry is informative: general health entered first (the strongest predictor), followed by mental health days, exercise, income, age, and sleep. This ordering reflects each variable’s marginal contribution given the variables already in the model. The convergence of forward and backward methods on the same model increases our confidence in this particular subset, though this convergence is not guaranteed in general.
3.4 Stepwise Selection
Stepwise selection combines forward and backward: after adding a variable, it checks whether any previously entered variable should now be removed. This addresses a limitation of pure forward selection, where a variable that was useful early on may become redundant after other variables enter.
mod_stepwise <- step(mod_null,
scope = list(lower = mod_null, upper = mod_max),
direction = "both", trace = 1)## Start: AIC=20865.24
## physhlth_days ~ 1
##
## Df Sum of Sq RSS AIC
## + gen_health 4 115918 208518 18663
## + menthlth_days 1 29743 294693 20386
## + exercise 1 19397 305038 20559
## + income_cat 1 19104 305332 20564
## + education 3 5906 318530 20779
## + age 1 4173 320263 20802
## + bmi 1 4041 320395 20805
## + sleep_hrs 1 3717 320719 20810
## <none> 324435 20865
## + sex 1 7 324429 20867
##
## Step: AIC=18662.93
## physhlth_days ~ gen_health
##
## Df Sum of Sq RSS AIC
## + menthlth_days 1 6395 202123 18509
## + exercise 1 1652 206865 18625
## + income_cat 1 1307 207211 18634
## + sleep_hrs 1 756 207762 18647
## + bmi 1 91 208427 18663
## <none> 208518 18663
## + sex 1 38 208479 18664
## + age 1 32 208486 18664
## + education 3 145 208373 18666
## - gen_health 4 115918 324435 20865
##
## Step: AIC=18509.19
## physhlth_days ~ gen_health + menthlth_days
##
## Df Sum of Sq RSS AIC
## + exercise 1 1651 200472 18470
## + income_cat 1 818 201305 18491
## + age 1 465 201658 18500
## + sleep_hrs 1 258 201865 18505
## + bmi 1 91 202032 18509
## <none> 202123 18509
## + sex 1 3 202120 18511
## + education 3 112 202011 18512
## - menthlth_days 1 6395 208518 18663
## - gen_health 4 92570 294693 20386
##
## Step: AIC=18470.19
## physhlth_days ~ gen_health + menthlth_days + exercise
##
## Df Sum of Sq RSS AIC
## + income_cat 1 509 199963 18460
## + age 1 334 200139 18464
## + sleep_hrs 1 253 200219 18466
## <none> 200472 18470
## + bmi 1 21 200451 18472
## + sex 1 11 200462 18472
## + education 3 27 200445 18476
## - exercise 1 1651 202123 18509
## - menthlth_days 1 6393 206865 18625
## - gen_health 4 78857 279330 20121
##
## Step: AIC=18459.48
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat
##
## Df Sum of Sq RSS AIC
## + age 1 322 199641 18453
## + sleep_hrs 1 250 199713 18455
## <none> 199963 18460
## + bmi 1 28 199935 18461
## + sex 1 27 199936 18461
## + education 3 27 199937 18465
## - income_cat 1 509 200472 18470
## - exercise 1 1342 201305 18491
## - menthlth_days 1 5988 205952 18605
## - gen_health 4 72713 272676 20002
##
## Step: AIC=18453.42
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age
##
## Df Sum of Sq RSS AIC
## + sleep_hrs 1 337 199305 18447
## <none> 199641 18453
## + sex 1 45 199596 18454
## + bmi 1 42 199599 18454
## + education 3 23 199619 18459
## - age 1 322 199963 18460
## - income_cat 1 497 200139 18464
## - exercise 1 1231 200873 18482
## - menthlth_days 1 6304 205945 18607
## - gen_health 4 68936 268577 19928
##
## Step: AIC=18446.98
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age + sleep_hrs
##
## Df Sum of Sq RSS AIC
## <none> 199305 18447
## + sex 1 42 199262 18448
## + bmi 1 34 199270 18448
## + education 3 25 199280 18452
## - sleep_hrs 1 337 199641 18453
## - age 1 409 199713 18455
## - income_cat 1 492 199797 18457
## - exercise 1 1214 200518 18475
## - menthlth_days 1 5882 205186 18590
## - gen_health 4 67980 267285 19906tidy(mod_stepwise, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Stepwise Selection Result (AIC-based)",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 3.1864 | 0.6663 | 4.7819 | 0.0000 | 1.8800 | 4.4927 |
| gen_healthVery good | 0.4617 | 0.2441 | 1.8914 | 0.0586 | -0.0169 | 0.9403 |
| gen_healthGood | 1.6368 | 0.2600 | 6.2953 | 0.0000 | 1.1271 | 2.1465 |
| gen_healthFair | 7.0787 | 0.3616 | 19.5735 | 0.0000 | 6.3697 | 7.7876 |
| gen_healthPoor | 20.5084 | 0.5423 | 37.8149 | 0.0000 | 19.4452 | 21.5716 |
| menthlth_days | 0.1461 | 0.0120 | 12.1352 | 0.0000 | 0.1225 | 0.1697 |
| exerciseYes | -1.2877 | 0.2336 | -5.5127 | 0.0000 | -1.7457 | -0.8298 |
| income_cat | -0.1657 | 0.0472 | -3.5115 | 0.0004 | -0.2582 | -0.0732 |
| age | 0.0174 | 0.0054 | 3.1981 | 0.0014 | 0.0067 | 0.0281 |
| sleep_hrs | -0.1951 | 0.0672 | -2.9038 | 0.0037 | -0.3269 | -0.0634 |
Interpretation: The stepwise procedure, which allows both addition and removal at each step, also converges on the identical model. In this dataset, no variable that was added early became redundant after later variables entered, so no removals were needed. This three-way agreement (backward = forward = stepwise) is reassuring but should not be taken as proof that this is the “correct” model. All three methods optimize the same criterion (AIC) on the same data.
3.5 Comparing All Selection Methods
method_comparison <- tribble(
~Method, ~`Variables selected`, ~`Adj. R²`, ~AIC, ~BIC,
"Maximum model",
length(coef(mod_max)) - 1,
round(glance(mod_max)$adj.r.squared, 4),
round(AIC(mod_max), 1),
round(BIC(mod_max), 1),
"Backward (AIC)",
length(coef(mod_backward)) - 1,
round(glance(mod_backward)$adj.r.squared, 4),
round(AIC(mod_backward), 1),
round(BIC(mod_backward), 1),
"Forward (AIC)",
length(coef(mod_forward)) - 1,
round(glance(mod_forward)$adj.r.squared, 4),
round(AIC(mod_forward), 1),
round(BIC(mod_forward), 1),
"Stepwise (AIC)",
length(coef(mod_stepwise)) - 1,
round(glance(mod_stepwise)$adj.r.squared, 4),
round(AIC(mod_stepwise), 1),
round(BIC(mod_stepwise), 1)
)
method_comparison |>
kable(caption = "Comparison of Variable Selection Methods") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Method | Variables selected | Adj. R² | AIC | BIC |
|---|---|---|---|---|
| Maximum model | 14 | 0.3843 | 32645.8 | 32750.1 |
| Backward (AIC) | 9 | 0.3846 | 32638.4 | 32710.1 |
| Forward (AIC) | 9 | 0.3846 | 32638.4 | 32710.1 |
| Stepwise (AIC) | 9 | 0.3846 | 32638.4 | 32710.1 |
Interpretation: All three automated methods selected the same model with 9 predictor terms (Adjusted R² = 0.385, AIC = 32,638.4, BIC = 32,710.1). This model has a lower AIC and BIC than the maximum model (AIC = 32,645.8, BIC = 32,750.1), confirming that removing sex, education, and BMI improved parsimony without sacrificing fit. The modest improvement in BIC (40 points lower) is more notable than the AIC improvement (7 points lower), consistent with BIC’s stronger preference for simpler models. In practice, the maximum model and the selected model would produce very similar predictions, but the selected model is preferred for its efficiency.
3.6 Cautions About Automated Selection
Use automated selection with extreme caution.
Automated methods (forward, backward, stepwise) have well-documented problems:
They ignore the research question. The algorithm selects variables based purely on statistical fit. If you are building an associative model and the exposure is not statistically significant, the algorithm will remove it, which defeats the purpose.
They inflate Type I error. The repeated testing involved in stepwise procedures inflates the probability of including spurious predictors.
They are path-dependent. Forward and backward selection can yield different final models because the order of variable entry/removal matters.
They ignore subject-matter knowledge. A variable may be a known confounder from the literature even if it is not statistically significant in your sample.
p-values and CIs from the final model are biased. Because the model was selected to optimize fit, the reported p-values are anti-conservative (too small).
Recommendation: Use automated selection as an exploratory tool to generate candidate models, but make final decisions based on substantive knowledge, confounding assessment, and parsimony.
4. Model Selection for Associative Models
4.1 A Different Philosophy
In associative modeling, the exposure variable is always in the model. It is never a candidate for removal, regardless of its p-value. The question is which covariates to include alongside it.
The standard epidemiological approach to covariate selection:
- Identify the exposure(s) of interest (these stay in the model)
- Identify candidate confounders from the literature and bivariate analyses
- Use the 10% change-in-estimate rule to determine which confounders to retain
4.2 The 10% Change-in-Estimate Procedure
Recall from the Confounding lecture: a covariate is a confounder if removing it changes the exposure coefficient by more than 10%.
The systematic procedure:
- Fit the maximum model (exposure + all candidate confounders) and note the exposure \(\hat{\beta}\)
- Compute the 10% interval: \(\hat{\beta} \pm 0.10 \times |\hat{\beta}|\)
- Remove one candidate confounder at a time
- If the exposure \(\hat{\beta}\) stays within the 10% interval, the removed variable is not a confounder (drop it)
- If the exposure \(\hat{\beta}\) moves outside the interval, the variable is a confounder (keep it)
- Repeat until all covariates have been evaluated
# Exposure: exercise; Outcome: physhlth_days
# Maximum associative model
mod_assoc_max <- lm(physhlth_days ~ exercise + menthlth_days + sleep_hrs + age +
sex + education + income_cat + bmi,
data = brfss_ms)
b_exposure_max <- coef(mod_assoc_max)["exerciseYes"]
interval_low <- b_exposure_max - 0.10 * abs(b_exposure_max)
interval_high <- b_exposure_max + 0.10 * abs(b_exposure_max)
cat("Exposure coefficient in maximum model:", round(b_exposure_max, 4), "\n")## Exposure coefficient in maximum model: -3.0688## 10% interval: ( -3.3757 , -2.7619 )# Systematically remove one covariate at a time
covariates_to_test <- c("menthlth_days", "sleep_hrs", "age", "sex",
"education", "income_cat", "bmi")
assoc_table <- map_dfr(covariates_to_test, \(cov) {
# Build formula without this covariate
remaining <- setdiff(covariates_to_test, cov)
form <- as.formula(paste("physhlth_days ~ exercise +", paste(remaining, collapse = " + ")))
mod_reduced <- lm(form, data = brfss_ms)
b_reduced <- coef(mod_reduced)["exerciseYes"]
pct_change <- (b_reduced - b_exposure_max) / abs(b_exposure_max) * 100
tibble(
`Removed covariate` = cov,
`Exercise β (max)` = round(b_exposure_max, 4),
`Exercise β (without)` = round(b_reduced, 4),
`% Change` = round(pct_change, 1),
`Within 10%?` = ifelse(abs(pct_change) <= 10, "Yes (drop)", "No (keep)"),
Confounder = ifelse(abs(pct_change) > 10, "Yes", "No")
)
})
assoc_table |>
kable(caption = "Associative Model: Systematic Confounder Assessment for Exercise") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) |>
column_spec(6, bold = TRUE)| Removed covariate | Exercise β (max) | Exercise β (without) | % Change | Within 10%? | Confounder |
|---|---|---|---|---|---|
| menthlth_days | -3.0688 | -3.3725 | -9.9 | Yes (drop) | No |
| sleep_hrs | -3.0688 | -3.0950 | -0.9 | Yes (drop) | No |
| age | -3.0688 | -3.4150 | -11.3 | No (keep) | Yes |
| sex | -3.0688 | -3.0534 | 0.5 | Yes (drop) | No |
| education | -3.0688 | -3.1036 | -1.1 | Yes (drop) | No |
| income_cat | -3.0688 | -3.4544 | -12.6 | No (keep) | Yes |
| bmi | -3.0688 | -3.2411 | -5.6 | Yes (drop) | No |
Interpretation: The exercise coefficient in the maximum associative model is -3.07, meaning exercisers report about 3 fewer physically unhealthy days after adjusting for all covariates. The systematic assessment identifies two confounders: age (11.3% change when removed) and income (12.6% change when removed). Removing age strengthens the exercise effect (to -3.42), suggesting that age positively confounds the association (older adults exercise less and have more unhealthy days, so ignoring age makes exercise look less protective). Removing income also strengthens the effect (to -3.45), with a similar confounding mechanism (higher income is associated with both more exercise and fewer unhealthy days). The remaining covariates (mental health days, sleep, sex, education, BMI) all produce changes within the 10% interval, so they are not confounders and could be dropped from the associative model. The final associative model would include exercise, age, and income.
4.3 Associative Models with Interactions
If a statistically significant interaction is present (from the previous lecture), the approach changes:
- Stratify by the effect modifier
- Within each stratum, apply the 10% change-in-estimate procedure separately
- Report stratum-specific results
For example, if age \(\times\) exercise is significant:
- For exercisers: fit
physhlth_days ~ age + [confounders]and assess confounding - For non-exercisers: fit the same and assess confounding
- The set of confounders may differ between strata
4.4 Predictive vs. Associative: Side-by-Side Comparison
| Feature | Predictive | Associative |
|---|---|---|
| Exposure variable | No fixed exposure | Always in the model |
| Covariate selection | Based on statistical fit | Based on confounding assessment |
| Automated methods | Useful (with caution) | Generally inappropriate |
| 10% change rule | Not used | Primary tool |
| Interaction terms | Include if improves prediction | Include if effect modification is present |
| Primary criterion | Adj. R², AIC, BIC | Validity of exposure β |
| Parsimony | Fewer variables = less overfitting | Fewer variables = more efficient, if not confounders |
5. Cross-Validation: Evaluating Model Reliability
5.1 Why Cross-Validate?
A model that fits the training data well may perform poorly on new data (overfitting). Cross-validation estimates how well the model would perform on data it has not seen.
The simplest approach is k-fold cross-validation:
- Randomly split the data into \(k\) equally sized folds
- For each fold, train the model on the other \(k - 1\) folds and predict on the held-out fold
- Average the prediction error across all \(k\) folds
# 10-fold cross-validation comparison
set.seed(1220)
n <- nrow(brfss_ms)
k_folds <- 10
fold_id <- sample(rep(1:k_folds, length.out = n))
# Compare a small model, medium model, and full model
cv_results <- map_dfr(1:k_folds, \(fold) {
train <- brfss_ms[fold_id != fold, ]
test <- brfss_ms[fold_id == fold, ]
# Small model
m_small <- lm(physhlth_days ~ menthlth_days + gen_health, data = train)
pred_small <- predict(m_small, newdata = test)
# Medium model
m_med <- lm(physhlth_days ~ menthlth_days + gen_health + exercise + age + sleep_hrs,
data = train)
pred_med <- predict(m_med, newdata = test)
# Full model
m_full <- lm(physhlth_days ~ menthlth_days + sleep_hrs + age + sex + education +
exercise + gen_health + income_cat + bmi, data = train)
pred_full <- predict(m_full, newdata = test)
tibble(
fold = fold,
RMSE_small = sqrt(mean((test$physhlth_days - pred_small)^2)),
RMSE_medium = sqrt(mean((test$physhlth_days - pred_med)^2)),
RMSE_full = sqrt(mean((test$physhlth_days - pred_full)^2))
)
})
cv_summary <- cv_results |>
summarise(
across(starts_with("RMSE"), \(x) round(mean(x), 3))
)
tribble(
~Model, ~Predictors, ~`CV RMSE`,
"Small", "menthlth_days + gen_health", cv_summary$RMSE_small,
"Medium", "+ exercise + age + sleep_hrs", cv_summary$RMSE_medium,
"Full", "All 9 predictors", cv_summary$RMSE_full
) |>
kable(caption = "10-Fold Cross-Validation: Out-of-Sample RMSE") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | Predictors | CV RMSE |
|---|---|---|
| Small | menthlth_days + gen_health | 6.362 |
| Medium |
|
6.334 |
| Full | All 9 predictors | 6.334 |
Interpretation: RMSE is the average prediction error in the units of the outcome (days). A lower CV RMSE indicates better out-of-sample prediction. If the full model has a similar CV RMSE to the medium model, the additional predictors are not improving prediction and may represent overfitting.
Summary of Key Concepts
| Concept | Key Point |
|---|---|
| Maximum model | Start with all candidate predictors from literature and research question |
| Overfitting vs. underfitting | Overfitting = more variance; underfitting = bias |
| Parsimony | Simplest model that captures the important relationships |
| \(R^2\) | Always increases with more variables; useless alone for comparison |
| Adjusted \(R^2\) | Penalizes complexity; maximize it |
| AIC | Balances fit and complexity; minimize it |
| BIC | Heavier penalty than AIC; favors simpler models; minimize it |
| Partial F-test | Compares reduced to maximum model |
| Best subsets | Exhaustive search; leaps::regsubsets() |
| Backward elimination | Start full, remove highest p-value;
step(direction = "backward") |
| Forward selection | Start empty, add lowest p-value;
step(direction = "forward") |
| Stepwise | Forward + backward at each step;
step(direction = "both") |
| Caution | Automated methods ignore research questions and inflate Type I error |
| Associative models | Exposure stays in model; use 10% change-in-estimate for covariates |
| Cross-validation | Estimates out-of-sample performance; protects against overfitting |
Part 2: In-Class Lab Activity
EPI 553 — Model Selection Lab Due: End of class, March 24, 2026
Instructions
In this lab, you will practice both predictive and associative model selection using the BRFSS 2020 dataset. Work through each task systematically. You may discuss concepts with classmates, but your written answers and R code must be your own.
Submission: Knit your .Rmd to HTML and upload to Brightspace by end of class.
Data for the Lab
Use the saved analytic dataset from today’s lecture.
| Variable | Description | Type |
|---|---|---|
physhlth_days |
Physically unhealthy days in past 30 | Continuous (0–30) |
menthlth_days |
Mentally unhealthy days in past 30 | Continuous (0–30) |
sleep_hrs |
Sleep hours per night | Continuous (1–14) |
age |
Age in years (capped at 80) | Continuous |
sex |
Sex (Male/Female) | Factor |
education |
Education level (4 categories) | Factor |
exercise |
Any physical activity (Yes/No) | Factor |
gen_health |
General health status (5 categories) | Factor |
income_cat |
Household income (1–8 ordinal) | Numeric |
bmi |
Body mass index | Continuous |
library(tidyverse)
library(broom)
library(knitr)
library(kableExtra)
library(car)
library(leaps)
library(MASS)
library(gridExtra)
# Load the pre-saved analytic dataset
brfss_ms <- readRDS("/Users/emmanuelsmac/Desktop/brfss_ms_2020.rds")Task 1: Maximum Model and Criteria Comparison (15 points)
1a. Fit the Maximum Model
Fit the maximum model predicting physhlth_days from all
9 candidate predictors and report \(R^2\), Adjusted \(R^2\), AIC, and BIC.
# Maximum model: all 9 candidate predictors
mod_max_lab <- lm(physhlth_days ~ menthlth_days + sleep_hrs + age + sex +
education + exercise + gen_health + income_cat + bmi,
data = brfss_ms)
# Extract and display fit criteria
glance(mod_max_lab) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(across(everything(), \(x) round(x, 4))) |>
rename(`R²` = r.squared, `Adj. R²` = adj.r.squared) |>
kable(caption = "Task 1a — Maximum Model Fit Criteria") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| R² | Adj. R² | AIC | BIC |
|---|---|---|---|
| 0.386 | 0.3843 | 32645.79 | 32750.06 |
Interpretation: The maximum model with all 9 candidate predictors yields an \(R^2\) of approximately 0.386, indicating that the model explains about 38.6% of the variance in physically unhealthy days. The Adjusted \(R^2\) is 0.384, which accounts for model complexity. The AIC and BIC values serve as baselines for comparison with simpler models below.
1b. Fit the Minimal Model
Fit a minimal model using only menthlth_days and
age. Compare the four criteria to the maximum model.
# Minimal model: only two predictors
mod_min_lab <- lm(physhlth_days ~ menthlth_days + age, data = brfss_ms)
# Side-by-side comparison
bind_rows(
glance(mod_max_lab) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(Model = "Maximum (9 predictors)"),
glance(mod_min_lab) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(Model = "Minimal (menthlth_days + age)")
) |>
dplyr::select(Model, r.squared, adj.r.squared, AIC, BIC) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
rename(`R²` = r.squared, `Adj. R²` = adj.r.squared) |>
kable(caption = "Task 1b — Maximum vs. Minimal Model Comparison") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | R² | Adj. R² | AIC | BIC |
|---|---|---|---|---|
| Maximum (9 predictors) | 0.386 | 0.3843 | 32645.79 | 32750.06 |
| Minimal (menthlth_days + age) | 0.115 | 0.1146 | 34449.78 | 34475.85 |
Interpretation: The minimal model (mental health days + age) has a substantially lower \(R^2\) (~0.09) and Adjusted \(R^2\) (~0.09) compared to the maximum model (~0.386 and ~0.384, respectively), indicating it explains far less variance. Both AIC and BIC are considerably higher for the minimal model, confirming it is a much worse fit. The large drops in all four criteria when moving to the maximum model suggest that the additional predictors — particularly general health status and exercise — carry substantial explanatory power.
1c. Why \(R^2\) Alone Is a Poor Comparison Criterion
Answer: \(R^2\) is a poor criterion for comparing these two models because it always increases (or stays the same) when additional predictors are added to a model, regardless of whether those predictors carry genuine explanatory value. In this case, the maximum model has a higher \(R^2\) than the minimal model partly because it contains 9 predictors rather than 2 — some of those predictors may contribute very little information but will still increase \(R^2\) trivially. This makes \(R^2\) unable to distinguish between a model that is genuinely better and one that is merely larger.
Adjusted \(R^2\), AIC, and BIC are better choices because all three penalize for model complexity:
- Adjusted \(R^2\) incorporates a degrees-of-freedom penalty so that adding an uninformative predictor can decrease it. This prevents the criterion from always rewarding larger models.
- AIC penalizes by \(2k\) (where \(k\) is the number of estimated parameters), balancing goodness of fit against complexity. It can decrease when a variable adds more noise than signal.
- BIC applies an even heavier penalty of \(k \log(n)\), which grows with sample size. In large samples like ours (\(n = 5,000\)), BIC strongly favors parsimony and is more likely to select the truly important predictors.
In short, raw \(R^2\) only tells you how well the model fits the training data; Adjusted \(R^2\), AIC, and BIC all provide an estimate of the trade-off between model fit and unnecessary complexity.
Task 2: Best Subsets Regression (20 points)
2a. Adjusted \(R^2\) vs. Number of Variables
Use leaps::regsubsets() with nvmax = 15.
Plot Adjusted \(R^2\) vs. number of
variables and identify where it plateaus.
# Best subsets regression
lab_subsets <- regsubsets(
physhlth_days ~ menthlth_days + sleep_hrs + age + sex + education +
exercise + gen_health + income_cat + bmi,
data = brfss_ms,
nvmax = 15,
method = "exhaustive"
)
lab_summary <- summary(lab_subsets)
# Build metrics tibble
lab_metrics <- tibble(
p = 1:length(lab_summary$adjr2),
adjr2 = lab_summary$adjr2,
bic = lab_summary$bic
)
# Plot Adjusted R²
ggplot(lab_metrics, aes(x = p, y = adjr2)) +
geom_line(linewidth = 1.1, color = "steelblue") +
geom_point(size = 3, color = "steelblue") +
geom_vline(xintercept = which.max(lab_summary$adjr2),
linetype = "dashed", color = "tomato", linewidth = 1) +
annotate("text",
x = which.max(lab_summary$adjr2) + 0.3,
y = min(lab_metrics$adjr2) + 0.01,
label = paste0("Max at p = ", which.max(lab_summary$adjr2)),
hjust = 0, color = "tomato", size = 4) +
labs(
title = "Best Subsets: Adjusted R² by Model Size",
subtitle = "Red dashed line = model with highest Adjusted R²",
x = "Number of Variables (predictor df)",
y = "Adjusted R²"
) +
theme_minimal(base_size = 13)
Task 2a — Best Subsets: Adjusted R² by Model Size
## Model size with highest Adjusted R²: 10 variablesInterpretation: Adjusted \(R^2\) rises steeply for the first several predictors — the largest single jump occurs when general health status enters — and then plateaus around 10–11 variables. Adding variables beyond that point produces negligible or no improvement in Adjusted \(R^2\), indicating those additional predictors do not contribute meaningful explanatory power after accounting for the degrees-of-freedom penalty. The plateau signals that the marginal gain from extra complexity is essentially zero.
2b. BIC vs. Number of Variables
Plot BIC vs. number of variables and identify which model size minimizes BIC.
ggplot(lab_metrics, aes(x = p, y = bic)) +
geom_line(linewidth = 1.1, color = "steelblue") +
geom_point(size = 3, color = "steelblue") +
geom_vline(xintercept = which.min(lab_summary$bic),
linetype = "dashed", color = "tomato", linewidth = 1) +
annotate("text",
x = which.min(lab_summary$bic) + 0.3,
y = max(lab_metrics$bic) - 10,
label = paste0("Min at p = ", which.min(lab_summary$bic)),
hjust = 0, color = "tomato", size = 4) +
labs(
title = "Best Subsets: BIC by Model Size",
subtitle = "Red dashed line = model minimizing BIC",
x = "Number of Variables (predictor df)",
y = "BIC"
) +
theme_minimal(base_size = 13)
Task 2b — Best Subsets: BIC by Model Size
## Model size minimizing BIC: 8 variablesInterpretation: BIC decreases sharply as the most important predictors enter the model, reaches its minimum at a smaller number of variables than Adjusted \(R^2\), then begins to increase again as additional predictors are added. This pattern reflects BIC’s heavier penalty for complexity (\(k \log n\) vs. \(2k\) for AIC). BIC therefore selects a more parsimonious model than Adjusted \(R^2\), preferring to drop variables whose marginal contribution is too small to justify the additional parameter.
2c. Identify BIC-Best Model and Fit with lm()
Identify the variables in the BIC-minimizing model, fit it explicitly, and report coefficients.
# Identify BIC-best model size and variables
best_bic_idx <- which.min(lab_summary$bic)
# regsubsets returns dummy-coded names (e.g. "exerciseYes", "gen_healthPoor").
# Map these back to the original variable names so lm() handles factors correctly.
all_predictors <- c("menthlth_days", "sleep_hrs", "age", "sex",
"education", "exercise", "gen_health", "income_cat", "bmi")
selected_cols <- names(which(lab_summary$which[best_bic_idx, -1]))
bic_vars <- all_predictors[sapply(all_predictors, function(v) {
any(grepl(paste0("^", v), selected_cols))
})]
cat("BIC-best model size:", best_bic_idx, "variables\n")## BIC-best model size: 8 variables## Original variables included:## menthlth_days
## sleep_hrs
## age
## exercise
## gen_health
## income_cat# Fit using original variable names (lm() expands factors automatically)
bic_formula <- as.formula(
paste("physhlth_days ~", paste(bic_vars, collapse = " + "))
)
# Fit the BIC-best model explicitly
mod_bic_best <- lm(bic_formula, data = brfss_ms)
tidy(mod_bic_best, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Task 2c — BIC-Best Model Coefficients",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 3.1864 | 0.6663 | 4.7819 | 0.0000 | 1.8800 | 4.4927 |
| menthlth_days | 0.1461 | 0.0120 | 12.1352 | 0.0000 | 0.1225 | 0.1697 |
| sleep_hrs | -0.1951 | 0.0672 | -2.9038 | 0.0037 | -0.3269 | -0.0634 |
| age | 0.0174 | 0.0054 | 3.1981 | 0.0014 | 0.0067 | 0.0281 |
| exerciseYes | -1.2877 | 0.2336 | -5.5127 | 0.0000 | -1.7457 | -0.8298 |
| gen_healthVery good | 0.4617 | 0.2441 | 1.8914 | 0.0586 | -0.0169 | 0.9403 |
| gen_healthGood | 1.6368 | 0.2600 | 6.2953 | 0.0000 | 1.1271 | 2.1465 |
| gen_healthFair | 7.0787 | 0.3616 | 19.5735 | 0.0000 | 6.3697 | 7.7876 |
| gen_healthPoor | 20.5084 | 0.5423 | 37.8149 | 0.0000 | 19.4452 | 21.5716 |
| income_cat | -0.1657 | 0.0472 | -3.5115 | 0.0004 | -0.2582 | -0.0732 |
glance(mod_bic_best) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(across(everything(), \(x) round(x, 4))) |>
rename(`R²` = r.squared, `Adj. R²` = adj.r.squared) |>
kable(caption = "Task 2c — BIC-Best Model Fit Statistics") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| R² | Adj. R² | AIC | BIC |
|---|---|---|---|
| 0.3857 | 0.3846 | 32638.37 | 32710.06 |
Interpretation: The BIC-best model retains the predictors most consistently identified as important: mental health days, sleep hours, age, exercise, general health status (multiple dummy levels), and income category. Variables such as sex, education levels, and BMI — which had the weakest marginal contributions in the maximum model — are excluded. All retained coefficients are statistically significant. General health status (particularly “Poor” and “Fair”) remains the dominant predictor, consistent with its behavior throughout the lecture analysis.
2d. BIC-Best vs. Adjusted \(R^2\)-Best: Comparison and Preference
# Identify Adj R²-best model
best_adjr2_idx <- which.max(lab_summary$adjr2)
# Map dummy names back to original variable names for the Adj R²-best model
selected_cols_adjr2 <- names(which(lab_summary$which[best_adjr2_idx, -1]))
adjr2_vars <- all_predictors[sapply(all_predictors, function(v) {
any(grepl(paste0("^", v), selected_cols_adjr2))
})]
cat("Adj R²-best model size:", best_adjr2_idx, "variables\n")## Adj R²-best model size: 10 variables## Variables:
## menthlth_days
## sleep_hrs
## age
## sex
## exercise
## gen_health
## income_cat## BIC-best model size: 8 variables## Variables:
## menthlth_days
## sleep_hrs
## age
## exercise
## gen_health
## income_cat# Fit Adj R²-best model using original variable names
adjr2_formula <- as.formula(
paste("physhlth_days ~", paste(adjr2_vars, collapse = " + "))
)
mod_adjr2_best <- lm(adjr2_formula, data = brfss_ms)
# Compare the two
bind_rows(
glance(mod_bic_best) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(Model = paste0("BIC-best (", best_bic_idx, " vars)")),
glance(mod_adjr2_best) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(Model = paste0("Adj R²-best (", best_adjr2_idx, " vars)"))
) |>
dplyr::select(Model, everything()) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
rename(`R²` = r.squared, `Adj. R²` = adj.r.squared) |>
kable(caption = "Task 2d — BIC-Best vs. Adjusted R²-Best Model Comparison") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | R² | Adj. R² | AIC | BIC |
|---|---|---|---|---|
| BIC-best (8 vars) | 0.3857 | 0.3846 | 32638.37 | 32710.06 |
| Adj R²-best (10 vars) | 0.3858 | 0.3846 | 32639.31 | 32717.51 |
Answer: The BIC-best and Adjusted \(R^2\)-best models are not identical. The Adjusted \(R^2\)-best model includes more variables because Adjusted \(R^2\) applies a lighter penalty for complexity than BIC. In large samples (\(n = 5,000\)), BIC’s penalty term (\(k \log n\)) is substantially larger than Adjusted \(R^2\)’s, so BIC aggressively trims predictors whose marginal \(R^2\) gain is small.
Which is preferable? For this predictive modeling context, I would prefer the BIC-best model. The Adjusted \(R^2\)-best and BIC-best models have nearly identical fit statistics — the difference in Adjusted \(R^2\) is negligible — yet the BIC-best model is more parsimonious. A simpler model is less likely to overfit, will generalize better to new data, and is more interpretable. When predictive performance is essentially equivalent, parsimony is the tiebreaker, and BIC is specifically designed to enforce that principle.
Task 3: Automated Selection Methods (20 points)
3a. Backward Elimination (AIC)
# Fit maximum model for lab (same as mod_max_lab above)
mod_null_lab <- lm(physhlth_days ~ 1, data = brfss_ms)
# AIC-based backward elimination
mod_backward_lab <- step(mod_max_lab, direction = "backward", trace = 1)## Start: AIC=18454.4
## physhlth_days ~ menthlth_days + sleep_hrs + age + sex + education +
## exercise + gen_health + income_cat + bmi
##
## Df Sum of Sq RSS AIC
## - education 3 29 199231 18449
## - bmi 1 32 199234 18453
## - sex 1 43 199245 18454
## <none> 199202 18454
## - sleep_hrs 1 329 199530 18461
## - age 1 434 199636 18463
## - income_cat 1 521 199722 18466
## - exercise 1 1174 200376 18482
## - menthlth_days 1 5898 205100 18598
## - gen_health 4 66437 265639 19886
##
## Step: AIC=18449.13
## physhlth_days ~ menthlth_days + sleep_hrs + age + sex + exercise +
## gen_health + income_cat + bmi
##
## Df Sum of Sq RSS AIC
## - bmi 1 32 199262 18448
## - sex 1 40 199270 18448
## <none> 199231 18449
## - sleep_hrs 1 327 199557 18455
## - age 1 439 199670 18458
## - income_cat 1 520 199751 18460
## - exercise 1 1151 200381 18476
## - menthlth_days 1 5929 205159 18594
## - gen_health 4 66459 265690 19880
##
## Step: AIC=18447.92
## physhlth_days ~ menthlth_days + sleep_hrs + age + sex + exercise +
## gen_health + income_cat
##
## Df Sum of Sq RSS AIC
## - sex 1 42 199305 18447
## <none> 199262 18448
## - sleep_hrs 1 334 199596 18454
## - age 1 427 199690 18457
## - income_cat 1 514 199776 18459
## - exercise 1 1222 200484 18476
## - menthlth_days 1 5921 205184 18592
## - gen_health 4 67347 266609 19896
##
## Step: AIC=18446.98
## physhlth_days ~ menthlth_days + sleep_hrs + age + exercise +
## gen_health + income_cat
##
## Df Sum of Sq RSS AIC
## <none> 199305 18447
## - sleep_hrs 1 337 199641 18453
## - age 1 409 199713 18455
## - income_cat 1 492 199797 18457
## - exercise 1 1214 200518 18475
## - menthlth_days 1 5882 205186 18590
## - gen_health 4 67980 267285 19906tidy(mod_backward_lab, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Task 3a — Backward Elimination Result (AIC)",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 3.1864 | 0.6663 | 4.7819 | 0.0000 | 1.8800 | 4.4927 |
| menthlth_days | 0.1461 | 0.0120 | 12.1352 | 0.0000 | 0.1225 | 0.1697 |
| sleep_hrs | -0.1951 | 0.0672 | -2.9038 | 0.0037 | -0.3269 | -0.0634 |
| age | 0.0174 | 0.0054 | 3.1981 | 0.0014 | 0.0067 | 0.0281 |
| exerciseYes | -1.2877 | 0.2336 | -5.5127 | 0.0000 | -1.7457 | -0.8298 |
| gen_healthVery good | 0.4617 | 0.2441 | 1.8914 | 0.0586 | -0.0169 | 0.9403 |
| gen_healthGood | 1.6368 | 0.2600 | 6.2953 | 0.0000 | 1.1271 | 2.1465 |
| gen_healthFair | 7.0787 | 0.3616 | 19.5735 | 0.0000 | 6.3697 | 7.7876 |
| gen_healthPoor | 20.5084 | 0.5423 | 37.8149 | 0.0000 | 19.4452 | 21.5716 |
| income_cat | -0.1657 | 0.0472 | -3.5115 | 0.0004 | -0.2582 | -0.0732 |
# Report what was removed
all_max_vars <- names(coef(mod_max_lab))[-1]
all_backward_vars <- names(coef(mod_backward_lab))[-1]
removed_vars <- setdiff(
attr(terms(mod_max_lab), "term.labels"),
attr(terms(mod_backward_lab), "term.labels")
)
cat("\nVariables REMOVED by backward elimination:\n",
paste(" ", removed_vars, collapse = "\n"), "\n")##
## Variables REMOVED by backward elimination:
## sex
## education
## bmicat("\nVariables RETAINED:\n",
paste(" ", attr(terms(mod_backward_lab), "term.labels"), collapse = "\n"), "\n")##
## Variables RETAINED:
## menthlth_days
## sleep_hrs
## age
## exercise
## gen_health
## income_catInterpretation: AIC-based backward elimination
removed the predictors with the weakest marginal contributions —
specifically sex, education, and
bmi — leaving a model with the remaining stronger
predictors. The eliminated variables had the highest p-values in the
maximum model and their removal did not meaningfully increase AIC,
confirming they add negligible information. The retained variables —
mental health days, sleep, age, exercise, general health, and income —
are all statistically significant and align with the best-subsets
results from Task 2.
3b. Forward Selection (AIC)
# AIC-based forward selection
mod_forward_lab <- step(mod_null_lab,
scope = list(lower = mod_null_lab, upper = mod_max_lab),
direction = "forward",
trace = 1)## Start: AIC=20865.24
## physhlth_days ~ 1
##
## Df Sum of Sq RSS AIC
## + gen_health 4 115918 208518 18663
## + menthlth_days 1 29743 294693 20386
## + exercise 1 19397 305038 20559
## + income_cat 1 19104 305332 20564
## + education 3 5906 318530 20779
## + age 1 4173 320263 20802
## + bmi 1 4041 320395 20805
## + sleep_hrs 1 3717 320719 20810
## <none> 324435 20865
## + sex 1 7 324429 20867
##
## Step: AIC=18662.93
## physhlth_days ~ gen_health
##
## Df Sum of Sq RSS AIC
## + menthlth_days 1 6394.9 202123 18509
## + exercise 1 1652.4 206865 18625
## + income_cat 1 1306.9 207211 18634
## + sleep_hrs 1 756.1 207762 18647
## + bmi 1 91.2 208427 18663
## <none> 208518 18663
## + sex 1 38.5 208479 18664
## + age 1 32.2 208486 18664
## + education 3 145.0 208373 18666
##
## Step: AIC=18509.19
## physhlth_days ~ gen_health + menthlth_days
##
## Df Sum of Sq RSS AIC
## + exercise 1 1650.52 200472 18470
## + income_cat 1 817.89 201305 18491
## + age 1 464.73 201658 18500
## + sleep_hrs 1 257.79 201865 18505
## + bmi 1 90.51 202032 18509
## <none> 202123 18509
## + sex 1 3.00 202120 18511
## + education 3 111.58 202011 18512
##
## Step: AIC=18470.19
## physhlth_days ~ gen_health + menthlth_days + exercise
##
## Df Sum of Sq RSS AIC
## + income_cat 1 509.09 199963 18460
## + age 1 333.74 200139 18464
## + sleep_hrs 1 253.06 200219 18466
## <none> 200472 18470
## + bmi 1 21.21 200451 18472
## + sex 1 10.74 200462 18472
## + education 3 26.94 200445 18476
##
## Step: AIC=18459.48
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat
##
## Df Sum of Sq RSS AIC
## + age 1 321.97 199641 18453
## + sleep_hrs 1 250.25 199713 18455
## <none> 199963 18460
## + bmi 1 27.98 199935 18461
## + sex 1 27.17 199936 18461
## + education 3 26.66 199937 18465
##
## Step: AIC=18453.42
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age
##
## Df Sum of Sq RSS AIC
## + sleep_hrs 1 336.79 199305 18447
## <none> 199641 18453
## + sex 1 45.31 199596 18454
## + bmi 1 42.00 199599 18454
## + education 3 22.62 199619 18459
##
## Step: AIC=18446.98
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age + sleep_hrs
##
## Df Sum of Sq RSS AIC
## <none> 199305 18447
## + sex 1 42.328 199262 18448
## + bmi 1 34.434 199270 18448
## + education 3 24.800 199280 18452tidy(mod_forward_lab, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Task 3b — Forward Selection Result (AIC)",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 3.1864 | 0.6663 | 4.7819 | 0.0000 | 1.8800 | 4.4927 |
| gen_healthVery good | 0.4617 | 0.2441 | 1.8914 | 0.0586 | -0.0169 | 0.9403 |
| gen_healthGood | 1.6368 | 0.2600 | 6.2953 | 0.0000 | 1.1271 | 2.1465 |
| gen_healthFair | 7.0787 | 0.3616 | 19.5735 | 0.0000 | 6.3697 | 7.7876 |
| gen_healthPoor | 20.5084 | 0.5423 | 37.8149 | 0.0000 | 19.4452 | 21.5716 |
| menthlth_days | 0.1461 | 0.0120 | 12.1352 | 0.0000 | 0.1225 | 0.1697 |
| exerciseYes | -1.2877 | 0.2336 | -5.5127 | 0.0000 | -1.7457 | -0.8298 |
| income_cat | -0.1657 | 0.0472 | -3.5115 | 0.0004 | -0.2582 | -0.0732 |
| age | 0.0174 | 0.0054 | 3.1981 | 0.0014 | 0.0067 | 0.0281 |
| sleep_hrs | -0.1951 | 0.0672 | -2.9038 | 0.0037 | -0.3269 | -0.0634 |
# Check if models are identical
backward_terms <- sort(attr(terms(mod_backward_lab), "term.labels"))
forward_terms <- sort(attr(terms(mod_forward_lab), "term.labels"))
cat("\nBackward and forward models identical?",
identical(backward_terms, forward_terms), "\n")##
## Backward and forward models identical? TRUEInterpretation: Forward selection begins with no predictors and adds the variable providing the greatest AIC reduction at each step. General health enters first (as expected, given its dominant explanatory power), followed by mental health days, exercise, income, age, and sleep. Forward selection arrives at the same final model as backward elimination — the same set of predictors is retained in both approaches. This convergence is reassuring: when two search strategies operating from opposite ends reach the same answer, it increases confidence that the selected model represents a stable, locally optimal subset.
3c. Stepwise and Summary Comparison Table
# Stepwise (both directions)
mod_stepwise_lab <- step(mod_null_lab,
scope = list(lower = mod_null_lab, upper = mod_max_lab),
direction = "both",
trace = 1)## Start: AIC=20865.24
## physhlth_days ~ 1
##
## Df Sum of Sq RSS AIC
## + gen_health 4 115918 208518 18663
## + menthlth_days 1 29743 294693 20386
## + exercise 1 19397 305038 20559
## + income_cat 1 19104 305332 20564
## + education 3 5906 318530 20779
## + age 1 4173 320263 20802
## + bmi 1 4041 320395 20805
## + sleep_hrs 1 3717 320719 20810
## <none> 324435 20865
## + sex 1 7 324429 20867
##
## Step: AIC=18662.93
## physhlth_days ~ gen_health
##
## Df Sum of Sq RSS AIC
## + menthlth_days 1 6395 202123 18509
## + exercise 1 1652 206865 18625
## + income_cat 1 1307 207211 18634
## + sleep_hrs 1 756 207762 18647
## + bmi 1 91 208427 18663
## <none> 208518 18663
## + sex 1 38 208479 18664
## + age 1 32 208486 18664
## + education 3 145 208373 18666
## - gen_health 4 115918 324435 20865
##
## Step: AIC=18509.19
## physhlth_days ~ gen_health + menthlth_days
##
## Df Sum of Sq RSS AIC
## + exercise 1 1651 200472 18470
## + income_cat 1 818 201305 18491
## + age 1 465 201658 18500
## + sleep_hrs 1 258 201865 18505
## + bmi 1 91 202032 18509
## <none> 202123 18509
## + sex 1 3 202120 18511
## + education 3 112 202011 18512
## - menthlth_days 1 6395 208518 18663
## - gen_health 4 92570 294693 20386
##
## Step: AIC=18470.19
## physhlth_days ~ gen_health + menthlth_days + exercise
##
## Df Sum of Sq RSS AIC
## + income_cat 1 509 199963 18460
## + age 1 334 200139 18464
## + sleep_hrs 1 253 200219 18466
## <none> 200472 18470
## + bmi 1 21 200451 18472
## + sex 1 11 200462 18472
## + education 3 27 200445 18476
## - exercise 1 1651 202123 18509
## - menthlth_days 1 6393 206865 18625
## - gen_health 4 78857 279330 20121
##
## Step: AIC=18459.48
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat
##
## Df Sum of Sq RSS AIC
## + age 1 322 199641 18453
## + sleep_hrs 1 250 199713 18455
## <none> 199963 18460
## + bmi 1 28 199935 18461
## + sex 1 27 199936 18461
## + education 3 27 199937 18465
## - income_cat 1 509 200472 18470
## - exercise 1 1342 201305 18491
## - menthlth_days 1 5988 205952 18605
## - gen_health 4 72713 272676 20002
##
## Step: AIC=18453.42
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age
##
## Df Sum of Sq RSS AIC
## + sleep_hrs 1 337 199305 18447
## <none> 199641 18453
## + sex 1 45 199596 18454
## + bmi 1 42 199599 18454
## + education 3 23 199619 18459
## - age 1 322 199963 18460
## - income_cat 1 497 200139 18464
## - exercise 1 1231 200873 18482
## - menthlth_days 1 6304 205945 18607
## - gen_health 4 68936 268577 19928
##
## Step: AIC=18446.98
## physhlth_days ~ gen_health + menthlth_days + exercise + income_cat +
## age + sleep_hrs
##
## Df Sum of Sq RSS AIC
## <none> 199305 18447
## + sex 1 42 199262 18448
## + bmi 1 34 199270 18448
## + education 3 25 199280 18452
## - sleep_hrs 1 337 199641 18453
## - age 1 409 199713 18455
## - income_cat 1 492 199797 18457
## - exercise 1 1214 200518 18475
## - menthlth_days 1 5882 205186 18590
## - gen_health 4 67980 267285 19906tribble(
~Method, ~`# Variables`, ~`Adj. R²`, ~AIC, ~BIC,
"Maximum model", length(coef(mod_max_lab)) - 1, round(glance(mod_max_lab)$adj.r.squared, 4), round(AIC(mod_max_lab), 1), round(BIC(mod_max_lab), 1),
"Backward (AIC)", length(coef(mod_backward_lab)) - 1, round(glance(mod_backward_lab)$adj.r.squared, 4), round(AIC(mod_backward_lab), 1), round(BIC(mod_backward_lab), 1),
"Forward (AIC)", length(coef(mod_forward_lab)) - 1, round(glance(mod_forward_lab)$adj.r.squared, 4), round(AIC(mod_forward_lab), 1), round(BIC(mod_forward_lab), 1),
"Stepwise (AIC)", length(coef(mod_stepwise_lab)) - 1, round(glance(mod_stepwise_lab)$adj.r.squared, 4), round(AIC(mod_stepwise_lab), 1), round(BIC(mod_stepwise_lab), 1)
) |>
kable(caption = "Task 3c — Comparison of Automated Selection Methods") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Method | # Variables | Adj. R² | AIC | BIC |
|---|---|---|---|---|
| Maximum model | 14 | 0.3843 | 32645.8 | 32750.1 |
| Backward (AIC) | 9 | 0.3846 | 32638.4 | 32710.1 |
| Forward (AIC) | 9 | 0.3846 | 32638.4 | 32710.1 |
| Stepwise (AIC) | 9 | 0.3846 | 32638.4 | 32710.1 |
Interpretation: All three automated methods — backward, forward, and stepwise — converge on the same model, with the same Adjusted \(R^2\), AIC, and BIC. Compared to the maximum model, the selected model has a lower AIC and BIC (reflecting better parsimony) while maintaining essentially the same Adjusted \(R^2\). The uniform convergence of all three methods on an identical result, combined with agreement with the best-subsets analysis in Task 2, strongly suggests this is a stable and well-supported variable subset for predictive purposes.
3d. Three Reasons Not to Blindly Trust Automated Selection
Answer:
Automated methods ignore the research question and subject-matter knowledge. The algorithm makes decisions based solely on statistical fit — it has no concept of causal pathways, known confounders from the literature, or biologically meaningful relationships. A variable that is a well-established confounder in the epidemiological literature may be dropped simply because it is not statistically significant in a particular sample, which would bias the results in an associative analysis.
Repeated testing inflates Type I error. Each step of a backward, forward, or stepwise procedure involves one or more significance tests. When many tests are conducted on the same data, the probability that at least one spurious predictor enters (or a true predictor is incorrectly removed) is substantially higher than the nominal \(\alpha\) level. This means the final p-values and confidence intervals from the selected model are anti-conservative and cannot be taken at face value.
Path-dependence and instability. Forward and backward methods can arrive at different final models because the order in which variables are added or removed affects subsequent decisions. The selected model is sensitive to small perturbations in the data (e.g., a slightly different random sample could lead to a different model), which undermines reproducibility and interpretability.
Most relevant concern for epidemiological research: The first concern — ignoring subject-matter knowledge and the research question — is the most critical in epidemiology. In associative (etiologic) research, the goal is to estimate the unconfounded effect of an exposure. Automated selection may drop important confounders (producing biased estimates) or keep irrelevant variables that introduce multicollinearity, neither of which serves the scientific goal. Epidemiological model building must be grounded in biological plausibility and causal reasoning, not statistical optimization alone.
Task 4: Associative Model Building (25 points)
Exposure: sleep_hrs |
Outcome: physhlth_days
4a. Crude Model
Fit the crude (unadjusted) model and report the sleep coefficient.
# Crude model
mod_crude_sleep <- lm(physhlth_days ~ sleep_hrs, data = brfss_ms)
tidy(mod_crude_sleep, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Task 4a — Crude Model: physhlth_days ~ sleep_hrs",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 7.9110 | 0.5959 | 13.2755 | 0 | 6.7428 | 9.0793 |
| sleep_hrs | -0.6321 | 0.0831 | -7.6104 | 0 | -0.7949 | -0.4693 |
b_crude_sleep <- coef(mod_crude_sleep)["sleep_hrs"]
cat("\nCrude sleep coefficient:", round(b_crude_sleep, 4), "\n")##
## Crude sleep coefficient: -0.6321cat("Interpretation: Each additional hour of sleep is associated with",
round(b_crude_sleep, 4), "physically unhealthy days (unadjusted).\n")## Interpretation: Each additional hour of sleep is associated with -0.6321 physically unhealthy days (unadjusted).Interpretation: The crude (unadjusted) sleep coefficient is negative, indicating that each additional hour of sleep per night is associated with fewer physically unhealthy days in the past 30 days. This bivariate association is statistically significant. However, the crude estimate may be confounded by variables such as age, general health, mental health, and socioeconomic status — people who sleep more or less may systematically differ on these important characteristics. We need to assess confounding before drawing any causal conclusions.
4b. Maximum Associative Model and 10% Confounder Assessment
Fit the full adjusted model, compute the 10% interval around the sleep coefficient, then remove each covariate one at a time to identify confounders.
# Maximum associative model: sleep + all other covariates (gen_health included)
mod_assoc_max_sleep <- lm(physhlth_days ~ sleep_hrs + menthlth_days + age + sex +
education + exercise + gen_health + income_cat + bmi,
data = brfss_ms)
b_sleep_max <- coef(mod_assoc_max_sleep)["sleep_hrs"]
interval_low <- b_sleep_max - 0.10 * abs(b_sleep_max)
interval_high <- b_sleep_max + 0.10 * abs(b_sleep_max)
cat("Sleep coefficient in maximum associative model:", round(b_sleep_max, 4), "\n")## Sleep coefficient in maximum associative model: -0.193## 10% interval: ( -0.2123 , -0.1737 )# Systematic one-at-a-time removal
covariates_sleep <- c("menthlth_days", "age", "sex", "education",
"exercise", "gen_health", "income_cat", "bmi")
assoc_sleep_table <- map_dfr(covariates_sleep, \(cov) {
remaining <- setdiff(covariates_sleep, cov)
form <- as.formula(
paste("physhlth_days ~ sleep_hrs +", paste(remaining, collapse = " + "))
)
mod_r <- lm(form, data = brfss_ms)
b_r <- coef(mod_r)["sleep_hrs"]
pct_change <- (b_r - b_sleep_max) / abs(b_sleep_max) * 100
tibble(
`Removed Covariate` = cov,
`Sleep β (max model)` = round(b_sleep_max, 4),
`Sleep β (without)` = round(b_r, 4),
`% Change` = round(pct_change, 1),
`Within 10%?` = ifelse(abs(pct_change) <= 10, "Yes — drop", "No — keep"),
`Confounder` = ifelse(abs(pct_change) > 10, "Yes", "No")
)
})
assoc_sleep_table |>
kable(caption = "Task 4b — Systematic Confounder Assessment for Sleep (10% Rule)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) |>
column_spec(6, bold = TRUE)| Removed Covariate | Sleep β (max model) | Sleep β (without) | % Change | Within 10%? | Confounder |
|---|---|---|---|---|---|
| menthlth_days | -0.193 | -0.2894 | -50.0 | No — keep | Yes |
| age | -0.193 | -0.1646 | 14.7 | No — keep | Yes |
| sex | -0.193 | -0.1937 | -0.4 | Yes — drop | No |
| education | -0.193 | -0.1923 | 0.3 | Yes — drop | No |
| exercise | -0.193 | -0.1957 | -1.4 | Yes — drop | No |
| gen_health | -0.193 | -0.3593 | -86.2 | No — keep | Yes |
| income_cat | -0.193 | -0.1936 | -0.3 | Yes — drop | No |
| bmi | -0.193 | -0.1950 | -1.0 | Yes — drop | No |
Interpretation: The sleep coefficient in the maximum associative model reflects the effect of sleep on physically unhealthy days after adjusting for all other covariates. The 10% change-in-estimate rule is applied by removing each covariate individually and observing whether the sleep coefficient shifts outside the 10% boundary. Covariates that produce a percent change greater than 10% in absolute value are classified as confounders and must be retained in the final associative model. Those producing changes within the 10% window do not materially confound the sleep–physical health association and can be dropped for a more parsimonious final model.
4c. Final Associative Model
Fit the final model including only sleep_hrs and the
identified confounders. Report the sleep coefficient and 95% CI.
# Identify confirmed confounders (% change > 10%)
confounders_sleep <- assoc_sleep_table |>
filter(Confounder == "Yes") |>
pull(`Removed Covariate`)
cat("Confirmed confounders (|% change| > 10%):\n",
paste(" ", confounders_sleep, collapse = "\n"), "\n\n")## Confirmed confounders (|% change| > 10%):
## menthlth_days
## age
## gen_health# Build and fit the final associative model
final_assoc_formula <- as.formula(
paste("physhlth_days ~ sleep_hrs +", paste(confounders_sleep, collapse = " + "))
)
mod_final_assoc <- lm(final_assoc_formula, data = brfss_ms)
tidy(mod_final_assoc, conf.int = TRUE) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(
caption = "Task 4c — Final Associative Model Coefficients",
col.names = c("Term", "Estimate", "SE", "t", "p-value", "CI Lower", "CI Upper")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Term | Estimate | SE | t | p-value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 0.8151 | 0.5615 | 1.4516 | 0.1467 | -0.2857 | 1.9158 |
| sleep_hrs | -0.2026 | 0.0675 | -3.0003 | 0.0027 | -0.3349 | -0.0702 |
| menthlth_days | 0.1512 | 0.0120 | 12.5637 | 0.0000 | 0.1276 | 0.1748 |
| age | 0.0205 | 0.0054 | 3.7595 | 0.0002 | 0.0098 | 0.0312 |
| gen_healthVery good | 0.5113 | 0.2451 | 2.0860 | 0.0370 | 0.0308 | 0.9919 |
| gen_healthGood | 1.9151 | 0.2579 | 7.4255 | 0.0000 | 1.4095 | 2.4207 |
| gen_healthFair | 7.7686 | 0.3488 | 22.2693 | 0.0000 | 7.0847 | 8.4524 |
| gen_healthPoor | 21.4868 | 0.5266 | 40.8018 | 0.0000 | 20.4544 | 22.5192 |
# Extract the sleep coefficient and CI for reporting
sleep_final <- tidy(mod_final_assoc, conf.int = TRUE) |>
filter(term == "sleep_hrs")
cat("\nFinal sleep coefficient:", round(sleep_final$estimate, 4),
"\n95% CI: (", round(sleep_final$conf.low, 4), ",",
round(sleep_final$conf.high, 4), ")\n")##
## Final sleep coefficient: -0.2026
## 95% CI: ( -0.3349 , -0.0702 )Interpretation: After adjusting for the identified confounders, the sleep coefficient represents the independent association between sleep duration and physically unhealthy days. Each additional hour of sleep per night is associated with a statistically significant change in the number of physically unhealthy days in the past 30 days, holding the confounders constant. The 95% confidence interval excludes zero, confirming the association is unlikely to be due to chance. This adjusted estimate is preferred over the crude estimate because it accounts for variables that distort the sleep–physical health relationship.
4d. Response to Reviewer: Why Not Stepwise Selection?
Answer: Automated stepwise selection is
inappropriate for this associative analysis for several important
reasons. First and most critically, stepwise procedures treat all
variables — including our exposure of interest, sleep_hrs —
as candidates for removal; if sleep were not statistically significant
at some step, the algorithm might eliminate it, which would be
scientifically meaningless because our explicit goal is to
estimate the sleep effect, not to decide whether it “belongs”
in the model. Second, the 10% change-in-estimate approach we use is
driven by whether a covariate confounds the sleep coefficient —
a substantive criterion — whereas stepwise selection is driven purely by
whether a covariate improves statistical fit, which will sometimes
retain non-confounders and drop true confounders. Third, because
stepwise methods perform many sequential hypothesis tests on the same
data, the reported p-values and confidence intervals for the final model
are anti-conservative, making inference about the sleep coefficient
unreliable. In associative epidemiological research, the validity of the
exposure estimate — not the overall model fit — is the primary concern,
and that requires a theory-guided, confounder-focused approach rather
than an algorithmic one.
Task 5: Synthesis (20 points)
5a. Predictive vs. Associative Model Comparison
# Predictive best model = BIC-best from Task 2c (or backward from Task 3a — same)
# Associative best model = final model from Task 4c
# Display the two models side by side
cat("=== PREDICTIVE MODEL (BIC-best / Backward) ===\n")## === PREDICTIVE MODEL (BIC-best / Backward) ===## Formula: physhlth_days ~ menthlth_days + sleep_hrs + age + exercise + gen_health + income_cat## === ASSOCIATIVE MODEL (sleep + confounders) ===## Formula: physhlth_days ~ sleep_hrs + menthlth_days + age + gen_health# Summarize key differences
pred_terms <- attr(terms(mod_backward_lab), "term.labels")
assoc_terms <- attr(terms(mod_final_assoc), "term.labels")
cat("Variables in PREDICTIVE model only:", paste(setdiff(pred_terms, assoc_terms), collapse = ", "), "\n")## Variables in PREDICTIVE model only: exercise, income_catcat("Variables in ASSOCIATIVE model only:", paste(setdiff(assoc_terms, pred_terms), collapse = ", "), "\n")## Variables in ASSOCIATIVE model only:cat("Variables in BOTH models:", paste(intersect(pred_terms, assoc_terms), collapse = ", "), "\n\n")## Variables in BOTH models: menthlth_days, sleep_hrs, age, gen_health# Sleep coefficient in each model
sleep_pred <- tidy(mod_backward_lab) |> filter(term == "sleep_hrs") |> pull(estimate)
sleep_assoc <- tidy(mod_final_assoc) |> filter(term == "sleep_hrs") |> pull(estimate)
cat("Sleep coefficient — Predictive model: ", round(sleep_pred, 4), "\n")## Sleep coefficient — Predictive model: -0.1951## Sleep coefficient — Associative model: -0.2026# Fit statistics comparison
bind_rows(
glance(mod_backward_lab) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(Model = "Predictive (BIC-best/Backward)"),
glance(mod_final_assoc) |>
dplyr::select(r.squared, adj.r.squared, AIC, BIC) |>
mutate(Model = "Associative (sleep + confounders)")
) |>
dplyr::select(Model, everything()) |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
rename(`R²` = r.squared, `Adj. R²` = adj.r.squared) |>
kable(caption = "Task 5a — Predictive vs. Associative Model Comparison") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | R² | Adj. R² | AIC | BIC |
|---|---|---|---|---|
| Predictive (BIC-best/Backward) | 0.3857 | 0.3846 | 32638.37 | 32710.06 |
| Associative (sleep + confounders) | 0.3796 | 0.3787 | 32684.06 | 32742.71 |
Answer: The predictive and associative models do not include the exact same variables, and they differ in a conceptually important way.
The predictive model (from best subsets / backward AIC) includes all variables that improve out-of-sample prediction performance — mental health days, sleep, age, exercise, general health, and income — selected purely because they maximize predictive accuracy by the AIC/BIC criterion. It does not include sex, education, or BMI because those variables do not improve prediction enough to justify their inclusion.
The associative model includes
sleep_hrs as a fixed exposure (never removed regardless of
its p-value) plus only those covariates confirmed as confounders by the
10% change-in-estimate rule. The set of confounders may be a subset of
the predictive model’s variables, because not every good predictor of
the outcome is also a confounder of the specific sleep–physical health
association.
The sleep coefficient is similar across both models because both adjust for the most important confounders; however, it is not identical. The predictive model adjusts for all retained variables simultaneously, whereas the associative model adjusts only for confirmed confounders. Small differences in the sleep coefficient between the two models arise because the predictive model may include additional variables that are partial confounders or precision variables, and each set of adjusters shifts the estimate slightly. The similarity in the sleep coefficients across both models actually provides reassurance that the association is robust to different model specifications.
5b. Public Health Summary Paragraph
Answer:
Using data from a nationally representative survey of 5,000 U.S. adults, we examined the relationship between nightly sleep duration and the number of physically unhealthy days reported in the past month. After accounting for key factors that can distort this relationship — including mental health burden, age, general health status, physical activity, and household income — we found that each additional hour of sleep per night was associated with a meaningful reduction in physically unhealthy days. In other words, adults who sleep more tend to report fewer days each month when their physical health was “not good,” even after considering that people who are older, poorer, or in worse overall health may sleep differently. The strength of this association suggests that promoting adequate sleep could be a valuable component of public health strategies aimed at improving physical well-being. It is important to note, however, that because this is a cross-sectional study — meaning we measured sleep and health at the same point in time — we cannot determine whether poor sleep caused worse physical health, or whether having poor physical health caused people to sleep differently; longitudinal research would be needed to confirm a causal direction.
End of Lab Activity