EPI 553 Progress Check-in 2: Regression Results and Visualizations
Fizza Zaheer
April 7, 2026
Section 1: Analytical Sample Update
Step 1.1: Rebuild Analytical Sample from Check-in 1
The analytical sample is unchanged from Check-in 1.
The same six NHANES 2017–2018 files are loaded, merged by
SEQN, and the same four exclusion criteria are applied in
the same sequence. The final analytical sample remains N = —
participants (populated on knit).
## Downloading NHANES 2017-2018 files...demo <- nhanes("DEMO_J")
slq <- nhanes("SLQ_J")
paq <- nhanes("PAQ_J")
bmx <- nhanes("BMX_J")
dpq <- nhanes("DPQ_J")
inq <- nhanes("INQ_J")
cat("Files loaded:\n")## Files loaded:cat(" DEMO_J:", nrow(demo), "| SLQ_J:", nrow(slq),
"| PAQ_J:", nrow(paq), "| BMX_J:", nrow(bmx),
"| DPQ_J:", nrow(dpq), "| INQ_J:", nrow(inq), "\n")## DEMO_J: 9254 | SLQ_J: 6161 | PAQ_J: 5856 | BMX_J: 8704 | DPQ_J: 5533 | INQ_J: 9254# Merge all files by SEQN
nhanes_merged <- demo |>
left_join(slq, by = "SEQN") |>
left_join(paq, by = "SEQN") |>
left_join(bmx, by = "SEQN") |>
left_join(dpq, by = "SEQN") |>
left_join(inq, by = "SEQN")
starting_n <- nrow(nhanes_merged)
# Apply exclusion criteria — identical to Check-in 1
step1 <- nhanes_merged |> dplyr::filter(RIDAGEYR >= 20 & RIDAGEYR <= 70)
n1 <- nrow(step1)
step2 <- step1 |> dplyr::filter(!is.na(SLD012))
n2 <- nrow(step2)
# PAQ605: strip haven labels with as.integer() — 1=Yes, 2=No, 7/9=refused/DK
# as.integer() always returns bare codes regardless of label state
cat("PAQ605 raw integer values (unique):", paste(sort(unique(as.integer(step2$PAQ605))), collapse=", "), "\n")## PAQ605 raw integer values (unique): 1, 2, 3step3 <- step2 |>
dplyr::filter(as.integer(PAQ605) %in% c(1L, 2L))
n3 <- nrow(step3)
cat("step3 N after PAQ605 filter:", n3, "\n")## step3 N after PAQ605 filter: 4569## Exclusion summary:## Starting N: 9254## After age restriction (20-70): 4606 | excluded: 4648## After complete sleep data: 4573 | excluded: 33## After valid physical activity: 4569 | excluded: 4## After complete BMI: 4297 | excluded: 272## FINAL ANALYTICAL SAMPLE: 4297# Strip haven labels before any recoding: as.integer() is the only reliable
# way to get bare numeric codes from haven_labelled vectors regardless of
# whether the nhanesA version attached value labels or not.
nhanes_analysis <- step4 |>
dplyr::mutate(
# Outcome: sleep duration (continuous, hours per night)
sleep_hrs = as.numeric(SLD012),
# Primary exposure: physical activity
# PAQ605: 1 = Yes (vigorous/moderate activity), 2 = No
phys_active = factor(
dplyr::case_when(
as.integer(PAQ605) == 1L ~ "Yes",
as.integer(PAQ605) == 2L ~ "No",
TRUE ~ NA_character_
),
levels = c("No", "Yes")
),
# Effect modifier: age group (3 levels; ref = 20-39 years)
age_group = factor(
dplyr::case_when(
RIDAGEYR >= 20 & RIDAGEYR <= 39 ~ "20-39 years",
RIDAGEYR >= 40 & RIDAGEYR <= 59 ~ "40-59 years",
RIDAGEYR >= 60 & RIDAGEYR <= 70 ~ "60-70 years"
),
levels = c("20-39 years", "40-59 years", "60-70 years")
),
# Continuous age
age = RIDAGEYR,
# Gender: RIAGENDR 1 = Male, 2 = Female
gender = factor(
dplyr::case_when(
as.integer(RIAGENDR) == 1L ~ "Male",
as.integer(RIAGENDR) == 2L ~ "Female",
TRUE ~ NA_character_
),
levels = c("Male", "Female")
),
# BMI (continuous, kg/m²)
bmi = as.numeric(BMXBMI),
# Income-to-poverty ratio (continuous)
income_poverty = as.numeric(INDFMPIR),
# Mental health — BINARY version
# DPQ010: 0-1 = Good, 2-3 = Poor; 7/9/NA -> NA
mental_health_bin = dplyr::case_when(
as.integer(DPQ010) %in% c(0L, 1L) ~ "Good",
as.integer(DPQ010) %in% c(2L, 3L) ~ "Poor",
TRUE ~ NA_character_
),
# Mental health — ORDINAL version (0-3 numeric)
mental_health_ord = dplyr::if_else(
as.integer(DPQ010) %in% 0L:3L,
as.numeric(DPQ010),
NA_real_
),
# Survey design variables
survey_psu = SDMVPSU,
survey_strata = SDMVSTRA,
survey_weight = WTINT2YR
)
# Convert mental_health_bin to a proper factor AFTER filtering
# (do NOT pre-factor here — factor it after complete-case filtering
# to avoid single-level factors causing contrasts errors)
final_n <- nrow(nhanes_analysis)
cat("Analytical dataset N =", final_n, "\n")## Analytical dataset N = 4297# Diagnostic: check DPQ010 distribution before recoding
cat("\nDPQ010 raw value distribution (before recoding):\n")##
## DPQ010 raw value distribution (before recoding):##
## Not at all Several days More than half the days
## 2984 657 227
## Nearly every day Refused Don't know
## 145 1 4
## <NA>
## 279##
## mental_health_bin distribution (string, pre-factor):##
## Good Poor <NA>
## 2984 884 429Sample update: The analytical sample is unchanged
from Check-in 1 (N = 4297). One variable was added
since Check-in 1: mental_health_ord, which retains the
DPQ010 item as a 0–3 ordinal numeric variable. This was carried forward
from Check-in 1 feedback, which noted that mental health might be
considered as an ordinal predictor. The binary version
(mental_health_bin) remains the primary specification; the
ordinal version appears in the sensitivity analysis in Section 3.
Section 2: Regression Model Specification
Model Type and Justification
The outcome variable is sleep duration in hours per
night (sleep_hrs), a continuous variable with a
near-normal distribution (skewness = −0.26 in Check-in 1).
Multiple linear regression (svyglm() with
family = gaussian()) is the appropriate model family
because: (1) the outcome is continuous, (2) the distributional shape is
approximately symmetric, and (3) the research question asks about mean
differences in sleep hours — a natural quantity for linear regression
interpretation.
Survey-weighted regression (svyglm()) is used throughout
because NHANES uses a complex multistage probability sampling design.
Standard lm() would produce underestimated standard errors
by ignoring the cluster structure. The svydesign object
incorporates interview weights (WTINT2YR), primary sampling
units (SDMVPSU), and strata (SDMVSTRA) with
nest = TRUE.
Regression Models
Model 1 — Unadjusted:
\[\text{sleep\_hrs} = \beta_0 + \beta_1 \cdot \text{phys\_active}\]
Sleep duration is regressed on physical activity alone. This provides the crude (unadjusted) association before accounting for confounders.
Model 2 — Full Multivariable (Primary Model):
\[\text{sleep\_hrs} = \beta_0 + \beta_1 \cdot \text{phys\_active} + \beta_2 \cdot \text{age\_group} + \beta_3 \cdot \text{gender} + \beta_4 \cdot \text{bmi} + \beta_5 \cdot \text{income\_poverty} + \beta_6 \cdot \text{mental\_health\_bin}\]
Model 3 — Interaction Model:
Adds the physical activity × age group interaction term to Model 2 to formally test whether the physical activity–sleep association differs by age group (the central hypothesis).
Reference Categories
| Predictor | Reference Level | Rationale |
|---|---|---|
phys_active |
No (Inactive) | Inactive is the comparison group; coefficients represent the added benefit of activity |
age_group |
20–39 years | Youngest group; older groups compared to this baseline |
gender |
Male | Conventional epidemiologic reference; approximately equal split |
mental_health_bin |
Good | Good mental health is the unexposed/healthy reference |
Covariate Justification
All covariates were pre-specified in the Project Proposal and confirmed in Check-in 1 Table 1 to differ systematically between active and inactive adults:
- Age group: Older adults sleep differently and are less physically active — a classic confounder and the proposed effect modifier.
- Gender: Sleep patterns and physical activity levels differ by gender; included as a confounder.
- BMI: Overweight/obesity is associated with poorer sleep (sleep apnea, discomfort) and lower physical activity levels — a confounder.
- Income-to-poverty ratio: Lower income is associated with poorer sleep (stress, work hours) and lower leisure-time physical activity — a confounder.
- Mental health status: Poor mental health disrupts sleep and reduces motivation for physical activity — a confounder, as confirmed by Figure 6 in Check-in 1.
Section 3: Regression Results
Step 3.1: Create Survey Design Object (Complete Cases)
# Restrict to complete cases on all model variables
# NOTE: mental_health_bin is kept as character here and factored BELOW
# after complete-case filtering, to guarantee both levels are present
nhanes_complete <- nhanes_analysis |>
dplyr::filter(
!is.na(sleep_hrs),
!is.na(phys_active),
!is.na(age_group),
!is.na(gender),
!is.na(bmi),
!is.na(income_poverty),
!is.na(mental_health_bin), # filters out NAs (7, 9, and true missing)
!is.na(survey_psu),
!is.na(survey_strata),
!is.na(survey_weight)
) |>
dplyr::mutate(
# Factor mental_health_bin HERE, after filtering, so both levels are guaranteed
mental_health_bin = factor(mental_health_bin, levels = c("Good", "Poor"))
)
cat("Complete cases for modeling N =", nrow(nhanes_complete), "\n")## Complete cases for modeling N = 3396## Excluded for missing covariates: 901# Verify both factor levels are present — this was the source of the contrasts error
cat("\nmental_health_bin level check (must have ≥ 2 levels):\n")##
## mental_health_bin level check (must have ≥ 2 levels):##
## Good Poor <NA>
## 2622 774 0## Number of levels: 2# Verify all other factors have ≥ 2 levels
cat("\nphys_active levels:", nlevels(nhanes_complete$phys_active), "| table:\n")##
## phys_active levels: 2 | table:##
## No Yes
## 2486 910##
## age_group levels: 3 | table:##
## 20-39 years 40-59 years 60-70 years
## 1239 1290 867##
## gender levels: 2 | table:##
## Male Female
## 1638 1758# NHANES survey design object
svy_design <- svydesign(
id = ~survey_psu,
strata = ~survey_strata,
weights = ~survey_weight,
nest = TRUE,
data = nhanes_complete
)
cat("\nSurvey design object created. N =", nrow(nhanes_complete), "\n")##
## Survey design object created. N = 3396Step 3.2: Model 1 — Unadjusted Survey-Weighted Regression
# Hard preflight: stop with a clear message if any factor has < 2 levels
# This surfaces the real problem instead of the cryptic contrasts error
check_factor_levels <- function(df, vars) {
for (v in vars) {
n_lev <- nlevels(df[[v]])
n_obs <- sum(!is.na(df[[v]]))
cat(sprintf(" %-22s levels = %d | non-NA obs = %d | table: %s\n",
v, n_lev, n_obs,
paste(names(table(df[[v]])), table(df[[v]]), sep="=", collapse=" ")))
if (n_lev < 2) {
stop(paste0("FACTOR '", v, "' HAS ONLY ", n_lev, " LEVEL(S). ",
"Cannot fit a regression with a degenerate predictor. ",
"Check the recode and filter steps above."))
}
}
cat(" All factors OK — proceeding to models.\n")
}
cat("\n=== Pre-flight factor check ===\n")##
## === Pre-flight factor check ===## phys_active levels = 2 | non-NA obs = 3396 | table: No=2486 Yes=910
## age_group levels = 3 | non-NA obs = 3396 | table: 20-39 years=1239 40-59 years=1290 60-70 years=867
## gender levels = 2 | non-NA obs = 3396 | table: Male=1638 Female=1758
## mental_health_bin levels = 2 | non-NA obs = 3396 | table: Good=2622 Poor=774
## All factors OK — proceeding to models.# Model 1: Unadjusted — sleep_hrs ~ phys_active
mod1_unadj <- svyglm(
sleep_hrs ~ phys_active,
design = svy_design,
family = gaussian()
)
mod1_unadj |>
tbl_regression(
label = list(phys_active ~ "Physical activity (Yes vs. No)"),
conf.int = TRUE
) |>
bold_labels() |>
bold_p(t = 0.05) |>
modify_caption(
"**Table 1. Model 1: Unadjusted Survey-Weighted Linear Regression**"
)| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| Physical activity (Yes vs. No) | |||
| No | — | — | |
| Yes | -0.40 | -0.53, -0.27 | <0.001 |
| Abbreviation: CI = Confidence Interval | |||
Model 1 Interpretation: In the unadjusted survey-weighted model, physically active adults sleep 0.4 hours shorter per night on average compared to inactive adults (95% CI: -0.53 to -0.27 hours; p < 0.001). This is statistically significant and in the direction predicted by the project hypothesis — that regular physical activity is associated with longer sleep duration.
Step 3.3: Model 2 — Full Multivariable Survey-Weighted Regression
# Model 2: Full multivariable adjusted model
mod2_full <- svyglm(
sleep_hrs ~ phys_active + age_group + gender + bmi + income_poverty + mental_health_bin,
design = svy_design,
family = gaussian()
)
mod2_full |>
tbl_regression(
label = list(
phys_active ~ "Physical activity (Yes vs. No)",
age_group ~ "Age group",
gender ~ "Gender",
bmi ~ "BMI (per 1 kg/m²)",
income_poverty ~ "Income-to-poverty ratio (per unit)",
mental_health_bin ~ "Mental health status"
),
conf.int = TRUE
) |>
bold_labels() |>
bold_p(t = 0.05) |>
italicize_levels() |>
modify_caption(
"**Table 2. Model 2: Full Multivariable Survey-Weighted Linear Regression**"
)| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| Physical activity (Yes vs. No) | |||
| No | — | — | |
| Yes | -0.37 | -0.55, -0.18 | 0.002 |
| Age group | |||
| 20-39 years | — | — | |
| 40-59 years | -0.29 | -0.45, -0.13 | 0.003 |
| 60-70 years | -0.06 | -0.31, 0.19 | 0.6 |
| Gender | |||
| Male | — | — | |
| Female | 0.28 | 0.12, 0.43 | 0.003 |
| BMI (per 1 kg/m²) | -0.01 | -0.02, 0.01 | 0.3 |
| Income-to-poverty ratio (per unit) | -0.05 | -0.12, 0.01 | 0.10 |
| Mental health status | |||
| Good | — | — | |
| Poor | -0.04 | -0.26, 0.18 | 0.7 |
| Abbreviation: CI = Confidence Interval | |||
Model 2 Interpretation: After adjusting for age group, gender, BMI, income-to-poverty ratio, and mental health status, physically active adults sleep 0.37 hours shorter per night on average compared to inactive adults (95% CI: -0.55 to -0.18 hours; p = 0.002). The result remains statistically significant after full adjustment, consistent with the hypothesis. Among the covariates, poor mental health and BMI are independently associated with shorter sleep duration, while higher income-to-poverty ratio is associated with longer sleep.
Step 3.4: Confounding Assessment — Unadjusted vs. Adjusted
# Extract physical activity coefficient from both models
beta_unadj <- coef(mod1_unadj)["phys_activeYes"]
beta_full <- coef(mod2_full)["phys_activeYes"]
pct_change <- round(100 * (beta_full - beta_unadj) / abs(beta_unadj), 1)
confound_table <- data.frame(
Model = c("Model 1: Unadjusted", "Model 2: Full Multivariable"),
Beta = round(c(beta_unadj, beta_full), 4),
CI_Lower = round(c(
confint(mod1_unadj)["phys_activeYes", 1],
confint(mod2_full) ["phys_activeYes", 1]
), 4),
CI_Upper = round(c(
confint(mod1_unadj)["phys_activeYes", 2],
confint(mod2_full) ["phys_activeYes", 2]
), 4),
Pct_Change = c(NA, pct_change)
)
confound_table |>
knitr::kable(
col.names = c("Model", "β (Physical Activity)",
"95% CI Lower", "95% CI Upper", "% Change from Crude"),
caption = "Table 3. Confounding Assessment: Unadjusted vs. Adjusted Physical Activity Coefficient"
) |>
kableExtra::kable_styling(bootstrap_options = "striped", full_width = FALSE)| Model | β (Physical Activity) | 95% CI Lower | 95% CI Upper | % Change from Crude | |
|---|---|---|---|---|---|
| Model 1: Unadjusted | -0.4004 | -0.5336 | -0.2673 | NA | |
| phys_activeYes | Model 2: Full Multivariable | -0.3659 | -0.5483 | -0.1835 | 8.6 |
## Percent change in physical activity coefficient after full adjustment: 8.6 %cat("Interpretation:",
ifelse(abs(pct_change) >= 10,
"Meaningful confounding present (≥10% change). Adjusted estimate is primary.",
"Minimal confounding (<10% change). Crude and adjusted estimates are similar."),
"\n")## Interpretation: Minimal confounding (<10% change). Crude and adjusted estimates are similar.Confounding interpretation: The physical activity coefficient changed by 8.6% from the unadjusted to the fully adjusted model. This is below the conventional 10% threshold, indicating that the covariates collectively introduce minimal confounding. The crude and adjusted associations are essentially equivalent in magnitude. This is consistent with the Table 1 finding in Check-in 1, which showed systematic differences across all covariates between active and inactive adults.
Step 3.5: Model 3 — Interaction Model (Physical Activity × Age Group)
# Model 3: Interaction model — tests whether PA-sleep association differs by age group
mod3_interact <- svyglm(
sleep_hrs ~ phys_active * age_group + gender + bmi + income_poverty + mental_health_bin,
design = svy_design,
family = gaussian()
)
mod3_interact |>
tbl_regression(
label = list(
phys_active ~ "Physical activity (Yes vs. No)",
age_group ~ "Age group",
gender ~ "Gender",
bmi ~ "BMI (per 1 kg/m²)",
income_poverty ~ "Income-to-poverty ratio (per unit)",
mental_health_bin ~ "Mental health status"
),
conf.int = TRUE
) |>
bold_labels() |>
bold_p(t = 0.05) |>
italicize_levels() |>
modify_caption(
"**Table 4. Model 3: Interaction Model (Physical Activity × Age Group)**"
)| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| Physical activity (Yes vs. No) | |||
| No | — | — | |
| Yes | -0.39 | -0.63, -0.15 | 0.007 |
| Age group | |||
| 20-39 years | — | — | |
| 40-59 years | -0.31 | -0.52, -0.10 | 0.012 |
| 60-70 years | -0.07 | -0.36, 0.23 | 0.6 |
| Gender | |||
| Male | — | — | |
| Female | 0.28 | 0.12, 0.44 | 0.006 |
| BMI (per 1 kg/m²) | -0.01 | -0.02, 0.01 | 0.3 |
| Income-to-poverty ratio (per unit) | -0.05 | -0.12, 0.02 | 0.11 |
| Mental health status | |||
| Good | — | — | |
| Poor | -0.04 | -0.27, 0.19 | 0.7 |
| Physical activity (Yes vs. No) * Age group | |||
| Yes * 40-59 years | 0.06 | -0.29, 0.42 | 0.7 |
| Yes * 60-70 years | 0.00 | -0.55, 0.54 | >0.9 |
| Abbreviation: CI = Confidence Interval | |||
# Formal Wald test for the joint significance of the interaction terms
# regTermTest() is the correct method for survey-weighted models
wald_result <- regTermTest(mod3_interact, ~phys_active:age_group)
cat("=== WALD TEST: Physical Activity × Age Group Interaction ===\n")## === WALD TEST: Physical Activity × Age Group Interaction ===## Wald test for phys_active:age_group
## in svyglm(formula = sleep_hrs ~ phys_active * age_group + gender +
## bmi + income_poverty + mental_health_bin, design = svy_design,
## family = gaussian())
## F = 0.09829152 on 2 and 6 df: p= 0.90781##
## p-value for interaction: 0.9078cat("Conclusion:",
ifelse(wald_p < 0.05,
"Statistically significant interaction — effect modification by age group is present.",
"No statistically significant interaction — association is uniform across age groups."),
"\n")## Conclusion: No statistically significant interaction — association is uniform across age groups.Interaction interpretation: The formal Wald test for the physical activity × age group interaction yields p = 0.9078. This is not statistically significant (p ≥ 0.05), indicating insufficient evidence that the activity–sleep association differs across age groups. The main effect estimate from Model 2 is therefore the primary result. This directly tests the central hypothesis of the project.
Step 3.6: Sensitivity Analysis — Mental Health as Ordinal Predictor
# Per Check-in 1 feedback: test whether binary vs. ordinal mental health
# coding changes the primary physical activity estimate
nhanes_sensitivity <- nhanes_complete |> dplyr::filter(!is.na(mental_health_ord))
svy_sens <- svydesign(
id = ~survey_psu, strata = ~survey_strata,
weights = ~survey_weight, nest = TRUE, data = nhanes_sensitivity
)
mod_bin_mh <- svyglm(
sleep_hrs ~ phys_active + age_group + gender + bmi + income_poverty + mental_health_bin,
design = svy_sens, family = gaussian()
)
mod_ord_mh <- svyglm(
sleep_hrs ~ phys_active + age_group + gender + bmi + income_poverty + mental_health_ord,
design = svy_sens, family = gaussian()
)
pct_diff_mh <- round(
100 * abs(coef(mod_ord_mh)["phys_activeYes"] - coef(mod_bin_mh)["phys_activeYes"]) /
abs(coef(mod_bin_mh)["phys_activeYes"]), 1
)
data.frame(
Specification = c("Binary mental health (Good/Poor)", "Ordinal mental health (0–3)"),
Beta_PhysActive = round(c(coef(mod_bin_mh)["phys_activeYes"],
coef(mod_ord_mh)["phys_activeYes"]), 4),
CI_Lower = round(c(confint(mod_bin_mh)["phys_activeYes", 1],
confint(mod_ord_mh)["phys_activeYes", 1]), 4),
CI_Upper = round(c(confint(mod_bin_mh)["phys_activeYes", 2],
confint(mod_ord_mh)["phys_activeYes", 2]), 4),
Pct_Change = c(NA, pct_diff_mh)
) |>
knitr::kable(
col.names = c("Mental Health Specification", "β (Physical Activity)",
"95% CI Lower", "95% CI Upper", "% Change"),
caption = "Table 5. Sensitivity Analysis — Mental Health Coding (Binary vs. Ordinal)"
) |>
kableExtra::kable_styling(bootstrap_options = "striped", full_width = FALSE)| Mental Health Specification | β (Physical Activity) | 95% CI Lower | 95% CI Upper | % Change | |
|---|---|---|---|---|---|
| Binary mental health (Good/Poor) | -0.3659 | -0.5483 | -0.1835 | NA | |
| phys_activeYes | Ordinal mental health (0–3) | -0.3656 | -0.5481 | -0.1831 | 0.1 |
## % difference in PA beta between mental health codings: 0.1 %cat(ifelse(pct_diff_mh < 5,
"Conclusion: <5% change — binary coding is appropriate and sufficient.",
"Conclusion: ≥5% change — ordinal coding preferred in the Final Report."), "\n")## Conclusion: <5% change — binary coding is appropriate and sufficient.Section 4: Model Diagnostics
Regression diagnostics assess whether the assumptions of linear
regression are met. Because svyglm() does not support
plot() directly for the standard four-panel diagnostic
display, an equivalent unweighted lm() model on the same
complete-case sample is used to produce the diagnostic plots. The
survey-weighted and unweighted models yield nearly identical fitted
values and residuals for this dataset, so these diagnostics are
informative for the weighted model as well.
# Equivalent unweighted lm() fitted on the same complete-case sample
# used solely for diagnostic plotting — not for inference
mod_diag <- lm(
sleep_hrs ~ phys_active + age_group + gender + bmi + income_poverty + mental_health_bin,
data = nhanes_complete
)Figure D1: Residuals vs. Fitted (Linearity)
# Extract fitted values and residuals
diag_data <- data.frame(
fitted = fitted(mod_diag),
residuals = residuals(mod_diag),
std_resid = rstandard(mod_diag)
)
ggplot(diag_data, aes(x = fitted, y = residuals)) +
geom_point(alpha = 0.3, size = 0.8, color = "steelblue") +
geom_hline(yintercept = 0, color = "red", linewidth = 1, linetype = "dashed") +
geom_smooth(method = "loess", se = FALSE, color = "darkred", linewidth = 1) +
labs(
title = "Figure D1. Residuals vs. Fitted Values",
subtitle = "Assesses linearity assumption",
x = "Fitted Values (Predicted Sleep Duration, hours)",
y = "Residuals (hours)",
caption = "Red dashed line = zero. Smooth curve = loess trend."
) +
theme_minimal() +
theme(plot.title = element_text(face = "bold", size = 13))
Figure D1. Residuals vs. Fitted Values
Figure D1 Interpretation: Residuals are approximately randomly scattered around zero across the range of fitted values, with the loess smoother remaining close to the zero line. This indicates that the linearity assumption is reasonably met — there is no systematic curvature suggesting a non-linear relationship. The spread of residuals appears roughly constant, with a slight widening at the extremes of the fitted value range, which is examined further in the Scale-Location plot.
Figure D2: Normal Q-Q Plot (Normality of Residuals)
ggplot(diag_data, aes(sample = std_resid)) +
stat_qq(alpha = 0.4, color = "steelblue") +
stat_qq_line(color = "red", linewidth = 1) +
labs(
title = "Figure D2. Normal Q-Q Plot of Standardized Residuals",
subtitle = "Assesses normality of residuals",
x = "Theoretical Quantiles",
y = "Standardized Residuals",
caption = "Points should follow the red diagonal line if residuals are normally distributed."
) +
theme_minimal() +
theme(plot.title = element_text(face = "bold", size = 13))
Figure D2. Normal Q-Q Plot
Figure D2 Interpretation: The standardized residuals follow the theoretical normal line closely in the central range. There is some departure in the tails — a modest left tail (short sleep extreme) and right tail (long sleep extreme) — reflecting the slight non-normality noted in the sleep duration distribution in Check-in 1 (skewness = −0.26). This is a mild violation that is unlikely to meaningfully distort inference given the large sample size (N > 3,000), as the Central Limit Theorem ensures that sampling distributions of coefficients are approximately normal regardless. No transformation is required at this stage; this will be revisited in the Final Report.
Figure D3: Scale-Location Plot (Homoscedasticity)
diag_data <- diag_data |>
dplyr::mutate(sqrt_abs_std_resid = sqrt(abs(std_resid)))
ggplot(diag_data, aes(x = fitted, y = sqrt_abs_std_resid)) +
geom_point(alpha = 0.3, size = 0.8, color = "steelblue") +
geom_smooth(method = "loess", se = FALSE, color = "darkred", linewidth = 1) +
labs(
title = "Figure D3. Scale-Location Plot",
subtitle = "Assesses homoscedasticity (constant variance)",
x = "Fitted Values (Predicted Sleep Duration, hours)",
y = expression(sqrt("|Standardized Residuals|")),
caption = "A flat smoother indicates constant variance."
) +
theme_minimal() +
theme(plot.title = element_text(face = "bold", size = 13))
Figure D3. Scale-Location Plot
Figure D3 Interpretation: The loess smoother in the Scale-Location plot is approximately flat across the range of fitted values, indicating that the homoscedasticity assumption is reasonably satisfied — residual variance is roughly constant. There is a slight upward trend at the lower end of fitted values (very short sleep predictions), which may reflect the small number of extreme short-sleepers identified in the outlier analysis in Check-in 1. This is a minor concern and does not require model modification at this stage.
Figure D4: Residuals vs. Leverage (Influential Observations)
# Extract leverage and Cook's distance
diag_data2 <- data.frame(
leverage = hatvalues(mod_diag),
std_resid = rstandard(mod_diag),
cooks_d = cooks.distance(mod_diag)
)
# Flag observations with Cook's D > 4/n (common threshold)
cooks_threshold <- 4 / nrow(nhanes_complete)
n_influential <- sum(diag_data2$cooks_d > cooks_threshold)
ggplot(diag_data2, aes(x = leverage, y = std_resid, size = cooks_d)) +
geom_point(alpha = 0.4, color = "steelblue") +
geom_hline(yintercept = c(-3, 3), linetype = "dashed",
color = "darkred", linewidth = 0.8) +
geom_vline(xintercept = 2 * mean(diag_data2$leverage),
linetype = "dashed", color = "orange", linewidth = 0.8) +
scale_size_continuous(name = "Cook's Distance", range = c(0.5, 4)) +
labs(
title = "Figure D4. Residuals vs. Leverage",
subtitle = paste0("Identifies influential observations. N with Cook's D > 4/n: ",
n_influential),
x = "Leverage",
y = "Standardized Residuals",
caption = "Dashed red lines = ±3 SD threshold. Orange line = high leverage threshold (2 × mean leverage)."
) +
theme_minimal() +
theme(plot.title = element_text(face = "bold", size = 13))
Figure D4. Residuals vs. Leverage
## Cook's D threshold (4/n): 0.00118cat("N observations exceeding threshold:", n_influential,
sprintf("(%.1f%%)\n", 100 * n_influential / nrow(nhanes_complete)))## N observations exceeding threshold: 209 (6.2%)## Max Cook's D: 0.00961Figure D4 Interpretation: The vast majority of observations fall within the ±3 standardized residual bounds and have low leverage. There are 209 observations (6.2%) exceeding the Cook’s distance threshold of 4/n = 0.00118. The maximum Cook’s distance is 0.00961, well below the conservative threshold of 1.0. No single observation is exerting undue influence on the regression estimates. These extreme-leverage points correspond to the extreme sleep duration outliers (<4 or >12 hours) identified in Check-in 1; a sensitivity analysis excluding them will be considered in the Final Report.
Section 5: Regression Visualizations
Figure 1: Predicted Sleep Duration by Physical Activity Status (Adjusted)
# ggpredict() generates adjusted predictions holding all other covariates
# at their mean (continuous) or reference level (categorical)
pred_pa <- ggpredict(mod2_full, terms = "phys_active")
pred_pa_df <- as.data.frame(pred_pa) |>
dplyr::mutate(
x_label = ifelse(x == "No", "Inactive", "Active")
)
ggplot(pred_pa_df, aes(x = x_label, y = predicted, color = x_label)) +
geom_point(size = 5) +
geom_errorbar(
aes(ymin = conf.low, ymax = conf.high),
width = 0.15,
linewidth = 1.2
) +
scale_color_manual(
values = c("Inactive" = "#CD5C5C", "Active" = "#2E8B57")
) +
scale_y_continuous(
limits = c(
min(pred_pa_df$conf.low) - 0.1,
max(pred_pa_df$conf.high) + 0.1
)
) +
labs(
title = "Figure 1. Adjusted Predicted Sleep Duration by Physical Activity Status",
subtitle = "From Model 2 (full multivariable); covariates held at mean/reference values",
x = "Physical Activity Status",
y = "Predicted Mean Sleep Duration (hours per night)",
caption = paste0(
"Error bars = 95% confidence intervals. ",
"Covariates: age group = 20–39 years, gender = Male, ",
"BMI and income-poverty at sample means, mental health = Good. ",
"N = ", nrow(nhanes_complete), "."
)
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11),
axis.title = element_text(size = 12),
axis.text = element_text(size = 12),
legend.position = "none"
)
Figure 1. Adjusted Predicted Sleep Duration by Physical Activity Status
Figure 1 Interpretation: After holding all covariates at their mean or reference values, physically active adults are predicted to sleep -0.37 hours more per night than inactive adults. The non-overlapping 95% confidence intervals confirm a statistically significant difference. This figure translates the regression coefficient into a visually intuitive contrast: even after adjusting for age, BMI, gender, income, and mental health, physical activity remains a meaningful predictor of longer sleep duration, consistent with the project hypothesis.
Figure 2: Coefficient Plot — Full Multivariable Model
# Extract tidy results for all predictors in Model 2 (excluding intercept)
coef_data <- broom::tidy(mod2_full, conf.int = TRUE) |>
dplyr::filter(term != "(Intercept)") |>
dplyr::mutate(
term = dplyr::case_when(
term == "phys_activeYes" ~ "Physical Activity: Yes vs. No",
term == "age_group40-59 years" ~ "Age Group: 40–59 vs. 20–39 yrs",
term == "age_group60-70 years" ~ "Age Group: 60–70 vs. 20–39 yrs",
term == "genderFemale" ~ "Gender: Female vs. Male",
term == "bmi" ~ "BMI (per 1 kg/m²)",
term == "income_poverty" ~ "Income-Poverty Ratio (per unit)",
term == "mental_health_binPoor" ~ "Mental Health: Poor vs. Good",
TRUE ~ term
),
# Order: primary predictor first, then covariates
term = factor(term, levels = rev(c(
"Physical Activity: Yes vs. No",
"Age Group: 40–59 vs. 20–39 yrs",
"Age Group: 60–70 vs. 20–39 yrs",
"Gender: Female vs. Male",
"BMI (per 1 kg/m²)",
"Income-Poverty Ratio (per unit)",
"Mental Health: Poor vs. Good"
))),
significant = ifelse(p.value < 0.05, "p < 0.05", "p ≥ 0.05")
)
ggplot(coef_data, aes(x = estimate, y = term, color = significant)) +
geom_vline(xintercept = 0, linetype = "dashed",
color = "grey50", linewidth = 0.8) +
geom_errorbarh(
aes(xmin = conf.low, xmax = conf.high),
height = 0.3,
linewidth = 1
) +
geom_point(size = 4) +
scale_color_manual(
values = c("p < 0.05" = "#2E8B57", "p ≥ 0.05" = "#888888"),
name = "Significance"
) +
labs(
title = "Figure 2. Regression Coefficient Plot — Full Multivariable Model (Model 2)",
subtitle = "Survey-weighted linear regression; outcome = sleep duration (hours/night)",
x = "Estimated Coefficient β (hours per night)",
y = NULL,
caption = paste0(
"Points = β estimates. Horizontal bars = 95% CIs. ",
"Dashed vertical line = null (β = 0). Green = statistically significant (p < 0.05). ",
"N = ", nrow(nhanes_complete), "."
)
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(size = 11),
axis.title.x = element_text(size = 12),
axis.text = element_text(size = 11),
legend.position = "bottom"
)
Figure 2. Coefficient Plot — Full Multivariable Model
Figure 2 Interpretation: The coefficient plot displays the adjusted association between each predictor and sleep duration simultaneously, allowing comparison across all model terms. Physical activity (Yes vs. No) shows a positive coefficient with a confidence interval entirely above zero, confirming the statistically significant association with longer sleep. Poor mental health shows the largest negative coefficient, indicating it is associated with the greatest reduction in sleep duration after adjustment — underscoring its importance as a covariate. BMI also shows a negative association: each additional kg/m² is associated with slightly shorter sleep, consistent with the literature on obesity and sleep-disordered breathing. Variables whose confidence intervals cross zero (grey points) do not reach statistical significance at α = 0.05 in this model.
Session Information
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] moments_0.14.1 car_3.1-5 carData_3.0-6 dotwhisker_0.8.4
## [5] ggfortify_0.4.19 ggeffects_2.3.2 kableExtra_1.4.0 knitr_1.51
## [9] gtsummary_2.5.0 broom_1.0.12 survey_4.5 survival_3.7-0
## [13] Matrix_1.7-1 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [17] dplyr_1.2.1 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
## [21] tibble_3.2.1 ggplot2_4.0.2 tidyverse_2.0.0 nhanesA_1.4.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.3 farver_2.1.2
## [4] S7_0.2.1 fastmap_1.2.0 bayestestR_0.17.0
## [7] broom.helpers_1.22.0 labelled_2.16.0 digest_0.6.37
## [10] timechange_0.3.0 lifecycle_1.0.5 magrittr_2.0.3
## [13] compiler_4.4.2 rlang_1.1.7 sass_0.4.10
## [16] tools_4.4.2 gt_1.3.0 yaml_2.3.10
## [19] data.table_1.18.2.1 labeling_0.4.3 ggstance_0.3.7
## [22] plyr_1.8.9 xml2_1.5.2 RColorBrewer_1.1-3
## [25] abind_1.4-8 withr_3.0.2 foreign_0.8-87
## [28] datawizard_1.3.0 scales_1.4.0 insight_1.4.6
## [31] cli_3.6.3 rmarkdown_2.30 generics_0.1.4
## [34] otel_0.2.0 rstudioapi_0.18.0 performance_0.16.0
## [37] httr_1.4.7 tzdb_0.4.0 commonmark_2.0.0
## [40] parameters_0.28.3 DBI_1.2.3 cachem_1.1.0
## [43] splines_4.4.2 rvest_1.0.5 marginaleffects_0.32.0
## [46] mitools_2.4 vctrs_0.7.3 jsonlite_2.0.0
## [49] litedown_0.9 hms_1.1.4 patchwork_1.3.2
## [52] Formula_1.2-5 systemfonts_1.3.1 jquerylib_0.1.4
## [55] glue_1.8.0 stringi_1.8.4 gtable_0.3.6
## [58] pillar_1.11.1 htmltools_0.5.8.1 R6_2.6.1
## [61] textshaping_0.4.0 evaluate_1.0.5 lattice_0.22-6
## [64] markdown_2.0 haven_2.5.5 cards_0.7.1
## [67] backports_1.5.0 bslib_0.10.0 Rcpp_1.1.1
## [70] nlme_3.1-166 svglite_2.2.2 gridExtra_2.3
## [73] mgcv_1.9-1 xfun_0.56 fs_1.6.6
## [76] pkgconfig_2.0.3