R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Data for the Lab

Use the same BRFSS 2020 dataset from the guided practice.

# Load packages and data
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.5.2
## Warning: package 'tibble' was built under R version 4.5.2
## Warning: package 'tidyr' was built under R version 4.5.2
## Warning: package 'readr' was built under R version 4.5.2
## Warning: package 'purrr' was built under R version 4.5.2
## Warning: package 'dplyr' was built under R version 4.5.2
## Warning: package 'stringr' was built under R version 4.5.2
## Warning: package 'forcats' was built under R version 4.5.2
## Warning: package 'lubridate' was built under R version 4.5.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom)
library(knitr)
## Warning: package 'knitr' was built under R version 4.5.2
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.5.2
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(car)
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.5.2
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(gtsummary)
## Warning: package 'gtsummary' was built under R version 4.5.2
brfss_mlr <- readRDS(
  "C:\\Users\\safwa\\OneDrive - University at Albany - SUNY\\EPI 553\\Labs\\brfss dataset\\brfss_mlr_2020.rds")

Round 1: Student A Drives, Student B Navigates

Task 1: Fitting Models and ANOVA Tables (15 points)

1a. (5 pts) Fit the following model:

#1a
# Fit the three-predictor model
m_task1 <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age,
              data = brfss_mlr)

# Display coefficients with confidence intervals
tidy(m_task1, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Table 1a. Three-Predictor Model Coefficients",
    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)
Table 1a. Three-Predictor Model Coefficients
Term Estimate Std. Error t-statistic p-value 95% CI Lower 95% CI Upper
(Intercept) 10.6684 0.6102 17.4844 0 9.4722 11.8646
physhlth_days 0.3182 0.0131 24.2179 0 0.2924 0.3439
sleep_hrs -0.5063 0.0760 -6.6630 0 -0.6553 -0.3573
age -0.0800 0.0059 -13.4532 0 -0.0916 -0.0683

\[\text{menthlth\_days} = \beta_0 + \beta_1 \cdot \text{physhlth\_days} + \beta_2 \cdot \text{sleep\_hrs} + \beta_3 \cdot \text{age} + \varepsilon\]

#1a
# Extract coefficients for the fitted equation
coefs <- round(coef(m_task1), 4)
cat("Fitted Equation:\n")
## Fitted Equation:
cat(sprintf(
  "menthlth_days = %.4f + %.4f * physhlth_days + %.4f * sleep_hrs + %.4f * age\n",
  coefs["(Intercept)"],
  coefs["physhlth_days"],
  coefs["sleep_hrs"],
  coefs["age"]
))
## menthlth_days = 10.6684 + 0.3182 * physhlth_days + -0.5063 * sleep_hrs + -0.0800 * age

\[\widehat{\text{MentalHealthDays}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \text{physhlth\_days} + \hat{\beta}_2 \cdot \text{sleep\_hrs} + \hat{\beta}_3 \cdot \text{age}\] Interpretation of key coefficients:

physhlth days: For each additional day of physical illness in the past 30 days, mentally unhealthy days increase by approximately 0.32 mentally unhealthy days, holding sleep and age constant. This positive association is expected — physical and mental health are closely linked.

sleep hrs: each additional hour of sleep is associated with a decrease of about 0.51 mentally unhealthy days, holding the other predictors constant, suggesting that more sleep is associated with better mental health.

age: Each additional year of age is associated with a decrease of about 0.08 mentally unhealthy days, adjusting for physical health and sleep.

1b. (5 pts) Use anova() on this model to obtain the Type I (sequential) sums of squares. Present the ANOVA table and verify that the sum of all predictor Type I SS equals the model SSR. Show this calculation explicitly.

# 1b
#Type I (sequential) sums of squares using anova()
anova_task1 <- anova(m_task1)

anova_task1 %>%
  as.data.frame() %>%
  rownames_to_column("Source") %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Table 1b. Type I (Sequential) ANOVA Table — anova(m_task1)",
    col.names = c("Source", "df", "Sum of Sq", "Mean Sq", "F value", "p-value")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  footnote(general = "Variables tested sequentially in the order entered: physhlth_days first, then sleep_hrs, then age.")
Table 1b. Type I (Sequential) ANOVA Table — anova(m_task1)
Source df Sum of Sq Mean Sq F value p-value
physhlth_days 1 29474.746 29474.7462 576.0009 0
sleep_hrs 1 3323.040 3323.0401 64.9395 0
age 1 9261.475 9261.4748 180.9894 0
Residuals 4996 255652.080 51.1714 NA NA
Note:
Variables tested sequentially in the order entered: physhlth_days first, then sleep_hrs, then age. 
#1b
# Verify: predictor Type I SS sum = model SSR
ss_phys  <- anova_task1$`Sum Sq`[1]   # SS(physhlth_days)
ss_sleep <- anova_task1$`Sum Sq`[2]   # SS(sleep_hrs | physhlth_days)
ss_age   <- anova_task1$`Sum Sq`[3]   # SS(age | physhlth_days, sleep_hrs)
sse      <- anova_task1$`Sum Sq`[4]   # Residual SS
ssr      <- ss_phys + ss_sleep + ss_age
ssy      <- ssr + sse

cat("--- Verification: Type I SS Sum = Model SSR ---\n\n")
## --- Verification: Type I SS Sum = Model SSR ---
cat("SS(physhlth_days):                          ", round(ss_phys, 2), "\n")
## SS(physhlth_days):                           29474.75
cat("SS(sleep_hrs | physhlth_days):              ", round(ss_sleep, 2), "\n")
## SS(sleep_hrs | physhlth_days):               3323.04
cat("SS(age | physhlth_days, sleep_hrs):         ", round(ss_age, 2), "\n")
## SS(age | physhlth_days, sleep_hrs):          9261.47
cat("---------------------------------------------------\n")
## ---------------------------------------------------
cat("SSR (sum of predictor Type I SS):           ", round(ssr, 2), "\n")
## SSR (sum of predictor Type I SS):            42059.26
cat("SSE (residual):                             ", round(sse, 2), "\n")
## SSE (residual):                              255652.1
cat("SSY (total = SSR + SSE):                    ", round(ssy, 2), "\n")
## SSY (total = SSR + SSE):                     297711.3
cat("R² = SSR / SSY:                             ", round(ssr / ssy, 4), "\n")
## R² = SSR / SSY:                              0.1413
cat("R² from model summary:                      ", round(summary(m_task1)$r.squared, 4), "\n")
## R² from model summary:                       0.1413
cat("Match:", isTRUE(all.equal(round(ssr/ssy, 4), round(summary(m_task1)$r.squared, 4))), "\n")
## Match: TRUE

Verification: The sum of all three predictor Type I sums of squares equals the total model SSR. This is a defining property of Type I SS — they partition the explained variance sequentially, and therefore must add up to the total explained variation. The resulting \(R^2\) matches the value from summary().

1c. (5 pts) Use car::Anova() with type = "III" on the same model to obtain the Type III sums of squares. Compare the Type I and Type III SS for each variable. Which variable’s SS is the same in both tables? Why?

# Type III sums of squares using car::Anova()
anova3_task1 <- Anova(m_task1, type = "III")

anova3_task1 %>%
  as.data.frame() %>%
  rownames_to_column("Source") %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Table 1c. Type III (Partial) ANOVA Table — car::Anova(type = 'III')",
    col.names = c("Source", "Sum of Sq", "df", "F value", "p-value")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  footnote(general = "Each variable is tested given ALL other variables in the model.")
Table 1c. Type III (Partial) ANOVA Table — car::Anova(type = ‘III’)
Source Sum of Sq df F value p-value
(Intercept) 15643.372 1 305.7057 0
physhlth_days 30012.259 1 586.5051 0
sleep_hrs 2271.809 1 44.3961 0
age 9261.475 1 180.9894 0
Residuals 255652.080 4996 NA NA
Note:
Each variable is tested given ALL other variables in the model. 
# Side-by-side comparison of Type I and Type III SS
comparison_1c <- tibble(
  Variable    = c("physhlth_days", "sleep_hrs", "age"),
  `Type I SS` = round(anova_task1$`Sum Sq`[1:3], 2),
  `Type III SS` = round(anova3_task1$`Sum Sq`[2:4], 2),
  `Difference` = round(anova_task1$`Sum Sq`[1:3] - anova3_task1$`Sum Sq`[2:4], 2),
  `Same?` = ifelse(
    abs(anova_task1$`Sum Sq`[1:3] - anova3_task1$`Sum Sq`[2:4]) < 0.01,
    "Yes ✓", "No ✗"
  )
)

comparison_1c %>%
  kable(caption = "Table 1c-2. Comparison of Type I vs. Type III Sums of Squares") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Table 1c-2. Comparison of Type I vs. Type III Sums of Squares
Variable Type I SS Type III SS Difference Same?
physhlth_days 29474.75 30012.26 -537.51 No ✗
sleep_hrs 3323.04 2271.81 1051.23 No ✗
age 9261.47 9261.47 0.00 Yes ✓ 
cat("\nType I SS for 'age' (last variable entered): ",
    round(anova_task1$`Sum Sq`[3], 4), "\n")
## 
## Type I SS for 'age' (last variable entered):  9261.475
cat("Type III SS for 'age':                        ",
    round(anova3_task1$`Sum Sq`[4], 4), "\n")
## Type III SS for 'age':                         9261.475
cat("Are they equal?",
    isTRUE(all.equal(anova_task1$`Sum Sq`[3], anova3_task1$`Sum Sq`[4],
                     tolerance = 0.01)), "\n")
## Are they equal? TRUE

Intepretation: the last variable entered, age, has identical Type I and Type III sums of squares (9,261.47). This occurs because Type I sums of squares evaluate each variable sequentially, conditioning on the variables that were entered earlier. By the time age is added to the model, both physhlth_days and sleep_hrs have already been accounted for. Since Type III sums of squares also adjust for all other variables in the model, the conditioning is the same for the last variable. As a result, the Type I and Type III sums of squares for age are identical.

Overall, this comparison demonstrates that Type I sums of squares depend on the order of variable entry, whereas Type III sums of squares evaluate each predictor after adjusting for all other variables in the model.

Task 2: Type I vs. Type III Sums of Squares (15 points)

2a. (5 pts) Fit the same model from Task 1 but reverse the variable order:

#2a
# Fit same model with REVERSED variable order
m_task1_rev <- lm(menthlth_days ~ age + sleep_hrs + physhlth_days,
                  data = brfss_mlr)

# Type I SS for reversed model
anova_rev <- anova(m_task1_rev)

anova_rev %>%
  as.data.frame() %>%
  rownames_to_column("Source") %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Table 2a. Type I ANOVA — Reversed Order (age → sleep_hrs → physhlth_days)",
    col.names = c("Source", "df", "Sum of Sq", "Mean Sq", "F value", "p-value")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Table 2a. Type I ANOVA — Reversed Order (age → sleep_hrs → physhlth_days)
Source df Sum of Sq Mean Sq F value p-value
age 1 7249.097 7249.0969 141.6632 0
sleep_hrs 1 4797.905 4797.9051 93.7615 0
physhlth_days 1 30012.259 30012.2591 586.5051 0
Residuals 4996 255652.080 51.1714 NA NA 
# Compare Type I SS between original and reversed ordering
comparison_2a <- tibble(
  Variable = c("physhlth_days", "sleep_hrs", "age"),
  `Type I SS (Original Order)`  = round(anova_task1$`Sum Sq`[1:3], 2),
  `Type I SS (Reversed Order)`  = c(
    round(anova_rev$`Sum Sq`[3], 2),   # physhlth_days is now last
    round(anova_rev$`Sum Sq`[2], 2),   # sleep_hrs is now middle
    round(anova_rev$`Sum Sq`[1], 2)    # age is now first
  ),
  Changed = c("Yes ✗", "Yes ✗", "Yes ✗")
)

comparison_2a %>%
  kable(caption = "Table 2a-2. Type I SS: Original vs. Reversed Variable Order") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  footnote(general = "ALL Type I SS values change when variable order changes — confirming order-dependence.")
Table 2a-2. Type I SS: Original vs. Reversed Variable Order
Variable Type I SS (Original Order) Type I SS (Reversed Order) Changed
physhlth_days 29474.75 30012.26 Yes ✗
sleep_hrs 3323.04 4797.91 Yes ✗
age 9261.47 7249.10 Yes ✗
Note:
ALL Type I SS values change when variable order changes — confirming order-dependence.

Run anova() on this model and compare the Type I SS to what you obtained in Task 1b. Which values changed and which stayed the same?

Interpretation: All Type I SS values changed when the variable order was reversed. This demonstrates the fundamental property of Type I SS: they are entirely order-dependent. In the original order, physhlth_days receives a large Type I SS because it is tested alone (unadjusted). In the reversed order, physhlth_days is entered last, so its Type I SS represents its unique contribution after adjusting for age and sleep_hrs — a much smaller value. The only value that would remain the same across any ordering is the last variable’s Type I SS, because in both orderings the last variable is always conditioned on all others. However, since we completely reversed the order, each variable now occupies a different position, so all values differ.

\[\text{menthlth\_days} = \beta_0 + \beta_1 \cdot \text{age} + \beta_2 \cdot \text{sleep\_hrs} + \beta_3 \cdot \text{physhlth\_days} + \varepsilon\]

2b. (5 pts) Run car::Anova(type = "III") on this reordered model. Did the Type III SS change compared to Task 1c? Explain why or why not.

#2b
# Type III SS on reversed model
anova3_rev <- Anova(m_task1_rev, type = "III")

# Compare Type III from original vs. reversed model
comparison_2b <- tibble(
  Variable = c("physhlth_days", "sleep_hrs", "age"),
  `Type III SS (Original Order)` = round(anova3_task1$`Sum Sq`[2:4], 4),
  `Type III SS (Reversed Order)` = c(
    round(anova3_rev$`Sum Sq`[4], 4),
    round(anova3_rev$`Sum Sq`[3], 4),
    round(anova3_rev$`Sum Sq`[2], 4)
  ),
  `Same?` = "Yes ✓"
)

comparison_2b %>%
  kable(caption = "Table 2b. Type III SS: Original vs. Reversed Order — Order Invariant") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Table 2b. Type III SS: Original vs. Reversed Order — Order Invariant
Variable Type III SS (Original Order) Type III SS (Reversed Order) Same?
physhlth_days 30012.259 30012.259 Yes ✓
sleep_hrs 2271.809 2271.809 Yes ✓
age 9261.475 9261.475 Yes ✓ 
cat("Verification — physhlth_days Type III SS:\n")
## Verification — physhlth_days Type III SS:
cat("  Original order:", round(anova3_task1$`Sum Sq`[2], 4), "\n")
##   Original order: 30012.26
cat("  Reversed order:", round(anova3_rev$`Sum Sq`[4], 4), "\n")
##   Reversed order: 30012.26
cat("  Equal?",
    isTRUE(all.equal(anova3_task1$`Sum Sq`[2], anova3_rev$`Sum Sq`[4],
                     tolerance = 0.01)), "\n")
##   Equal? TRUE

Intepretation: The Type III sums of squares did not change when the order of the variables in the model was reversed. The values for physhlth_days (30012.259), sleep_hrs (2271.809), and age (9261.475) remained exactly the same in both models. This occurs because Type III sums of squares evaluate the contribution of each predictor after adjusting for all other variables in the model, regardless of the order in which the variables are entered. Type III sums of squares are order-independent. As a result, reversing the order of the predictors does not change the Type III SS values.

2c. (5 pts) In 2–3 sentences, explain when an epidemiologist should prefer Type III SS over Type I SS. Give a concrete example from public health research where the choice matters. An epidemiologist should prefer Type III sums of squares when the goal is to estimate the independent effect of each predictor while controlling for all other variables in the model, especially when predictors may be correlated. For example, in a study examining the relationship between income, physical activity, BMI, and diabetes prevalence, using Type I SS with income entered first would attribute shared variation with physical activity and BMI to income, inflating its apparent effect. Type III SS instead estimates the contribution of income after controlling for physical activity and BMI, providing the adjusted association that epidemiologic analyses typically require.

Task 3: Partial F-Tests for Individual Variables (20 points)

3a. (10 pts) Conduct a partial F-test to determine whether age adds significantly to the prediction of mental health days, given that physhlth_days and sleep_hrs are already in the model. Do this by:

  1. Fitting a reduced model (without age)
  2. Fitting the full model (with age)
  3. Using anova(reduced, full) to compare them

State the null hypothesis, report the F-statistic and p-value, and state your conclusion at \(\alpha = 0.05\).

#3a
# Reduced model: without age
m_reduced_3a <- lm(menthlth_days ~ physhlth_days + sleep_hrs,
                   data = brfss_mlr)

# Full model: with age
m_full_3a <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age,
                data = brfss_mlr)

# Partial F-test
pf_test_3a <- anova(m_reduced_3a, m_full_3a)

pf_test_3a %>%
  as.data.frame() %>%
  rownames_to_column("Model") %>%
  mutate(
    Model = c("Reduced (physhlth + sleep only)",
              "Full (physhlth + sleep + age)"),
    across(where(is.numeric), ~ round(., 4))
  ) %>%
  kable(
    caption = "Table 3a. Partial F-Test: Does age add to the model?",
    col.names = c("Model", "Res. df", "RSS", "df", "Sum of Sq", "F", "p-value")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Table 3a. Partial F-Test: Does age add to the model?
Model Res. df RSS df Sum of Sq F p-value
Reduced (physhlth + sleep only) 4997 264913.6 NA NA NA NA
Full (physhlth + sleep + age) 4996 255652.1 1 9261.475 180.9894 0 
#3a
f_val  <- pf_test_3a$F[2]
p_val  <- pf_test_3a$`Pr(>F)`[2]
df1    <- 1
df2    <- pf_test_3a$Res.Df[2]
f_crit <- qf(0.95, df1, df2)

cat("F-statistic:", round(f_val, 4), "\n")
## F-statistic: 180.9894
cat("df (numerator):", df1, "\n")
## df (numerator): 1
cat("df (denominator):", df2, "\n")
## df (denominator): 4996
cat("Critical value F(1,", df2, ", 0.95):", round(f_crit, 4), "\n")
## Critical value F(1, 4996 , 0.95): 3.8433
cat("p-value:", format.pval(p_val, digits = 4), "\n")
## p-value: < 2.2e-16
cat("Reject H0?", ifelse(p_val < 0.05, "Yes — age is significant", "No — age is not significant"), "\n")
## Reject H0? Yes — age is significant

Interpretation: The F-statistic is greater than the critical value, and the p-value is below 0.05 (2.2 x 10^(-16)), indicating that we can reject null hypothesis. This shows that age contributes statistically significant predictive information for mentally unhealthy days, above and beyond what is explained by physical health days and sleep hours. Therefore, age should remain in the model.

3b. (10 pts) Now verify your result from 3a manually. Using the anova() output from the full model (Task 1b), identify \(SS(\text{age} \mid \text{physhlth\_days}, \text{sleep\_hrs})\) from the Type I table. Compute the F-statistic as:

\[F = \frac{SS(\text{age} \mid \text{physhlth\_days}, \text{sleep\_hrs}) / 1}{MSE(\text{full model})}\]

# Get Type I ANOVA output for the full 3-predictor model
# (entered as physhlth_days, sleep_hrs, age — so age is last)
anova_full_3b <- anova(m_full_3a)

# Extract needed components
ss_age_typeI <- anova_full_3b$`Sum Sq`[3]   # SS(age | physhlth, sleep) — Type I, entered last
mse_full     <- anova_full_3b$`Mean Sq`[4]  # MSE of full model
df2_manual   <- anova_full_3b$Df[4]         # Residual df

# Compute F-statistic manually
F_manual <- (ss_age_typeI / 1) / mse_full

# Critical value
f_crit_manual <- qf(0.95, 1, df2_manual)

cat("=== Manual Partial F-Test for age ===\n\n")
## === Manual Partial F-Test for age ===
cat("SS(age | physhlth_days, sleep_hrs) from Type I table:", round(ss_age_typeI, 2), "\n")
## SS(age | physhlth_days, sleep_hrs) from Type I table: 9261.47
cat("MSE (full model):                                    ", round(mse_full, 4), "\n")
## MSE (full model):                                     51.1714
cat("Numerator df:                                        ", 1, "\n")
## Numerator df:                                         1
cat("Denominator df (residual):                           ", df2_manual, "\n")
## Denominator df (residual):                            4996
cat("---------------------------------------------------\n")
## ---------------------------------------------------
cat("F = SS(age) / 1 / MSE(full) =",
    round(ss_age_typeI, 2), "/", round(mse_full, 4), "=",
    round(F_manual, 4), "\n")
## F = SS(age) / 1 / MSE(full) = 9261.47 / 51.1714 = 180.9894
cat("Critical value F(1,", df2_manual, ", 0.95):         ", round(f_crit_manual, 4), "\n")
## Critical value F(1, 4996 , 0.95):          3.8433
cat("F > critical value?", ifelse(F_manual > f_crit_manual, "Yes — Reject H0", "No — Fail to reject H0"), "\n")
## F > critical value? Yes — Reject H0
cat("\nComparison with anova() result:\n")
## 
## Comparison with anova() result:
cat("  F from anova(reduced, full):", round(pf_test_3a$F[2], 4), "\n")
##   F from anova(reduced, full): 180.9894
cat("  F from manual calculation:  ", round(F_manual, 4), "\n")
##   F from manual calculation:   180.9894
cat("  Match?",
    isTRUE(all.equal(round(pf_test_3a$F[2], 2), round(F_manual, 2))), "\n")
##   Match? TRUE

Compare to the critical value \(F_{1, n-p-1, 0.95}\). Does your manual calculation agree with the anova() comparison from 3a? Intepretation:Since age is the last variable entered in the model, its Type I sum of squares is identical to its Type III sum of squares (i.e., SS(age∣physhlth_days,sleep_hrs)). The F-statistic calculated manually from the Type I table matches exactly the F-statistic obtained from anova(reduced, full). Both values exceed the critical F value at the 0.95 confidence level and have p<0.05, therefore, we will reject null H0. —

Round 2: Student B Drives, Student A Navigates

⟳ Switch roles now!

Task 4: T-Tests and the F-Test Equivalence (15 points)

4a. (5 pts) Using the full model from Task 1 (menthlth_days ~ physhlth_days + sleep_hrs + age), run summary() and extract the t-statistics and p-values for each coefficient.

#4a
# Full model from Task 1
m_task4 <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age,
              data = brfss_mlr)

# Extract t-statistics and p-values
t_results <- tidy(m_task4) %>%
  filter(term != "(Intercept)") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  mutate(across(where(is.numeric), ~ round(., 4)))

t_results %>%
  kable(
    caption = "Table 4a. T-Statistics and P-Values from summary()",
    col.names = c("Predictor", "Estimate (β̂)", "Std. Error", "t-statistic", "p-value")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Table 4a. T-Statistics and P-Values from summary()
Predictor Estimate (β̂
Std. Erro
physhlth_days 0.3182 0.0131 24.2179 0
sleep_hrs -0.5063 0.0760 -6.6630 0
age -0.0800 0.0059 -13.4532 0

Interpretation: t = -13.45, p < 0.001: All three predictors (physhlth_days, sleep_hrs, age) have p-values less than 0.05, indicating each independently contributes to predicting mentally unhealthy days. This effect is highly significant, indicating that older adults report slightly fewer mentally unhealthy day

4b. (5 pts) For each predictor, compute \(t^2\) and compare it to the Type III F-statistic from Task 1c.

#4b
# Get Type III F-statistics
f_results <- Anova(m_task4, type = "III") %>%
  as.data.frame() %>%
  rownames_to_column("term") %>%
  filter(!term %in% c("(Intercept)", "Residuals")) %>%
  select(term, f_stat = `F value`, f_pvalue = `Pr(>F)`)

# Combine and compute t²
equivalence_table <- t_results %>%
  select(term, t_stat = statistic, t_pvalue = p.value) %>%
  left_join(f_results, by = "term") %>%
  mutate(
    `t²`         = round(t_stat^2, 4),
    `F (Type III)` = round(f_stat, 4),
    `t² = F?`    = ifelse(abs(t_stat^2 - f_stat) < 0.01, "Yes ✓", "No ✗"),
    `p (t-test)` = round(t_pvalue, 6),
    `p (F-test)` = round(f_pvalue, 6),
    `p-values equal?` = ifelse(abs(t_pvalue - f_pvalue) < 0.0001, "Yes ✓", "No ✗")
  ) %>%
  select(term, t_stat, `t²`, `F (Type III)`, `t² = F?`,
         `p (t-test)`, `p (F-test)`, `p-values equal?`) %>%
  mutate(t_stat = round(t_stat, 4))

equivalence_table %>%
  kable(caption = "Table 4b. Equivalence of t-Tests and Type III Partial F-Tests") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  column_spec(5, bold = TRUE, color = "darkgreen") %>%
  column_spec(8, bold = TRUE, color = "darkgreen")
Table 4b. Equivalence of t-Tests and Type III Partial F-Tests
term t_stat t² F (Type III) t² = F? p (t-test) p (F-test) p-values equal?
physhlth_days 24.2179 586.5067 586.5051 Yes ✓ 0 0 Yes ✓
sleep_hrs -6.6630 44.3956 44.3961 Yes ✓ 0 0 Yes ✓
age -13.4532 180.9886 180.9894 Yes ✓ 0 0 Yes ✓

Create a table showing the t-statistic, \(t^2\), the Type III F-statistic, and both p-values. Are they equivalent? Intepretation: The table demonstrates the equivalence of t-tests and Type III partial F-tests in multiple regression: for each individual predictor, t^2=𝐹 and the p-values are identical. This confirms that both approaches give the same inference regarding the statistical significance of each variable.

4c. (5 pts) In your own words, explain why the t-test and the Type III partial F-test give the same result. What is the fundamental relationship between the t-distribution and the F-distribution that makes this true? The t-test is significantly different from zero, while the Type III partial F-test tests the same hypothesis by comparing the variance explained by a model with and without that predictor, adjusting for all other variables. The reason they give the same result for a single predictor is that the F-distribution is the square of the t-distribution when the numerator degrees of freedom is 1 as F1,df=tdf2. This fundamental relationship means that whether you compute a t-test for a single coefficient or a Type III F-test for the same variable, the statistic and the p-value are identical. —

Task 5: Chunk Test — Testing Groups of Variables (20 points)

5a. (10 pts) Now consider the full 6-predictor model:

\[\text{menthlth\_days} = \beta_0 + \beta_1 \cdot \text{physhlth\_days} + \beta_2 \cdot \text{sleep\_hrs} + \beta_3 \cdot \text{age} + \beta_4 \cdot \text{income\_cat} + \beta_5 \cdot \text{sex} + \beta_6 \cdot \text{exercise} + \varepsilon\]

Test whether income_cat, sex, and exercise — as a group — significantly add to the prediction of mental health days, given that physhlth_days, sleep_hrs, and age are already in the model.

# Reduced model: physhlth_days + sleep_hrs + age only
m_reduced_5a <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age,
                   data = brfss_mlr)

# Full model: all six predictors
m_full_5a <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age +
                  income_cat + sex + exercise,
                data = brfss_mlr)

# Chunk test
chunk_test_5a <- anova(m_reduced_5a, m_full_5a)

chunk_test_5a %>%
  as.data.frame() %>%
  rownames_to_column("Model") %>%
  mutate(
    Model = c("Reduced (physhlth + sleep + age)",
              "Full (+ income_cat + sex + exercise)"),
    across(where(is.numeric), ~ round(., 4))
  ) %>%
  kable(
    caption = "Table 5a. Chunk Test: Do income_cat, sex, and exercise collectively add to the model?",
    col.names = c("Model", "Res. df", "RSS", "df", "Sum of Sq", "F", "p-value")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Table 5a. Chunk Test: Do income_cat, sex, and exercise collectively add to the model?
Model Res. df RSS df Sum of Sq F p-value
Reduced (physhlth + sleep + age) 4996 255652.1 NA NA NA NA
Full (+ income_cat + sex + exercise) 4993 250992.8 3 4659.276 30.8957 0 
f_chunk <- chunk_test_5a$F[2]
p_chunk <- chunk_test_5a$`Pr(>F)`[2]
df1_chunk <- chunk_test_5a$Df[2]
df2_chunk <- chunk_test_5a$Res.Df[2]
f_crit_chunk <- qf(0.95, df1_chunk, df2_chunk)

cat("F-statistic:          ", round(f_chunk, 4), "\n")
## F-statistic:           30.8957
cat("df (numerator):       ", df1_chunk, "(3 additional variables)\n")
## df (numerator):        3 (3 additional variables)
cat("df (denominator):     ", df2_chunk, "\n")
## df (denominator):      4993
cat("Critical value F(3,", df2_chunk, ", 0.95):", round(f_crit_chunk, 4), "\n")
## Critical value F(3, 4993 , 0.95): 2.6067
cat("p-value:              ", format.pval(p_chunk, digits = 4), "\n")
## p-value:               < 2.2e-16
cat("Reject H0?",
    ifelse(p_chunk < 0.05,
           "Yes — the group of variables is collectively significant",
           "No — the group does not add significantly"), "\n")
## Reject H0? Yes — the group of variables is collectively significant

State the null hypothesis (in both words and mathematical notation), conduct the test, and state your conclusion. The null hypothesis: the group of predictors—physhlth_days, sleep_hrs, and age—does not collectively explain any variation in mentally unhealthy days beyond what is expected by chance.

Mathematical notation: 𝐻0:𝛽physhlth_days=𝛽sleep_hrs=𝛽age=0

5b. (5 pts) Compute the chunk test F-statistic manually using:

\[F = \frac{\{SSR(\text{full}) - SSR(\text{reduced})\} / (df_{\text{full}} - df_{\text{reduced}})}{MSE(\text{full})}\]

Show all intermediate values.

# Compute SSR for full and reduced models
anova_full_5b    <- anova(m_full_5a)
anova_reduced_5b <- anova(m_reduced_5a)

ssr_full_5b    <- sum(anova_full_5b$`Sum Sq`[1:6])   # 6 predictors
ssr_reduced_5b <- sum(anova_reduced_5b$`Sum Sq`[1:3]) # 3 predictors
mse_full_5b    <- anova_full_5b$`Mean Sq`[7]          # MSE of full model
df_diff_5b     <- 3                                    # 3 additional variables added

# Manual F-statistic
F_manual_chunk <- ((ssr_full_5b - ssr_reduced_5b) / df_diff_5b) / mse_full_5b

cat("=== Manual Chunk F-Test ===\n\n")
## === Manual Chunk F-Test ===
cat("SSR (full model, 6 predictors):    ", round(ssr_full_5b, 2), "\n")
## SSR (full model, 6 predictors):     46718.54
cat("SSR (reduced model, 3 predictors): ", round(ssr_reduced_5b, 2), "\n")
## SSR (reduced model, 3 predictors):  42059.26
cat("Difference (SSR_full - SSR_red):   ", round(ssr_full_5b - ssr_reduced_5b, 2), "\n")
## Difference (SSR_full - SSR_red):    4659.28
cat("Number of additional variables:    ", df_diff_5b, "\n")
## Number of additional variables:     3
cat("MSE (full model):                  ", round(mse_full_5b, 4), "\n")
## MSE (full model):                   50.2689
cat("--------------------------------------------------\n")
## --------------------------------------------------
cat("F = [(SSR_full - SSR_red) / df_diff] / MSE(full)\n")
## F = [(SSR_full - SSR_red) / df_diff] / MSE(full)
cat("F = [", round(ssr_full_5b - ssr_reduced_5b, 2),
    "/", df_diff_5b, "] /", round(mse_full_5b, 4),
    "=", round(F_manual_chunk, 4), "\n\n")
## F = [ 4659.28 / 3 ] / 50.2689 = 30.8957
cat("F from anova(reduced, full):", round(chunk_test_5a$F[2], 4), "\n")
## F from anova(reduced, full): 30.8957
cat("F from manual calculation:  ", round(F_manual_chunk, 4), "\n")
## F from manual calculation:   30.8957
cat("Match?",
    isTRUE(all.equal(round(chunk_test_5a$F[2], 2),
                     round(F_manual_chunk, 2))), "\n")
## Match? TRUE

Does your manual computation match the anova() result? The manual calculation confirms the ANOVA function results, showing that the three additional predictors collectively explain a statistically significant portion of the variation in mentally unhealthy days.The formula partitions the additional variation explained by the three new variables (numerator) relative to the unexplained variation in the full model (denominator). The numerator degrees of freedom correspond to the number of constraints imposed under 𝐻0 in this case, 3, one for each variable in the group.

5c. (5 pts) Note that exercise was not individually significant in the Type III table, yet it is part of a group that is collectively significant. In 2–3 sentences, explain how this is possible and what it means for model building in epidemiology. This occurs because Type III tests and chunk (block) F-tests address different questions: the Type III test evaluates whether a single variable adds unique predictive value after adjusting for all others, while the chunk test assesses whether a group of variables collectively explains additional variation beyond a reduced model. In epidemiologic model-building, this means that individual non-significance does not automatically justify dropping a variable, especially if it is theoretically or clinically important. Decisions about variable retention should be guided by subject-matter knowledge and confounding structure, not solely by p-values. —

Task 6: Synthesis and Interpretation (15 points)

6a. (5 pts) Based on the full model, which predictors are statistically significant at \(\alpha = 0.05\)? List them and briefly state the direction of each association (positive or negative).

# Full six-predictor model Type III results
type3_full_6 <- Anova(m_full_5a, type = "III") %>%
  as.data.frame() %>%
  rownames_to_column("Predictor") %>%
  filter(!Predictor %in% c("(Intercept)", "Residuals")) %>%
  mutate(
    `Significant (α=0.05)` = ifelse(`Pr(>F)` < 0.05, "Yes ✓", "No ✗"),
    across(where(is.numeric), ~ round(., 4))
  )

# Get coefficient signs from the model
coef_signs <- tidy(m_full_5a) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    Direction = case_when(
      estimate > 0 ~ "Positive (+)",
      estimate < 0 ~ "Negative (−)",
      TRUE ~ "Zero"
    )
  ) %>%
  select(term, Direction)

# Build final table
final_table_6a <- type3_full_6 %>%
  left_join(coef_signs, by = c("Predictor" = "term")) %>%
  select(Predictor, `Sum Sq`, `Df`, `F value`, `Pr(>F)`,
         `Significant (α=0.05)`, Direction)

final_table_6a %>%
  kable(
    caption = "Table 6a. Type III Tests — Full Six-Predictor Model",
    col.names = c("Predictor", "Sum of Sq", "df", "F value",
                  "p-value", "Significant?", "Direction of Association")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  column_spec(6, bold = TRUE)
Table 6a. Type III Tests — Full Six-Predictor Model
Predictor Sum of Sq df F value p-value Significant? Direction of Association
physhlth_days 23189.1197 1 461.3012 0.0000 Yes ✓ Positive (+)
sleep_hrs 2295.4257 1 45.6629 0.0000 Yes ✓ Negative (−)
age 9673.9387 1 192.4437 0.0000 Yes ✓ Negative (−)
income_cat 1918.5269 1 38.1653 0.0000 Yes ✓ Negative (−)
sex 1903.4532 1 37.8654 0.0000 Yes ✓ NA
exercise 92.1241 1 1.8326 0.1759 No ✗ NA

Physical health days (physhlth_days) Sum of Squares = 23,189.12, F = 461.30, p < 0.001, Statistically significant. Direction: Positive (+): more physically unhealthy days are associated with more mentally unhealthy days.

Sleep hours (sleep_hrs) Sum of Squares = 2,295.43, F = 45.66, p < 0.001, Statistically significant. Direction: Negative (−) — more sleep is associated with fewer mentally unhealthy days.

Age (age) Sum of Squares = 9,673.94, F = 192.44, p < 0.001, Statistically significant. Direction: Negative (−) — older individuals tend to report fewer mentally unhealthy days.

Income category (income_cat) Sum of Squares = 1,918.53, F = 38.17, p < 0.001, Statistically significant. Direction: Negative (−) — higher income is associated with fewer mentally unhealthy days.

Sex (sex) Sum of Squares = 1,903.45, F = 37.87, p < 0.001, Statistically significant. Direction: NA — significant difference exists, but “direction” is not meaningful for categorical variables coded as factor levels.

Exercise (exercise) Sum of Squares = 92.12, F = 1.83, p = 0.176, Not statistically significant. Direction: NA — no evidence of association in this model.

6b. (5 pts) A colleague argues: “We should drop exercise from the model because it’s not significant.” Do you agree? Write a 2–3 sentence response explaining your reasoning. Consider the chunk test results and epidemiologic rationale. I would not prefer to drop exercise from the model. Although exercise is not statistically significant on its own (p = 0.176), the chunk test results may show that it contributes meaningfully as part of a group of variables, and epidemiologic knowledge indicates that exercise is an established determinant of mental health. Removing it solely based on the p-value could introduce bias, so variable retention should consider theoretical importance and confounding structure, not just statistical significance.

6c. (5 pts) Write a 3–4 sentence summary of the hypothesis testing results for a non-statistical audience (e.g., a public health program manager). Your summary should convey which factors were identified as independently associated with mental health days and which were not, without using jargon like “p-value,” “F-test,” or “sums of squares.” Our analysis examined which individual and lifestyle factors were independently linked to the number of mentally unhealthy days reported by U.S. adults. After accounting for all factors simultaneously, we found that physical health, sleep, age, income, and sex each had meaningful associations with mental health days. Adults who experienced more physically unhealthy days, slept less, had lower income, or belonged to a particular sex tended to report more mentally unhealthy days, even when all other factors were considered. Regular exercise was not independently associated with mental health days once the other factors were included, although it contributed as part of a broader set of related factors. These findings suggest that interventions focusing on improving sleep, supporting physical health, and addressing income disparities may be especially effective in promoting better mental health at the population level. ### Task 7: Reflection (Not Graded, strongly recommended) In 2–3 sentences, reflect on the process of working with a partner in this lab. What did you find helpful about the Driver-Navigator approach? How did it affect your understanding of the material? Working with a partner using the Driver-Navigator approach was very helpful, it encouraged active discussion, clarified points of confusion in real time, and helped catch mistakes early. —

End of Lab Activity