Multiple Linear Regression

EPI 553 — Principles of Statistical Inference II (Spring 2026)

---

Introduction

In the previous lectures, we used Simple Linear Regression (SLR) to model a continuous outcome as a function of a single predictor. But real-world health phenomena are never driven by one variable alone. A person’s mental health is shaped by how much they sleep, how physically ill they are, how old they are, their income, their sex, and much more — simultaneously.

Multiple Linear Regression (MLR) extends SLR to accommodate this complexity. It allows us to:

• Adjust for confounding — estimate the association between an exposure and outcome after accounting for other covariates
• Improve prediction — use multiple pieces of information to better forecast an outcome
• Model effect modification — examine whether the effect of one variable changes across levels of another
• Gain precision — by explaining additional variance in the outcome, we can reduce uncertainty in our estimates

In epidemiology, MLR is our primary tool for multivariable analysis of continuous outcomes.

---

Setup and Data

library(tidyverse)
library(haven)
library(janitor)
library(knitr)
library(kableExtra)
library(plotly)
library(broom)
library(ggeffects)
library(gtsummary)
library(GGally)
library(car)
library(lmtest)
library(corrplot)
library(rlang)

The BRFSS 2020 Dataset

We continue using the Behavioral Risk Factor Surveillance System (BRFSS) 2020 — a nationally representative telephone survey of U.S. adults conducted by the CDC. This dataset contains over 400,000 respondents with detailed information on health behaviors, chronic conditions, and social determinants of health.

Research question for today:

What individual, behavioral, and socioeconomic factors are associated with the number of mentally unhealthy days in the past 30 days among U.S. adults?

This is a quintessential epidemiological question: we have a health outcome (mental health days) and want to understand its multivariable determinants while appropriately controlling for confounders.

brfss_class <- read_xpt(
  '/Users/samriddhi/Downloads/LLCP2020.XPT '
) %>%
  clean_names()

# ── Recode and clean ──────────────────────────────────────────────────────────
brfss_mlr <- brfss_class %>%
  mutate(
    # Outcome: mentally unhealthy days in past 30 (88 = none = 0; 77/99 = DK/refused = NA)
    menthlth_days = case_when(
      menthlth == 88              ~ 0,
      menthlth >= 1 & menthlth <= 30 ~ as.numeric(menthlth),
      TRUE                        ~ NA_real_
    ),
    # Physical health days (key predictor)
    physhlth_days = case_when(
      physhlth == 88              ~ 0,
      physhlth >= 1 & physhlth <= 30 ~ as.numeric(physhlth),
      TRUE                        ~ NA_real_
    ),
    # Sleep hours (practical cap at 14)
    sleep_hrs = case_when(
      sleptim1 >= 1 & sleptim1 <= 14 ~ as.numeric(sleptim1),
      TRUE                        ~ NA_real_
    ),
    # Age (capped at 80 per BRFSS coding)
    age = age80,
    # Income category (ordinal 1–8)
    income_cat = case_when(
      income2 %in% 1:8 ~ as.numeric(income2),
      TRUE             ~ NA_real_
    ),
    # Sex
    sex = factor(sexvar, levels = c(1, 2), labels = c("Male", "Female")),
    # Exercise in past 30 days (any physical activity)
    exercise = factor(case_when(
      exerany2 == 1 ~ "Yes",
      exerany2 == 2 ~ "No",
      TRUE          ~ NA_character_
    ), levels = c("No", "Yes")),
    # BMI (stored as integer × 100 in BRFSS)
    bmi = ifelse(bmi5 > 0, bmi5 / 100, NA_real_),
    # Income as labeled factor (for display)
    income_f = factor(income2, levels = 1:8,
      labels = c("<$10k", "$10-15k", "$15-20k", "$20-25k",
                 "$25-35k", "$35-50k", "$50-75k", ">$75k"))
  ) %>%
  filter(
    !is.na(menthlth_days),
    !is.na(physhlth_days),
    !is.na(sleep_hrs),
    !is.na(age),   age >= 18,
    !is.na(income_cat),
    !is.na(sex),
    !is.na(exercise)
  )

# ── Analytic sample (reproducible random sample) ────────────────────────────
set.seed(553)
brfss_mlr <- brfss_mlr %>%
  select(menthlth_days, physhlth_days, sleep_hrs, age,
         income_cat, income_f, sex, exercise, bmi) %>%
  slice_sample(n = 5000)

# Save for lab activity
saveRDS(brfss_mlr,
  "/Users/samriddhi/Downloads/brfss_mlr_2020.rds")

tibble(Metric = c("Observations", "Variables"),
       Value  = c(nrow(brfss_mlr), ncol(brfss_mlr))) %>%
  kable(caption = "Analytic Dataset Dimensions") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Analytic Dataset Dimensions

Metric	Value
Observations	5000
Variables	9

Descriptive Statistics

brfss_mlr %>%
  select(menthlth_days, physhlth_days, sleep_hrs, age,
         income_f, sex, exercise) %>%
  tbl_summary(
    label = list(
      menthlth_days ~ "Mentally unhealthy days (past 30)",
      physhlth_days ~ "Physically unhealthy days (past 30)",
      sleep_hrs     ~ "Sleep (hours/night)",
      age           ~ "Age (years)",
      income_f      ~ "Household income",
      sex           ~ "Sex",
      exercise      ~ "Any physical activity (past 30 days)"
    ),
    statistic = list(
      all_continuous()  ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 1,
    missing = "no"
  ) %>%
  add_n() %>%
  bold_labels() %>%
  modify_caption("**Table 1. Descriptive Statistics — BRFSS 2020 Analytic Sample (n = 5,000)**")

Table 1. Descriptive Statistics — BRFSS 2020 Analytic Sample (n = 5,000)

Characteristic	N	N = 5,0001
Mentally unhealthy days (past 30)	5,000	3.8 (7.7)
Physically unhealthy days (past 30)	5,000	3.3 (7.8)
Sleep (hours/night)	5,000	7.1 (1.3)
Age (years)	5,000	54.3 (17.2)
Household income	5,000
<$10k		190 (3.8%)
$10-15k		169 (3.4%)
$15-20k		312 (6.2%)
$20-25k		434 (8.7%)
$25-35k		489 (9.8%)
$35-50k		683 (14%)
$50-75k		841 (17%)
>$75k		1,882 (38%)
Sex	5,000
Male		2,331 (47%)
Female		2,669 (53%)
Any physical activity (past 30 days)	5,000	3,874 (77%)
1 Mean (SD); n (%)

p_hist <- ggplot(brfss_mlr, aes(x = menthlth_days)) +
  geom_histogram(binwidth = 1, fill = "steelblue", color = "white", alpha = 0.85) +
  labs(
    title    = "Distribution of Mentally Unhealthy Days in the Past 30 Days",
    subtitle = "BRFSS 2020 Analytic Sample (n = 5,000)",
    x        = "Number of Mentally Unhealthy Days",
    y        = "Count"
  ) +
  theme_minimal(base_size = 13)

ggplotly(p_hist)

Distribution of Mentally Unhealthy Days (BRFSS 2020)

Note: The outcome is right-skewed — most respondents report zero mentally unhealthy days, with a long tail toward 30. This is common in count outcomes and is worth keeping in mind when we evaluate model assumptions later.

---

Part 1: Guided Practice — Multiple Linear Regression

---

1. From Simple to Multiple: The Big Picture

1.1 Extending the Model

Recall the SLR model:

\[Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i\]

Multiple linear regression simply adds more predictors:

\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_p X_{pi} + \varepsilon_i\]

Or equivalently, using summation notation:

\[Y_i = \beta_0 + \sum_{j=1}^{p} \beta_j X_{ji} + \varepsilon_i\]

Symbol	Meaning
\(Y_i\)	Outcome for subject \(i\) (mentally unhealthy days)
\(\beta_0\)	Intercept — expected value of \(Y\) when all \(X\)s = 0
\(\beta_j\)	Partial regression coefficient for predictor \(j\)
\(X_{ji}\)	Value of predictor \(j\) for subject \(i\)
\(p\)	Number of predictors in the model
\(\varepsilon_i\)	Random error for subject \(i\), \(\varepsilon_i \overset{iid}{\sim} N(0,\sigma^2)\)

The word partial is critical. Each \(\beta_j\) represents the estimated change in \(E(Y)\) for a one-unit increase in \(X_j\), holding all other predictors constant. This “holding constant” is the mathematical mechanism by which we adjust for confounding.

1.2 Three Goals of Regression in Epidemiology

There are three analytic goals. Understanding which goal is driving your analysis determines how you build and interpret your model:

Goal	Question	Example
Descriptive	How do the \(X\)s describe \(Y\) in this sample?	What variables are correlated with mental health days in BRFSS 2020?
Predictive	Do the \(X\)s help predict \(Y\) statistically?	Can we build a model to identify people likely to report many mentally unhealthy days?
Associative (Etiologic)	What is the effect of a specific \(X\) on \(Y\), controlling for confounders?	What is the effect of sleep on mental health days, after adjusting for age, income, and physical health?

In epidemiology, we are most often doing associative analysis — estimating causal or causal-adjacent effects after controlling for confounding.

1.3 Why MLR is More Challenging than SLR

• Model selection is non-trivial — there may be many candidate predictors, and the “right” model depends on the scientific question.
• You cannot visualize the relationship in \(p\)-dimensional space. We work with 2D projections (partial regression plots, predicted vs. observed plots).
• Interpretation requires discipline — the “holding all else constant” framing must always accompany coefficient interpretation.
• Coefficients are computed numerically — unlike SLR, there are no simple closed-form formulas we can compute by hand for \(p > 2\).

---

2. Assumptions of Multiple Linear Regression

The assumptions are a direct extension of the SLR assumptions — often summarized with the LINE mnemonic:

	Assumption	In MLR
L	Linearity	\(E(Y)\) is a linear function of the \(\beta_j\)s. Non-linear functions of \(X\) (e.g., \(X^2\), \(\log X\)) are allowed, as long as the model is linear in the coefficients.
I	Independence	Observations \(Y_i\) are independent of each other.
N	Normality	Errors \(\varepsilon_i \sim N(0, \sigma^2)\). Required for inference (t-tests, CIs). Not required for OLS estimation itself.
E	Equal variance	\(\text{Var}(\varepsilon_i) = \sigma^2\) for all \(i\) (homoscedasticity).

Important clarifications:

• The linearity assumption refers to linearity in the parameters (\(\beta\)s), not in the predictors. Including \(X^2\) as a predictor doesn’t violate linearity.
• Normality is only needed for valid inference. With large samples (like BRFSS), the Central Limit Theorem means mild violations cause minimal bias.
• Homoscedasticity matters most for standard errors and inference. Mild departures are tolerable; severe heteroscedasticity requires correction (robust SEs, WLS).

# Pairs plot of continuous predictors vs outcome
brfss_mlr %>%
  select(menthlth_days, physhlth_days, sleep_hrs, age) %>%
  rename(
    `Mental Health\nDays`  = menthlth_days,
    `Physical Health\nDays` = physhlth_days,
    `Sleep\n(hrs)`         = sleep_hrs,
    Age                    = age
  ) %>%
  ggpairs(
    lower  = list(continuous = wrap("points", alpha = 0.05, size = 0.5)),
    diag   = list(continuous = wrap("densityDiag", fill = "steelblue", alpha = 0.5)),
    upper  = list(continuous = wrap("cor", size = 4)),
    title  = "Pairwise Relationships Among Key Variables (BRFSS 2020)"
  ) +
  theme_minimal(base_size = 11)

Visualizing Key Predictors vs. Mental Health Days

---

3. Missing Data

Before fitting any model, we must address missing data. Let’s expand on it with our BRFSS context.

3.1 Why Missing Data Matters

Statistical software performs complete case analysis by default: observations with any missing value on any variable in the model are dropped. This has two consequences:

1. Loss of sample size — reduces power
2. Potential bias — if the missing cases differ systematically from observed cases

The direction and magnitude of bias depends on the mechanism of missingness.

3.2 The Three Mechanisms

Mechanism	Definition	Bias?	Example in BRFSS
MCAR — Missing Completely at Random	P(missing) is unrelated to observed or unobserved data	No	Lab equipment randomly fails; a beaker is dropped
MAR — Missing at Random	P(missing) depends on observed data, not on the missing value itself	Potentially, but adjustable	Men are more likely to skip the income question; once you account for sex, the missingness is independent of income
MNAR — Missing Not at Random	P(missing) depends on the unobserved value itself	Yes, difficult to correct	People with very low income refuse to report income because of the stigma

MCAR is rarely a safe assumption with human subjects data. In BRFSS, income non-response is strongly socioeconomically patterned — a classic MNAR concern.

3.3 Implications for Practice

• For MCAR and MAR, multiple imputation (using the mice package in R) can recover unbiased estimates.
• For MNAR, there is no general statistical fix. You must acknowledge the limitation and reason about the direction of potential bias.
• In this course, we will use complete case analysis but will always note the amount and likely pattern of missingness.

# Assess missingness in our analytic variables before the sample was taken
brfss_class %>%
  select(menthlth, physhlth, sleptim1, age80, income2, sexvar, exerany2, bmi5) %>%
  summarise(across(everything(), ~ sum(is.na(.) | . %in% c(77, 88, 99)) / n() * 100)) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Pct_Missing_or_DK") %>%
  mutate(Pct_Missing_or_DK = round(Pct_Missing_or_DK, 1)) %>%
  arrange(desc(Pct_Missing_or_DK)) %>%
  kable(
    caption = "Table 2. Missing or Don't Know/Refused (%) — BRFSS 2020 Full Sample",
    col.names = c("Variable", "% Missing / DK / Refused")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 2. Missing or Don’t Know/Refused (%) — BRFSS 2020 Full Sample

Variable	% Missing / DK / Refused
physhlth	71.5
menthlth	65.6
income2	19.9
bmi5	10.3
age80	1.4
sleptim1	1.2
sexvar	0.0
exerany2	0.0

---

4. Fitting the MLR Model in R

4.1 Model Specification

We will build our model in stages — a common and pedagogically powerful approach that mirrors what you would do in a real epidemiological analysis.

Model 1 — Unadjusted (SLR baseline): Mentally unhealthy days ~ Physical health days

Model 2 — Sleep added: + Sleep hours

Model 3 — Full model: + Age + Income + Sex + Exercise

This sequential approach lets us see how each block of variables changes the coefficients and model fit.

# Model 1: Unadjusted
m1 <- lm(menthlth_days ~ physhlth_days, data = brfss_mlr)

# Model 2: Add sleep
m2 <- lm(menthlth_days ~ physhlth_days + sleep_hrs, data = brfss_mlr)

# Model 3: Full multivariable model
m3 <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age + income_cat + sex + exercise,
         data = brfss_mlr)

4.2 The lm() Formula Interface

In R, the lm() formula has the structure:

lm(outcome ~ predictor1 + predictor2 + predictor3, data = dataset)

• The outcome is to the left of the ~
• Predictors are on the right, separated by +
• R automatically creates dummy variables for factors
• The intercept is included by default (suppress with - 1 if needed)

4.3 Examining Model 3 Output

summary(m3)

## 
## Call:
## lm(formula = menthlth_days ~ physhlth_days + sleep_hrs + age + 
##     income_cat + sex + exercise, data = brfss_mlr)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.9192  -3.4262  -1.7803   0.2948  30.0568 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   12.475489   0.716959  17.401  < 2e-16 ***
## physhlth_days  0.291657   0.013579  21.478  < 2e-16 ***
## sleep_hrs     -0.509160   0.075348  -6.757 1.57e-11 ***
## age           -0.082307   0.005933 -13.872  < 2e-16 ***
## income_cat    -0.321323   0.052012  -6.178 7.02e-10 ***
## sexFemale      1.245053   0.202333   6.153 8.17e-10 ***
## exerciseYes   -0.342685   0.253138  -1.354    0.176    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.09 on 4993 degrees of freedom
## Multiple R-squared:  0.1569, Adjusted R-squared:  0.1559 
## F-statistic: 154.9 on 6 and 4993 DF,  p-value: < 2.2e-16

tidy(m3, conf.int = TRUE) %>%
  mutate(
    term     = dplyr::recode(term,
      "(Intercept)" = "Intercept",
      "physhlth_days" = "Physical health days",
      "sleep_hrs"     = "Sleep (hours/night)",
      "age"           = "Age (years)",
      "income_cat"    = "Income (ordinal 1-8)",
      "sexFemale"     = "Sex: Female (ref = Male)",
      "exerciseYes"   = "Exercise: Yes (ref = No)"
    ),
    across(where(is.numeric), ~ round(., 4))
  ) %>%
  kable(
    caption  = "Table 3. Multiple Linear Regression: Mentally Unhealthy Days ~ Multiple Predictors (BRFSS 2020, n = 5,000)",
    col.names = c("Term", "Estimate (β̂)", "Std. Error", "t-statistic",
                  "p-value", "95% CI Lower", "95% CI Upper")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE) %>%
  row_spec(c(2, 3), background = "#EBF5FB")  # highlight key predictors

Table 3. Multiple Linear Regression: Mentally Unhealthy Days ~ Multiple Predictors (BRFSS 2020, n = 5,000)

Term	Estimate (β̂	Std. Erro
Intercept	12.4755	0.7170	17.4006	0.0000	11.0699	13.8810
Physical health days	0.2917	0.0136	21.4779	0.0000	0.2650	0.3183
Sleep (hours/night)	-0.5092	0.0753	-6.7574	0.0000	-0.6569	-0.3614
Age (years)	-0.0823	0.0059	-13.8724	0.0000	-0.0939	-0.0707
Income (ordinal 1-8)	-0.3213	0.0520	-6.1778	0.0000	-0.4233	-0.2194
Sex: Female (ref = Male)	1.2451	0.2023	6.1535	0.0000	0.8484	1.6417
Exercise: Yes (ref = No)	-0.3427	0.2531	-1.3537	0.1759	-0.8389	0.1536

---

5. Interpreting Coefficients

This is the most important and most commonly mis-stated skill in applied regression. Let’s do it carefully.

5.1 The General Rule

For predictor \(X_j\): \(\hat{\beta}_j\) is the estimated change in the mean of \(Y\) for a one-unit increase in \(X_j\), holding all other predictors in the model constant.

The phrase “holding all other predictors constant” (also called ceteris paribus or partial effect) is what distinguishes MLR from SLR. It is the mathematical expression of confounder adjustment.

b <- round(coef(m3), 3)
ci <- round(confint(m3), 3)

5.2 The Fitted Equation

\[\widehat{\text{Mental Health Days}} = 12.475 + 0.292(\text{Phys. Health Days}) + -0.509(\text{Sleep}) + -0.082(\text{Age}) + -0.321(\text{Income}) + 1.245(\text{Female}) + -0.343(\text{Exercise: Yes})\]

5.3 Interpreting Each Term

Physical health days (\(\hat{\beta}\) = 0.292): Each additional day of poor physical health is associated with an estimated 0.292 additional mentally unhealthy day on average, holding sleep, age, income, sex, and exercise constant (95% CI: 0.265 to 0.318). This is the strongest predictor in the model and is consistent with the bidirectional mind-body connection documented in the literature.

Sleep hours (\(\hat{\beta}\) = -0.509): Each additional hour of sleep per night is associated with an estimated 0.509 fewer mentally unhealthy days on average, adjusting for all other covariates (95% CI: -0.657 to -0.361). The negative sign indicates a protective association.

Age (\(\hat{\beta}\) = -0.082): Each additional year of age is associated with 0.082 fewer mentally unhealthy days on average (holding all else constant). This counterintuitive finding is well-documented — older adults often report fewer mental health difficulties, possibly due to better emotion regulation, survivor bias, or cohort effects.

Income category (\(\hat{\beta}\) = -0.321): Each one-unit increase in the income category (on the 1–8 ordinal scale) is associated with 0.321 fewer mentally unhealthy days on average, consistent with the well-established socioeconomic gradient in mental health.

Sex: Female (\(\hat{\beta}\) = 1.245): Compared to males (the reference group), females report an estimated 1.245 more mentally unhealthy days on average, holding all other variables constant.

Exercise: Yes (\(\hat{\beta}\) = -0.343): People who engaged in any physical activity in the past 30 days report an estimated 0.343 fewer mentally unhealthy days compared to those who did not exercise, adjusting for all other covariates.

Intercept (\(\hat{\beta}_0\) = 12.475): The estimated mean mental health days for a person with zero physically unhealthy days, zero sleep hours, age = 0, income category = 0, who is male and does not exercise. This is a mathematical artifact — not a meaningful quantity in context.

Epidemiological insight: Notice that the physical health coefficient changes when we add other covariates. Compare \(\hat{\beta}\) from Model 1 vs. Model 3. This change represents the confounding that age, income, sleep, sex, and exercise collectively introduce into the crude physical health–mental health association.

# Compare the physhlth_days coefficient across models
tribble(
  ~Model, ~`Physical Health Days β̂`, ~`95% CI`, ~`Adj. R²`,
  "M1 (unadjusted)", round(coef(m1)[2], 3),
    paste0("(", round(confint(m1)[2,1],3), ", ", round(confint(m1)[2,2],3), ")"),
    round(summary(m1)$adj.r.squared, 3),
  "M2 (+sleep)",     round(coef(m2)[2], 3),
    paste0("(", round(confint(m2)[2,1],3), ", ", round(confint(m2)[2,2],3), ")"),
    round(summary(m2)$adj.r.squared, 3),
  "M3 (full)",       round(coef(m3)[2], 3),
    paste0("(", round(confint(m3)[2,1],3), ", ", round(confint(m3)[2,2],3), ")"),
    round(summary(m3)$adj.r.squared, 3)
) %>%
  kable(caption = "Table 4. Physical Health Days Coefficient Across Sequential Models") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE)

Table 4. Physical Health Days Coefficient Across Sequential Models

Model	Physical Health Days β	|95% CI	Adj. R
M1 (unadjusted)	0.312	(0.286, 0.338)	0.099
M2 (+sleep)	0.300	(0.274, 0.326)	0.110
M3 (full)	0.292	(0.265, 0.318)	0.156

---

6. Hypothesis Testing for Individual Coefficients

6.1 The t-Test for Each \(\hat{\beta}_j\)

For each predictor \(X_j\), we test:

\[H_0: \beta_j = 0 \quad \text{vs.} \quad H_A: \beta_j \neq 0\]

(given all other predictors are in the model)

The test statistic is:

\[t = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim t_{n-p-1}\]

where \(n\) is sample size and \(p\) is the number of predictors. We compare this to a \(t\)-distribution with \(n - p - 1\) degrees of freedom.

Critically: This tests whether \(X_j\) has a statistically significant partial association with \(Y\), given the other predictors already in the model. A predictor might be significant alone (unadjusted) but non-significant in the multivariable model due to overlap with other predictors (multicollinearity or confounding).

6.2 Interpreting the p-Values

tidy(m3, conf.int = TRUE) %>%
  mutate(
    Significance = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.1   ~ ".",
      TRUE            ~ "ns"
    ),
    term = dplyr::recode(term,
      "(Intercept)"    = "Intercept",
      "physhlth_days"  = "Physical health days",
      "sleep_hrs"      = "Sleep (hours/night)",
      "age"            = "Age (years)",
      "income_cat"     = "Income (ordinal 1-8)",
      "sexFemale"      = "Sex: Female",
      "exerciseYes"    = "Exercise: Yes"
    )
  ) %>%
  select(term, estimate, std.error, statistic, p.value, conf.low, conf.high, Significance) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption   = "Table 5. Hypothesis Tests for Each Partial Coefficient (Model 3)",
    col.names = c("Predictor", "β̂", "SE", "t", "p-value", "CI Lower", "CI Upper", "Sig.")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 5. Hypothesis Tests for Each Partial Coefficient (Model 3)

Predictor	β	S
Intercept	12.4755	0.7170	17.4006	0.0000	11.0699	13.8810	***
Physical health days	0.2917	0.0136	21.4779	0.0000	0.2650	0.3183	***
Sleep (hours/night)	-0.5092	0.0753	-6.7574	0.0000	-0.6569	-0.3614	***
Age (years)	-0.0823	0.0059	-13.8724	0.0000	-0.0939	-0.0707	***
Income (ordinal 1-8)	-0.3213	0.0520	-6.1778	0.0000	-0.4233	-0.2194	***
Sex: Female	1.2451	0.2023	6.1535	0.0000	0.8484	1.6417	***
Exercise: Yes	-0.3427	0.2531	-1.3537	0.1759	-0.8389	0.1536	ns

tidy(m3, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    term = dplyr::recode(term,
      "physhlth_days" = "Physical health days",
      "sleep_hrs"     = "Sleep (hours/night)",
      "age"           = "Age (years)",
      "income_cat"    = "Income (1-8)",
      "sexFemale"     = "Sex: Female",
      "exerciseYes"   = "Exercise: Yes"
    ),
    term = fct_reorder(term, estimate),
    sig  = ifelse(p.value < 0.05, "Significant (p < 0.05)", "Non-significant")
  ) %>%
  ggplot(aes(x = estimate, y = term, color = sig)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.25, linewidth = 0.9) +
  geom_point(size = 3.5) +
  scale_color_manual(values = c("Significant (p < 0.05)" = "steelblue",
                                "Non-significant" = "tomato")) +
  labs(
    title    = "Partial Regression Coefficients with 95% Confidence Intervals",
    subtitle = "Outcome: Mentally Unhealthy Days (BRFSS 2020, n = 5,000)",
    x        = "Estimated Change in Mental Health Days (β̂)",
    y        = NULL,
    color    = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

Coefficient Plot: Partial Associations with Mental Health Days

---

7. The Overall F-Test and ANOVA Decomposition

7.1 Partitioning Variability

Just as in SLR, MLR partitions total variability in \(Y\):

\[\underbrace{SSY}_{\text{Total}} = \underbrace{SSR}_{\text{Model explains}} + \underbrace{SSE}_{\text{Model doesn't explain}}\]

Source	Sum of Squares	df	Mean Square	F
Model (Regression)	\(SSR\)	\(p\)	\(MSR = SSR/p\)	\(F = MSR/MSE\)
Error (Residual)	\(SSE\)	\(n-p-1\)	\(MSE = SSE/(n-p-1)\)
Total	\(SSY\)	\(n-1\)

With \(p\) predictors, the model uses \(p\) degrees of freedom (one for each \(\hat{\beta}_j\)), leaving \(n - p - 1\) residual degrees of freedom.

7.2 The Overall F-Test

The omnibus F-test asks:

\[H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0 \quad \text{(none of the predictors matter)}\] \[H_A: \text{At least one } \beta_j \neq 0\]

\[F = \frac{MSR}{MSE} \sim F_{p,\; n-p-1} \quad \text{under } H_0\]

A large F-value (small p-value) means the predictors collectively explain more variation than expected by chance.

anova_m3 <- anova(m3)
print(anova_m3)

## Analysis of Variance Table
## 
## Response: menthlth_days
##                 Df Sum Sq Mean Sq  F value    Pr(>F)    
## physhlth_days    1  29475 29474.7 586.3411 < 2.2e-16 ***
## sleep_hrs        1   3323  3323.0  66.1052 5.345e-16 ***
## age              1   9261  9261.5 184.2385 < 2.2e-16 ***
## income_cat       1   2649  2648.8  52.6918 4.503e-13 ***
## sex              1   1918  1918.4  38.1626 7.025e-10 ***
## exercise         1     92    92.1   1.8326    0.1759    
## Residuals     4993 250993    50.3                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

glance(m3) %>%
  select(r.squared, adj.r.squared, sigma, statistic, p.value, df, df.residual, nobs) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  pivot_longer(everything(), names_to = "Statistic", values_to = "Value") %>%
  mutate(Statistic = dplyr::recode(Statistic,
    "r.squared"     = "R²",
    "adj.r.squared" = "Adjusted R²",
    "sigma"         = "Residual Std. Error (Root MSE)",
    "statistic"     = "F-statistic",
    "p.value"       = "p-value (overall F-test)",
    "df"            = "Model df (p)",
    "df.residual"   = "Residual df (n − p − 1)",
    "nobs"          = "n (observations)"
  )) %>%
  kable(caption = "Table 6. Overall Model Summary — Model 3") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 6. Overall Model Summary — Model 3

Statistic	Value
R²	0.1569
Adjusted R²	0.1559
Residual Std. Error (Root MSE)	7.0901
F-statistic	154.8953
p-value (overall F-test)	0.0000
Model df (p)	6.0000
Residual df (n − p − 1)	4993.0000
n (observations)	5000.0000

---

8. Coefficient of Determination: \(R^2\) and Adjusted \(R^2\)

8.1 \(R^2\)

\[R^2 = \frac{SSR}{SSY} = 1 - \frac{SSE}{SSY}\]

\(R^2\) is the proportion of total variance in \(Y\) explained by the model. It ranges from 0 to 1.

Critical property: \(R^2\) always increases (or stays the same) when you add a predictor, even if that predictor is random noise. It cannot be used to compare models with different numbers of predictors.

8.2 Adjusted \(R^2\)

\[R^2_{adj} = 1 - \frac{(1 - R^2)(n - 1)}{n - p - 1}\]

Adjusted \(R^2\) penalizes for the number of predictors \(p\). It only increases when a new predictor improves the model by more than chance alone.

• \(R^2_{adj} \leq R^2\) always
• \(R^2_{adj}\) can decrease if you add a useless predictor
• Use \(R^2_{adj}\) when comparing models with different numbers of predictors

tribble(
  ~Model, ~Predictors, ~R2, ~`Adj. R²`,
  "M1: Physical health only",      1, round(summary(m1)$r.squared, 4), round(summary(m1)$adj.r.squared, 4),
  "M2: + Sleep",                   2, round(summary(m2)$r.squared, 4), round(summary(m2)$adj.r.squared, 4),
  "M3: Full (+ age, income, sex, exercise)", 6,
    round(summary(m3)$r.squared, 4), round(summary(m3)$adj.r.squared, 4)
) %>%
  kable(caption = "Table 7. R² and Adjusted R² Across Sequential Models") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(3, bold = TRUE, background = "#EBF5FB")

Table 7. R² and Adjusted R² Across Sequential Models

Model	Predictors	R2	Adj. R²
M1: Physical health only	1	0.0990	0.0988
M2: + Sleep	2	0.1102	0.1098
M3: Full (+ age, income, sex, exercise)	6	0.1569	0.1559

Interpretation: R² for mental health outcomes in population surveys is typically low (5–20%), and that’s expected. Mental health is shaped by hundreds of unmeasured factors — genetics, life events, social support, trauma history. A low R² does not mean the model is “wrong” or “useless” for estimating associations. Statistical significance (the t-tests) and the magnitude and direction of coefficients matter more for etiologic questions.

8.3 Root MSE

\[RMSE = \sqrt{MSE} = s_{Y|X}\]

The Root MSE (also called the Residual Standard Error in R output) is the estimated standard deviation of the errors. It tells you, on average, how many days the model’s prediction is off by.

rmse_m3 <- round(summary(m3)$sigma, 2)
cat("Root MSE (Model 3):", rmse_m3, "mentally unhealthy days\n")

## Root MSE (Model 3): 7.09 mentally unhealthy days

---

9. Visualizing Predicted Values and Effects

9.1 Predicted vs. Observed

augment(m3) %>%
  ggplot(aes(x = .fitted, y = menthlth_days)) +
  geom_point(alpha = 0.08, size = 0.8, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linewidth = 1, linetype = "dashed") +
  geom_smooth(method = "loess", color = "darkblue", linewidth = 1, se = FALSE) +
  labs(
    title    = "Predicted vs. Observed Mentally Unhealthy Days",
    subtitle = "Dashed = perfect prediction; Blue = LOESS smoother",
    x        = "Fitted (Predicted) Values (Ŷ)",
    y        = "Observed Mentally Unhealthy Days (Y)"
  ) +
  theme_minimal(base_size = 13)

Predicted vs. Observed Mental Health Days (Model 3)

9.2 Marginal Effects Plots (ggeffects)

The ggeffects package computes adjusted predictions — the predicted mean of \(Y\) across values of one predictor, with all other predictors held at their means (or reference categories). This is a powerful visualization of what MLR is actually estimating.

ggpredict(m3, terms = "physhlth_days") %>%
  plot() +
  labs(
    title    = "Marginal Effect of Physical Health Days on Mentally Unhealthy Days",
    subtitle = "Other covariates held at their means/modes",
    x        = "Number of Physically Unhealthy Days (past 30)",
    y        = "Predicted Mentally Unhealthy Days"
  ) +
  theme_minimal(base_size = 13)

Adjusted Predictions: Physical Health Days → Mental Health Days

ggpredict(m3, terms = "sleep_hrs") %>%
  plot() +
  labs(
    title    = "Marginal Effect of Sleep on Mentally Unhealthy Days",
    subtitle = "Other covariates held at their means/modes",
    x        = "Sleep Hours per Night",
    y        = "Predicted Mentally Unhealthy Days"
  ) +
  theme_minimal(base_size = 13)

Adjusted Predictions: Sleep → Mental Health Days

ggpredict(m3, terms = c("exercise", "sex")) %>%
  plot() +
  labs(
    title    = "Predicted Mental Health Days by Exercise Status and Sex",
    subtitle = "Adjusted for physical health, sleep, age, income",
    x        = "Any Physical Activity in Past 30 Days",
    y        = "Predicted Mentally Unhealthy Days",
    color    = "Sex", fill = "Sex"
  ) +
  theme_minimal(base_size = 13)

Adjusted Predictions by Sex and Exercise

---

10. Residual Diagnostics

We assess the LINE assumptions through residual plots. Use plot(model) in base R or augment() + ggplot2.

par(mfrow = c(2, 2))
plot(m3, which = 1:4, col = adjustcolor("steelblue", alpha.f = 0.3), pch = 16)

Standard Residual Diagnostic Plots

par(mfrow = c(1, 1))

How to read each plot:

Plot	What to Look For	Assumption Tested
Residuals vs. Fitted	Random scatter around 0; no pattern	Linearity + Homoscedasticity
Normal Q-Q	Points follow diagonal line	Normality of residuals
Scale-Location (√|e|)	Flat red line; constant spread	Homoscedasticity
Cook’s Distance / Residuals vs. Leverage	No points beyond 0.5 or 1.0 Cook’s D contour	Influential observations

augment(m3) %>%
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  geom_point(alpha = 0.07, size = 0.8, color = "steelblue") +
  geom_smooth(method = "loess", color = "darkblue", se = FALSE) +
  labs(
    title = "Residuals vs. Fitted Values",
    subtitle = "Ideally: random scatter around zero, no systematic pattern",
    x = "Fitted Values",
    y = "Residuals (Y − Ŷ)"
  ) +
  theme_minimal(base_size = 13)

Residuals vs. Fitted (ggplot2)

Note on our outcome: Because mental health days is bounded at 0 and 30 and is heavily right-skewed, perfect normality and homoscedasticity are not achievable. Our large sample size (n = 5,000) provides robustness. In practice, you might consider a log transformation or Poisson/negative binomial regression as alternatives — topics we’ll cover later in the course.

---

11. Main Effects and Interaction Terms

11.1 Higher-Order Terms

The MLR framework accommodates higher-order terms (polynomials, interactions) while remaining “linear in the parameters.” Some key rules:

• Quadratic terms: If you include \(X^2\), you must also include \(X\) itself.
• Interaction terms: If you include \(X_1 \times X_2\), you must include both \(X_1\) and \(X_2\) as main effects.

These are called the hierarchical principle or marginality rule in regression modeling.

11.2 Testing for an Interaction

Does the effect of sleep on mental health differ between males and females?

# Model with sleep × sex interaction
m4_interact <- lm(menthlth_days ~ physhlth_days + sleep_hrs * sex + age + income_cat + exercise,
                  data = brfss_mlr)

tidy(m4_interact, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Table 8. Model with Sleep × Sex Interaction",
    col.names = c("Term", "β̂", "SE", "t", "p-value", "CI Low", "CI High")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 8. Model with Sleep × Sex Interaction

Term	β	S
(Intercept)	12.3812	0.9203	13.4537	0.0000	10.5770	14.1853
physhlth_days	0.2917	0.0136	21.4763	0.0000	0.2651	0.3183
sleep_hrs	-0.4960	0.1101	-4.5072	0.0000	-0.7118	-0.2803
sexFemale	1.4178	1.0757	1.3180	0.1876	-0.6911	3.5267
age	-0.0823	0.0059	-13.8720	0.0000	-0.0940	-0.0707
income_cat	-0.3210	0.0521	-6.1673	0.0000	-0.4231	-0.2190
exerciseYes	-0.3422	0.2532	-1.3518	0.1765	-0.8386	0.1541
sleep_hrs:sexFemale	-0.0244	0.1495	-0.1635	0.8701	-0.3174	0.2686

ggpredict(m4_interact, terms = c("sleep_hrs", "sex")) %>%
  plot() +
  labs(
    title    = "Predicted Mental Health Days by Sleep Hours and Sex",
    subtitle = "Interaction model: slopes are allowed to differ by sex",
    x        = "Sleep Hours per Night",
    y        = "Predicted Mentally Unhealthy Days",
    color    = "Sex", fill = "Sex"
  ) +
  theme_minimal(base_size = 13)

Interaction: Effect of Sleep on Mental Health by Sex

# Likelihood ratio / F-test comparing models with and without interaction
anova(m3, m4_interact) %>%
  kable(caption = "Table 9. F-Test Comparing Models With and Without Interaction") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Table 9. F-Test Comparing Models With and Without Interaction

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
4993	250992.8	NA	NA	NA	NA
4992	250991.5	1	1.344174	0.0267344	0.8701261

Interpretation: If the interaction term is statistically significant, the slopes (effect of sleep) differ between males and females. If not significant, a model with only main effects is preferred (parsimony). We’ll cover interaction/effect modification formally in a dedicated lecture.

---

12. The Multivariate Normal Distribution (Conceptual Bridge)

Just as SLR is linked to the bivariate normal distribution, MLR is linked to the multivariate normal distribution. When all continuous variables follow a multivariate normal distribution, the conditional distribution of \(Y\) given any combination of \(X\)s is:

\[Y \mid X_1, X_2, \ldots, X_p \sim N(\beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p, \; \sigma^2_{Y|X})\]

This is the population model that our sample estimates are trying to recover. The bivariate normal from the correlation lecture is simply a special case with \(p = 1\).

---

13. Summary of Key Takeaways

Concept	Key Point
MLR formula	\(Y_i = \beta_0 + \sum_{j=1}^p \beta_j X_{ji} + \varepsilon_i\)
Partial coefficient	\(\hat{\beta}_j\) = change in \(E(Y)\) per 1-unit increase in \(X_j\), holding all else constant
Confounder adjustment	Achieved mathematically by including confounders as additional predictors
Assumptions	LINE: Linearity, Independence, Normality, Equal variance
Missing data	MCAR, MAR, MNAR — mechanism determines bias risk
F-test	Tests \(H_0: \beta_1 = \cdots = \beta_p = 0\) (omnibus test)
\(R^2\) vs. Adj. \(R^2\)	Use Adj. \(R^2\) to compare models of different sizes
Root MSE	Average prediction error in outcome units

---

---

Part 2: In-Class Lab Activity

EPI 553 — Multiple Linear Regression Lab Due: End of class, March 10, 2026

---

Instructions

In this lab, you will build and interpret multiple linear regression models using the BRFSS 2020 analytic 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. It contains 5,000 randomly sampled BRFSS 2020 respondents with the following variables:

Variable	Description	Type
menthlth_days	Mentally unhealthy days in past 30	Continuous (0–30)
physhlth_days	Physically 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
income_cat	Household income (1 = <$10k to 8 = >$75k)	Ordinal
sex	Sex (Male/Female)	Factor
exercise	Any physical activity past 30 days (Yes/No)	Factor

# Load the dataset
library(tidyverse)
library(broom)
library(knitr)
library(kableExtra)
library(gtsummary)
library(ggeffects)

brfss_mlr <- readRDS(
  "/Users/samriddhi/Downloads/brfss_mlr_2020.rds"
)

---

Task 1: Exploratory Data Analysis (15 points)

1a. (5 pts) Create a descriptive statistics table using tbl_summary() that includes all variables in the dataset. Include means (SD) for continuous variables and n (%) for categorical variables.

brfss_mlr %>%
  select(menthlth_days, physhlth_days, sleep_hrs, age, income_cat,
         income_f, sex, exercise, bmi) %>%
  tbl_summary(
    label = list(
      menthlth_days ~ "Mentally unhealthy days (past 30)",
      physhlth_days ~ "Physically unhealthy days (past 30)",
      sleep_hrs     ~ "Sleep (hours/night)",
      age           ~ "Age (years)",
      income_cat    ~ "Income category",
      income_f      ~ "Household income",
      sex           ~ "Sex",
      exercise      ~ "Any physical activity (past 30 days)",
      bmi           ~  "Body Mass Index (kg/m2)"
    ),
    statistic = list(
      all_continuous()  ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 1,
    missing = "no"
  ) %>%
  add_n() %>%
  bold_labels() %>%
  modify_caption("**Table 1. Descriptive Statistics — BRFSS 2020 Analytic Sample (n = 5,000)**")

Table 1. Descriptive Statistics — BRFSS 2020 Analytic Sample (n = 5,000)

Characteristic	N	N = 5,0001
Mentally unhealthy days (past 30)	5,000	3.8 (7.7)
Physically unhealthy days (past 30)	5,000	3.3 (7.8)
Sleep (hours/night)	5,000	7.1 (1.3)
Age (years)	5,000	54.3 (17.2)
Income category	5,000
1		190 (3.8%)
2		169 (3.4%)
3		312 (6.2%)
4		434 (8.7%)
5		489 (9.8%)
6		683 (14%)
7		841 (17%)
8		1,882 (38%)
Household income	5,000
<$10k		190 (3.8%)
$10-15k		169 (3.4%)
$15-20k		312 (6.2%)
$20-25k		434 (8.7%)
$25-35k		489 (9.8%)
$35-50k		683 (14%)
$50-75k		841 (17%)
>$75k		1,882 (38%)
Sex	5,000
Male		2,331 (47%)
Female		2,669 (53%)
Any physical activity (past 30 days)	5,000	3,874 (77%)
Body Mass Index (kg/m2)	4,706	28.4 (6.4)
1 Mean (SD); n (%)

1b. (5 pts) Create a histogram of menthlth_days. Describe the shape of the distribution. Is it symmetric, right-skewed, or left-skewed? What are the implications of this shape for regression modeling? The histogram plot is not symmetric and it is rightly skewed. We can see that the data points are on the extreme right side of the curve towards 30. It is important to consider it when evaluating the model assumption later.

p_hist <- ggplot(brfss_mlr, aes(x = menthlth_days)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "white", alpha = 0.85) +
  labs(
    title    = "Distribution of Mentally Unhealthy Days in the Past 30 Days",
    subtitle = "BRFSS 2020 Analytic Sample (n = 5,000)",
    x        = "Number of Mentally Unhealthy Days",
    y        = "Count"
  ) +
  theme_minimal(base_size = 13)

ggplotly(p_hist)

1c. (5 pts) Create a scatterplot matrix (using GGally::ggpairs() or similar) for the continuous variables: menthlth_days, physhlth_days, sleep_hrs, and age. Comment on the direction and strength of each pairwise correlation with the outcome. Mental health days and sleep is weekly negative correlation with age. This means every unit increase in age the mental healthy days and sleep decreases by 0.156 days and 0.140 hours respectively. While physical health is moderately positive correlation with mental health days, which implies that every unit increase in age the mental healthy days increases by 0.135

# Pairs plot of continuous predictors vs outcome
brfss_mlr %>%
  select(menthlth_days, physhlth_days, sleep_hrs, age) %>%
  rename(
    `Mental Health\nDays`  = menthlth_days,
    `Physical Health\nDays` = physhlth_days,
    `Sleep\n(hrs)`         = sleep_hrs,
    Age                    = age
  ) %>%
  ggpairs(
    lower  = list(continuous = wrap("points", alpha = 0.05, size = 0.5)),
    diag   = list(continuous = wrap("densityDiag", fill = "blue", alpha = 0.5)),
    upper  = list(continuous = wrap("cor", size = 4)),
    title  = "Pairwise Relationships Among Key Variables (BRFSS 2020)"
  ) +
  theme_minimal(base_size = 11)

---

Task 2: Unadjusted (Simple) Linear Regression (15 points)

2a. (5 pts) Fit a simple linear regression model regressing menthlth_days on sleep_hrs alone.

#fitted regression equation
model_mlr <- lm(menthlth_days ~ sleep_hrs, data = brfss_mlr)

# Summary output
summary(model_mlr)

## 
## Call:
## lm(formula = menthlth_days ~ sleep_hrs, data = brfss_mlr)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.670 -3.845 -3.040 -0.040 31.785 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.47429    0.57712   16.42   <2e-16 ***
## sleep_hrs   -0.80424    0.08025  -10.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.642 on 4998 degrees of freedom
## Multiple R-squared:  0.0197, Adjusted R-squared:  0.0195 
## F-statistic: 100.4 on 1 and 4998 DF,  p-value: < 2.2e-16

#coefficient table
tidy(model_mlr, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Simple Linear Regression: menthlth_days ~ sleep_hrs (BRFSS 2020)",
    col.names = c("Term", "Estimate", "Std. Error", "t-statistic",
                  "p-value", "95% CI Lower", "95% CI Upper"),
    align = "lrrrrrrr"
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE)

Simple Linear Regression: menthlth_days ~ sleep_hrs (BRFSS 2020)

Term	Estimate	Std. Error	t-statistic	p-value	95% CI Lower	95% CI Upper
(Intercept)	9.4743	0.5771	16.4165	0	8.3429	10.6057
sleep_hrs	-0.8042	0.0802	-10.0218	0	-0.9616	-0.6469

Write out the fitted regression equation. Poor mental health=9.4743+(−0.8042)× sleep hours

2b. (5 pts) Interpret the slope for sleep in a single sentence appropriate for a public health audience (no statistical jargon). Slope (𝑏1=−0.8042): Each additional hour of sleep was associated with approximately 0.8 fewer poor mental health days per month, on average, holding all else constant (though there are no other variables in this simple model).

2c. (5 pts) State the null and alternative hypotheses for the slope, report the t-statistic and p-value, and state your conclusion. What is the degree of freedom for this test? 𝐻0:𝛽1=0 no linear relationship between poor mental health and sleep hours 𝐻𝐴:𝛽1≠0 there is a linear relationship between poor mental health and sleep hours t-statistic: −10.02 p-value: < 2x 10(-16) Because the p-value is less than 0.05, we reject the null hypothesis. This suggests that sleep hours are associated with the number of mentally unhealthy days. Specifically, people who sleep more hours tend to report fewer mentally unhealthy days, with about 0.8 fewer days for each additional hour of sleep on average. Given n = 5000: df=5000-2=4998

---

Task 3: Building the Multivariable Model (25 points)

3a. (5 pts) Fit three models:

• Model A: menthlth_days ~ sleep_hrs
• Model B: menthlth_days ~ sleep_hrs + age + sex
• Model C: menthlth_days ~ sleep_hrs + age + sex + physhlth_days + income_cat + exercise

# Model 1: Unadjusted
m1 <- lm(menthlth_days ~ sleep_hrs, data = brfss_mlr)

# Model 2: Add sleep
m2 <- lm(menthlth_days ~  sleep_hrs + age + sex, data = brfss_mlr)

# Model 3: Full multivariable model
m3 <- lm(menthlth_days ~sleep_hrs + age + sex + physhlth_days + income_cat + exercise,
         data = brfss_mlr)

summary(m3)

## 
## Call:
## lm(formula = menthlth_days ~ sleep_hrs + age + sex + physhlth_days + 
##     income_cat + exercise, data = brfss_mlr)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.9192  -3.4262  -1.7803   0.2948  30.0568 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   12.475489   0.716959  17.401  < 2e-16 ***
## sleep_hrs     -0.509160   0.075348  -6.757 1.57e-11 ***
## age           -0.082307   0.005933 -13.872  < 2e-16 ***
## sexFemale      1.245053   0.202333   6.153 8.17e-10 ***
## physhlth_days  0.291657   0.013579  21.478  < 2e-16 ***
## income_cat    -0.321323   0.052012  -6.178 7.02e-10 ***
## exerciseYes   -0.342685   0.253138  -1.354    0.176    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.09 on 4993 degrees of freedom
## Multiple R-squared:  0.1569, Adjusted R-squared:  0.1559 
## F-statistic: 154.9 on 6 and 4993 DF,  p-value: < 2.2e-16

3b. (10 pts) Create a table comparing the sleep coefficient (\(\hat{\beta}\), SE, 95% CI, p-value) across Models A, B, and C. Does the sleep coefficient change substantially when you add more covariates? What does this suggest about confounding? Yes, the sleep coefficient changes when additional covariates are added, but the change is not very large. In the previous regression model, each additional hour of sleep was associated with 0.80 fewer mentally unhealthy days. After adjusting for other covariates (age, sex, physical unhealthy days, income, and exercise) in the multiple regression model, the mentally unhealthy days decreased to 0.51.This reduction suggests that some confounding is present, meaning that factors such as age, income, and physical health partially explain the relationship between sleep and mentally unhealthy days. However, the association remains statistically significant, indicating that sleep independently contributes to mental health outcomes.

tidy(m3, conf.int = TRUE) %>%
  mutate(
    term     = dplyr::recode(term,
      "(Intercept)" = "Intercept",
      "menthlth_days" = "Poor  mental health days",
      "sleep_hrs"     = "Sleep (hours/night)",
      "age"           = "Age (years)",
      "physhlth_days" = "Physical unealthy days",
      "income_cat"    = "Income (ordinal 1-8)",
      "sexFemale"     = "Sex: Female (ref = Male)",
      "exerciseYes"   = "Exercise: Yes (ref = No)"
    ),
    across(where(is.numeric), ~ round(., 4))
  ) %>%
  kable(
    caption  = "Table 3. Multiple Linear Regression: Mentally Unhealthy Days ~ Multiple Predictors (BRFSS 2020, n = 5,000)",
    col.names = c("Term", "Estimate (β̂)", "Std. Error", "t-statistic",
                  "p-value", "95% CI Lower", "95% CI Upper")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE) %>%
  row_spec(c(2, 3), background = "#EBF5FB")  # highlight key predictors

Table 3. Multiple Linear Regression: Mentally Unhealthy Days ~ Multiple Predictors (BRFSS 2020, n = 5,000)

Term	Estimate (β̂	Std. Erro
Intercept	12.4755	0.7170	17.4006	0.0000	11.0699	13.8810
Sleep (hours/night)	-0.5092	0.0753	-6.7574	0.0000	-0.6569	-0.3614
Age (years)	-0.0823	0.0059	-13.8724	0.0000	-0.0939	-0.0707
Sex: Female (ref = Male)	1.2451	0.2023	6.1535	0.0000	0.8484	1.6417
Physical unealthy days	0.2917	0.0136	21.4779	0.0000	0.2650	0.3183
Income (ordinal 1-8)	-0.3213	0.0520	-6.1778	0.0000	-0.4233	-0.2194
Exercise: Yes (ref = No)	-0.3427	0.2531	-1.3537	0.1759	-0.8389	0.1536

3c. (10 pts) For Model C, write out the full fitted regression equation and interpret every coefficient in plain language appropriate for a public health report. Mentally Unhealthy Days =12.48+−0.51(Sleep Hours)+−0.08(Age)+1.25(Female)+0.29(Physically Unhealthy Days)−0.32(Income)+−0.34(Exercise) Where:Female = 1 if female, 0 if male Exercise = 1 if exercises, 0 if does not exercise Income is an ordinal scale from 1 (lowest) to 8 (highest).

Sleep hours (𝛽̂ = -0.509): Each additional hour of sleep per night is associated with an estimated 0.509 fewer mentally unhealthy days on average, adjusting for covariates like age, sex, physical unhealthy days, income category and excersize (95% CI: -0.657 to -0.361). The negative sign indicates a protective association. Age (𝛽̂ = -0.082): Each additional year of age is associated with 0.082 (95% CI: -0.094 to -0.071) fewer mentally unhealthy days on average (holding all else constant). This finding is well-documented — older adults often report fewer mental health difficulties, possibly due to better emotion regulation, survivor bias, or cohort effects. Sex: Female (𝛽= 1.245): Compared to males (the reference group), females report an estimated 1.25 (95% CI:0.85 to 1.64) more mentally unhealthy days on average, holding all other variables constant. Physical unhealthy days (𝛽̂ = 0.292): Each additional day of poor physical health is associated with an estimated 0.292 (95% CI: 0.27 to 0.32) additional mentally unhealthy day on average. Income category (𝛽̂ = -0.321): Each one-unit increase in the income category (on the 1–8 ordinal scale) is associated with 0.321 (CI 95%: -0.42 to -0.22) fewer mentally unhealthy days on average, consistent with the well-established socioeconomic gradient in mental health. Exercise:(𝛽̂ = -0.343): People who engaged in any physical activity in the past 30 days report an estimated 0.343 (95% CI: -0.84 to 0.15) fewer mentally unhealthy days compared to those who did not exercise.

The model suggests that sleep duration, age, sex, physical health, and income are significant predictors of mentally unhealthy days. Individuals who sleep more, are older, and have higher income tend to report fewer mentally unhealthy days, while women and individuals experiencing more physical health problems report more mentally unhealthy days.

Task 4: Model Fit and ANOVA (20 points)

4a. (5 pts) Report \(R^2\) and Adjusted \(R^2\) for each of the three models (A, B, C). Create a table. Which model explains the most variance in mental health days? Model 3 has the highest Adjusted R² (0.156), meaning it explains about 15.6% of the variability in mentally unhealthy days. The fully adjusted model (Model 3) explains the most variance in mental health days and therefore provides the best overall fit among the three models.

# Compare the menthlth_days coefficient across models
tribble(
  ~Model, ~`menthlth_days β̂`, ~`95% CI`, ~`Adj. R²`,
  "M1 (unadjusted)", round(coef(m1)[2], 3),
    paste0("(", round(confint(m1)[2,1],3), ", ", round(confint(m1)[2,2],3), ")"),
    round(summary(m1)$adj.r.squared, 3),
  "M2 (+sleep)",     round(coef(m2)[2], 3),
    paste0("(", round(confint(m2)[2,1],3), ", ", round(confint(m2)[2,2],3), ")"),
    round(summary(m2)$adj.r.squared, 3),
  "M3 (full)",       round(coef(m3)[2], 3),
    paste0("(", round(confint(m3)[2,1],3), ", ", round(confint(m3)[2,2],3), ")"),
    round(summary(m3)$adj.r.squared, 3)
) %>%
  kable(caption = "Table 4. Mental UnHealthy Days Coefficient Across Sequential Models") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE)

Table 4. Mental UnHealthy Days Coefficient Across Sequential Models

Model	menthlth_days β	|95% CI	Adj. R
M1 (unadjusted)	-0.804	(-0.962, -0.647)	0.020
M2 (+sleep)	-0.734	(-0.889, -0.578)	0.050
M3 (full)	-0.509	(-0.657, -0.361)	0.156

4b. (5 pts) What is the Root MSE for Model C? Interpret it in practical terms — what does it tell you about prediction accuracy?

anova_m3 <- anova(m3)
print(anova_m3)

## Analysis of Variance Table
## 
## Response: menthlth_days
##                 Df Sum Sq Mean Sq  F value    Pr(>F)    
## sleep_hrs        1   5865  5864.8 116.6678 < 2.2e-16 ***
## age              1   6182  6182.2 122.9832 < 2.2e-16 ***
## sex              1   2947  2947.1  58.6266 2.274e-14 ***
## physhlth_days    1  29456 29455.5 585.9585 < 2.2e-16 ***
## income_cat       1   2177  2176.8  43.3031 5.169e-11 ***
## exercise         1     92    92.1   1.8326    0.1759    
## Residuals     4993 250993    50.3                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

glance(m3) %>%
  select(r.squared, adj.r.squared, sigma, statistic, p.value, df, df.residual, nobs) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  pivot_longer(everything(), names_to = "Statistic", values_to = "Value") %>%
  mutate(Statistic = dplyr::recode(Statistic,
    "r.squared"     = "R²",
    "adj.r.squared" = "Adjusted R²",
    "sigma"         = "Residual Std. Error (Root MSE)",
    "statistic"     = "F-statistic",
    "p.value"       = "p-value (overall F-test)",
    "df"            = "Model df (p)",
    "df.residual"   = "Residual df (n − p − 1)",
    "nobs"          = "n (observations)"
  )) %>%
  kable(caption = "Table 6. Overall Model Summary — Model C") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 6. Overall Model Summary — Model C

Statistic	Value
R²	0.1569
Adjusted R²	0.1559
Residual Std. Error (Root MSE)	7.0901
F-statistic	154.8953
p-value (overall F-test)	0.0000
Model df (p)	6.0000
Residual df (n − p − 1)	4993.0000
n (observations)	5000.0000

4c. (10 pts) Using the ANOVA output for Model C, fill in the following table manually (i.e., compute the values using the output from anova() or glance()):

Source	df	SS	MS	F
Model	6	46719	7786.50	154.86
Residual	4993	250993	50.3
Total	4999	297712

State the null hypothesis for the overall F-test and your conclusion. The null hypothesis states that all regression coefficients are equal to zero, meaning none of the predictors are associated with mentally unhealthy days. Since the overall F-test is statistically significant (p value=2.2 × 10⁻¹⁶), we reject the null hypothesis and conclude that at least one predictor in the model is significantly associated with mentally unhealthy days. —

Task 5: Residual Diagnostics (15 points)

5a. (5 pts) For Model C, produce the four standard diagnostic plots (Residuals vs. Fitted, Normal Q-Q, Scale-Location, Cook’s Distance). Comment on what each plot tells you about the LINE assumptions. Residuals vs Fitted Plot: The red line should be relatively flat and should stay near the 0. In the residual vs fitted plot, The residuals also shows a downward trend. There is a slight violation of the linearity.Therefore, we will conclude the linearity and equal variance assumptions may not be fully satisfied. Q-Q Plot (Normality of Residuals): For this plot we check that the data points should follow the dashed diagonal line. Here, the data points follow the diagonal line. However, there is a deviation at the upper tail, suggesting the data is rightly skewed. Therefore we conclude that residuals are not normally distributed. Also, the data is skewed rightly in the distribution. Scale–Location Plot (Equal Variance): This plot checks if the “spread” across the fitted values, the red line should be flat. For the plot below, the red line increases with fitted values.This indicates heteroscedasticity, meaning the variance of residuals is not constant. Thus, we conclude that the equal variance assumption is violated. Cook’s Distance Plot: Most observations should have very small Cook’s distance values. Most of the observations are close toCook’s distance. There are only few data points which are away from the cook’s distance. Hence, we conclude there are some potentially influential observations, but they do not appear to severely distort the model.

par(mfrow = c(2, 2))
plot(m3, which = 1:4, col = adjustcolor("blue", alpha.f = 0.3), pch = 16)

par(mfrow = c(1, 1))

5b. (5 pts) Given the nature of the outcome (mental health days, bounded at 0 and 30, heavily right-skewed), which assumptions are most likely to be violated? Does this invalidate the analysis? Explain. The assumptions most likely violated are normality of residuals and constant variance (homoscedasticity) due to the bounded, right-skewed count outcome. This does not completely invalidate the regression, especially with a large sample.

5c. (5 pts) Identify any observations with Cook’s Distance > 1. How many are there? What would you do with them in a real analysis? There are few observations which are distant from the crooks line. These are three of them (830, 1439, 1632).These values are far below the threshold of 1, meaning they are not strongly influential on the regression model. Hence, its is not required to remove these values.

---

Task 6: Interpretation for a Public Health Audience (10 points)

Suppose you are writing the results section of a public health paper. Write a 3–4 sentence paragraph summarizing the findings from Model C for a non-statistical audience. Your paragraph should:

• Identify which predictors were significantly associated with mental health days
• State the direction and approximate magnitude of the most important associations
• Appropriately caveat the cross-sectional nature of the data (no causal language)
• Not use any statistical jargon (no “significant,” “coefficient,” “p-value”)

Based on the results in Table 3, several factors were related to the number of mentally unhealthy days reported in the past month.

Based upon the study results. People who slept more hours per night tended to report fewer mentally unhealthy days. On average, each additional hour of sleep was linked with about half a day fewer of poor mental health in a month.

Age was also related to mental health days. Older adults reported fewer mentally unhealthy days than younger adults. The difference was relatively small for each additional year of age, but it could add up across many years.

Sex showed a difference as well. Women reported about 1.2 more mentally unhealthy days per month than men, on average.

The strongest relationship was with physical unhealthy days. People who reported more days of poor physical health also tended to report more mentally unhealthy days. For each additional day of poor physical health, there were about 0.29 more mentally unhealthy days. Income was also related to mental health. People in higher income categories tended to report fewer mentally unhealthy days, with roughly 0.3 fewer days for each step up in the income category. In contrast, exercise was not clearly related to mentally unhealthy days in this analysis.

It is important to note that these data come from a cross-sectional survey, meaning all information was collected at one snap shot of time. Because of this, the results cannot determine cause and effect. The findings only show that these factors tended to occur together with differences in mentally unhealthy days among the survey participants.

---

End of Lab Activity

---