knitr::opts_chunk$set(
  echo    = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width  = 10,
  fig.height = 6,
  cache  = FALSE
)

Lab Activity #8

Data for the Lab

Use the same BRFSS 2020 dataset from the guided practice.

# Load packages and data
library(tidyverse)
library(haven)
library(janitor)
library(knitr)
library(kableExtra)
library(plotly)
library(broom)
library(ggeffects)
library(ggstats)
library(gtsummary)
library(GGally)
library(car)
library(lmtest)
library(corrplot)
library(dplyr)

brfss_mlr <- readRDS(
  "C:/Users/graci/OneDrive/Documents/UA GRAD SCHOOL/2nd Semester/EPI553/BRFSS_mlr_2020.rds"
)
brfss_mlr %>%
  select(menthlth_days, physhlth_days, sleep_hrs, age, income_cat, sex, exercise) %>%
  tbl_summary(
    statistic = list(
      all_continuous()  ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    label = list(
      menthlth_days ~ "Mentally Unhealthy Days (past 30)",
      physhlth_days ~ "Physically Unhealthy Days (past 30)",
      sleep_hrs     ~ "Sleep Hours per Night",
      age           ~ "Age (years)",
      income_cat    ~ "Income Category (1–8)",
      sex           ~ "Sex",
      exercise      ~ "Any Exercise (past 30 days)"
    )
  ) %>%
  bold_labels()
Characteristic N = 5,0001
Mentally Unhealthy Days (past 30) 4 (8)
Physically Unhealthy Days (past 30) 3 (8)
Sleep Hours per Night 7.06 (1.35)
Age (years) 54 (17)
Income Category (1–8)
    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%)
Sex
    Male 2,331 (47%)
    Female 2,669 (53%)
Any Exercise (past 30 days) 3,874 (77%)
1 Mean (SD); n (%)

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

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

\[\text{menthlth\_days} = \beta_0 + \beta_1 \cdot \text{physhlth\_days} + \beta_2 \cdot \text{sleep\_hrs} + \beta_3 \cdot \text{age} + \varepsilon\] Menthlth_days= 10.6684+ 0.3182(physhlth_days)-0.5063(sleep_hrs)-0.0800(age)+e

Use tidy() with conf.int = TRUE to display the coefficients. Report the fitted equation with rounded coefficients.

m_full <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age,
             data = brfss_mlr)

tidy(m_full, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Table 1. Full 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 1. Full 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

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.

SSR=SS(physhlth_days)+SS(sleep_hrs)+SS(age) SSR=29474.75+3323.04+9261.47 SSR= 42,059.26

anova_type1 <- anova(m_full)

anova_type1 %>%
  as.data.frame() %>%
  rownames_to_column("Source") %>%
  mutate(across(where(is.numeric), ~ round(., 2))) %>%
  kable(
    caption = "Table 3. Type I (Sequential) Sums of Squares — anova()",
    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 are tested in the order they appear in the model formula.")
Table 3. Type I (Sequential) Sums of Squares — anova()
Source df Sum of Sq Mean Sq F value p-value
physhlth_days 1 29474.75 29474.75 576.00 0
sleep_hrs 1 3323.04 3323.04 64.94 0
age 1 9261.47 9261.47 180.99 0
Residuals 4996 255652.08 51.17 NA NA
Note:
Variables are tested in the order they appear in the model formula.

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?

The variable with the same SS is age. This happens because Type I SS depend on the order in which variables are entered. Changing the order changes the Type I SS for all but the last variable.

anova_type3 <- Anova(m_full, type = "III")

# Side-by-side comparison
comparison <- tibble(
  Variable = c("physhlth_days", "sleep_hrs", "age"),
  `Type I SS` = round(anova_type1$`Sum Sq`[1:3], 1),
  `Type III SS` = round(anova_type3$`Sum Sq`[2:4], 1)
)

comparison %>%
  kable(caption = "Table 4. Type I vs. Type III Sums of Squares") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  footnote(general = "Type I SS depends on variable entry order; Type III SS does not. The last variable (exercise) has identical values.")
Table 4. Type I vs. Type III Sums of Squares
Variable Type I SS Type III SS
physhlth_days 29474.7 30012.3
sleep_hrs 3323.0 2271.8
age 9261.5 9261.5
Note:
Type I SS depends on variable entry order; Type III SS does not. The last variable (exercise) has identical values.

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:

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

m_full <- lm(menthlth_days ~ age + sleep_hrs + physhlth_days,
             data = brfss_mlr)

tidy(m_full, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(
    caption = "Table 1. Full 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 1. Full 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
age -0.0800 0.0059 -13.4532 0 -0.0916 -0.0683
sleep_hrs -0.5063 0.0760 -6.6630 0 -0.6553 -0.3573
physhlth_days 0.3182 0.0131 24.2179 0 0.2924 0.3439

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?

The values that changed was age in type I. Physhlth_days and sleep_hrs stayed the same for both type I and type II SS.

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.

The type III SS did not change when compared to task 1c. This is because type III SS do not depend on the order of variable entry. They are the same regardless of how you specify the model formula.

anova_type3 <- Anova(m_full, type = "III")

# Side-by-side comparison
comparison <- tibble(
  Variable = c("age", "sleep_hrs", "physhlth_days"),
  `Type I SS` = round(anova_type1$`Sum Sq`[1:3], 1),
  `Type III SS` = round(anova_type3$`Sum Sq`[2:4], 1)
)

comparison %>%
  kable(caption = "Table 4. Type I vs. Type III Sums of Squares") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  footnote(general = "Type I SS depends on variable entry order; Type III SS does not. The last variable (exercise) has identical values.")
Table 4. Type I vs. Type III Sums of Squares
Variable Type I SS Type III SS
age 29474.7 9261.5
sleep_hrs 3323.0 2271.8
physhlth_days 9261.5 30012.3
Note:
Type I SS depends on variable entry order; Type III SS does not. The last variable (exercise) has identical 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.

Epidemiologist should prefer Type III when we want to know the effect of each variable after adjusting for all others. For example, if you are examining confounding factors like age, income, or physical activity on depression, these predictors may be correlated.Using type III SS allows researchers to determine if the confounding factors are contributing to mental health outcomes.

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
m_no_age <- lm(menthlth_days ~ physhlth_days + sleep_hrs,
                 data = brfss_mlr)

# Full model (includes age)
m_full <- lm(menthlth_days ~ physhlth_days + age + sleep_hrs,
             data = brfss_mlr)

# Compare using anova() — this performs the partial F-test
anova(m_no_age, m_full) %>%
  as.data.frame() %>%
  rownames_to_column("Model") %>%
  mutate(
    Model = c("Reduced (no age)", "Full (with age)"),
    across(where(is.numeric), ~ round(., 4))
  ) %>%
  kable(
    caption = "Table 6. 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 6. Partial F-Test: Does age add to the model?
Model Res. df RSS df Sum of Sq F p-value
Reduced (no age) 4997 264913.6 NA NA NA NA
Full (with age) 4996 255652.1 1 9261.475 180.9894 0

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

Null Hypothesis: Age does not improve prediction of mental health days after accounting for physical health days and sleep. H0:βage=0 F-Statistic: 180.9894 P-value: 0

Conclusion: The p-value was less than 0.05, we can reject the null hypothesis. Age significantly improves the prediction of mental health days after controlling for physical health days and sleep.

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:

m_no_age <- lm(menthlth_days ~ physhlth_days + sleep_hrs,
                 data = brfss_mlr)

# Full model (includes age)
m_full <- lm(menthlth_days ~ physhlth_days + age + sleep_hrs,
             data = brfss_mlr)

anova( m_no_age, m_full)
## Analysis of Variance Table
## 
## Model 1: menthlth_days ~ physhlth_days + sleep_hrs
## Model 2: menthlth_days ~ physhlth_days + age + sleep_hrs
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   4997 264914                                  
## 2   4996 255652  1    9261.5 180.99 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ssr_full    <- sum(anova(m_full)$`Sum Sq`[1:3])
ssr_reduced <- sum(anova(m_no_age)$`Sum Sq`[1:2])
mse_full    <- anova(m_full)$`Mean Sq`[4]  

F_stat <- (ssr_full - ssr_reduced) / 1 / mse_full

cat("SSR(full):", round(ssr_full, 2), "\n")
## SSR(full): 42059.26
cat("SSR(reduced):", round(ssr_reduced, 2), "\n")
## SSR(reduced): 32797.79
cat("SS(age | others):", round(ssr_full - ssr_reduced, 2), "\n")
## SS(age | others): 9261.47
cat("MSE(full):", round(mse_full, 2), "\n")
## MSE(full): 51.17
cat("F-statistic:", round(F_stat, 4), "\n")
## F-statistic: 180.9894
cat("Critical value F(1, 4996, 0.95):", round(qf(0.95, 1, 4996), 4), "\n")
## Critical value F(1, 4996, 0.95): 3.8433
cat("p-value:", format.pval(pf(F_stat, 1, 4996, lower.tail = FALSE)),"\n")
## p-value: < 2.22e-16

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

Compare to the critical value \(F_{1, n-p-1, 0.95}\). Does your manual calculation agree with the anova() comparison from 3a?

SSfull=29474.75+3323.04+9261.47=42059.26

SSreduced=29474.75+3323.04=32800.79

42059.26−32800.79=9258.47/51.17= 180.99

Yes the manual calculation agrees with the anova comparison from 3a and are both 180.99.

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.

t_stats <- tidy(m_full) %>%
  filter(term != "(Intercept)") %>%
  select(term, t_stat = statistic, t_pvalue = p.value)

# Get F-statistics from Type III
f_stats <- Anova(m_full, 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)`)

# Compare
left_join(t_stats, f_stats, by = "term") %>%
  mutate(
    `t²` = round(t_stat^2, 4),
    f_stat = round(f_stat, 4),
    `t² = F?` = ifelse(abs(t_stat^2 - f_stat) < 0.001, "✓", "✗"),
    t_pvalue = round(t_pvalue, 6),
    f_pvalue = round(f_pvalue, 6),
    `p-values equal?` = ifelse(abs(t_pvalue - f_pvalue) < 0.0001, "✓", "✗")
  ) %>%
  select(term, t_stat = t_stat, `t²`, `F (Type III)` = f_stat, `t² = F?`,
         `p (t-test)` = t_pvalue, `p (F-test)` = f_pvalue, `p-values equal?`) %>%
  mutate(t_stat = round(t_stat, 4)) %>%
  kable(caption = "Table 9. Equivalence of T-Tests and Type III Partial F-Tests") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Table 9. 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.5051 586.5051 ✓ 0 0 ✓
age -13.4532 180.9894 180.9894 ✓ 0 0 ✓
sleep_hrs -6.6630 44.3961 44.3961 ✓ 0 0 ✓

4b. (5 pts) For each predictor, compute \(t^2\) and compare it to the Type III F-statistic from Task 1c. Create a table showing the t-statistic, \(t^2\), the Type III F-statistic, and both p-values. Are they equivalent?

Type III F-Statistic from task 1c: age= 9261.5 sleep_hours=2271.8 Physhlth_days= 30012.3

T^2 for table 9: T^2 age=180.9894 T^2 sleep_hrs=44.3961 T^2 physhlth_days= 586.5051

Yes the values are equivalent. The square of the t-statistic from the coefficient test is exactly the value of the type III F-statistic.

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 and type III F-test give the same results because both tests are evaluating whether the regression coefficient equals zero. This results in both tests producing identical values and conclusions. In addition when, the numerator degrees of freedom equal 1, the F-test is equal to the square of the T-test.

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.

m_health_only <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age, data = brfss_mlr)

# Full model: health behaviors + demographics
m_full <- lm(menthlth_days ~ physhlth_days + sleep_hrs + age + income_cat + sex + exercise,
             data = brfss_mlr)

# Chunk test
anova(m_health_only, m_full) %>%
  as.data.frame() %>%
  rownames_to_column("Model") %>%
  mutate(
    Model = c("Reduced (physhlth + sleep +age)", "Full (+ demographics)"),
    across(where(is.numeric), ~ round(., 4))
  ) %>%
  kable(
    caption = "Table 10. Chunk Test: Do demographic variables 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 10. Chunk Test: Do demographic variables 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 (+ demographics) 4993 250992.8 3 4659.276 30.8957 0

State the null hypothesis (in both words and mathematical notation), conduct the test, and state your conclusion.

Null hypothesis: Income, sex, and exercise do not improve prediction of mental health days after accounting for physical health, sleep and age. H0:βincome=βsex=βexercise=0

Conclusion: The ANOVA model shows a significant F-statistic with a very small p-value. We can reject the null hypothesis and conclude that income, sex, and exercise improve mental health days.

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. Does your manual computation match the anova() result?

The manually computed F-statistic matched the value obtained in ANOVA and confirms the chunk test.

ssr_full    <- sum(anova(m_full)$`Sum Sq`[1:6])
ssr_reduced <- sum(anova(m_health_only)$`Sum Sq`[1:3])
mse_full    <- anova(m_full)$`Mean Sq`[7]
df_diff     <- 4  # 4 additional variables

F_chunk <- ((ssr_full - ssr_reduced) / df_diff) / mse_full

cat("SSR(full):", round(ssr_full, 2), "\n")
## SSR(full): 46718.54
cat("SSR(reduced):", round(ssr_reduced, 2), "\n")
## SSR(reduced): 42059.26
cat("Difference:", round(ssr_full - ssr_reduced, 2), "\n")
## Difference: 4659.28
cat("df (number of added variables):", df_diff, "\n")
## df (number of added variables): 4
cat("MSE(full):", round(mse_full, 2), "\n")
## MSE(full): 50.27
cat("F-statistic:", round(F_chunk, 4), "\n")
## F-statistic: 23.1717
cat("Critical value F(4, 4993, 0.95):", round(qf(0.95, df_diff, 4993), 4), "\n")
## Critical value F(4, 4993, 0.95): 2.3737
cat("p-value:", format.pval(pf(F_chunk, df_diff, 4993, lower.tail = FALSE)), "\n")
## p-value: < 2.22e-16

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.

A variable can not be significant on its own but when several related variables are tested together, the combined power is substantial enough to make the chunk test significant. This means that even if the single variable appears weak on its own, it may still be relevant in a broader context.

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).

Based on the full model, physhlth_days, sleep_hrs, and age are all statistically significant because all the p-values were less than 0.05. Physhlth_days had a positive association. Sleep_hrs and age had a negative association.

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 agree with my colleague because although it is not individually significant, the chunk test showed that the variables collectively improved the model including exercise. Variables may be retained for theoretical or confounding control purposes, even if their individual statisitcal signifcance is weak.

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.”

The results showed that physical health, sleep, and age are all important factors that are associated with the number of poor mental health days. Individuals with more physical health problems tended to report more poor mental health days. Individuals who slept more reported less poor mental health days. When sex, exercise, and income were looked at together also explained mental health differences. The results suggest that health behaviors and demographic factors play a role in mental health outcomes.