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

Setup and Data

In-Class Lab Activity

Data for the Lab

Variable Description Type
fmd Frequent mental distress (No/Yes) Factor (outcome)
menthlth_days Mentally unhealthy days (0–30) Numeric
physhlth_days Physically unhealthy days (0–30) Numeric
sleep_hrs Sleep hours per night (1–14) Numeric
age Age in years (capped at 80) Numeric
sex Sex (Male/Female) Factor
bmi Body mass index Numeric
exercise Exercised in past 30 days (No/Yes) Factor
income_cat Household income category (1–8) Numeric
smoker Former/Never vs. Current Factor 
library(tidyverse)
library(knitr)
library(kableExtra)
library(broom)
library(gtsummary)
library(car)
library(ggeffects)
library(ResourceSelection)  # for Hosmer-Lemeshow
library(pROC)               # for ROC/AUC
library(performance)        # for model performance metrics
library(sjPlot)
library(modelsummary)

brfss_logistic <- readRDS(
  "C:/Users/graci/OneDrive/Documents/UA GRAD SCHOOL/2nd Semester/EPI553/BRFSS_2020.rds"
)

Descriptive Overview

brfss_logistic |>
  tbl_summary(
    by = fmd,
    include = c(physhlth_days, sleep_hrs, age, sex, bmi, exercise,
                income_cat, smoker),
    type = list(
      c(physhlth_days, sleep_hrs, age, bmi, income_cat) ~ "continuous"
    ),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})"
    ),
    label = list(
      physhlth_days ~ "Physical unhealthy days",
      sleep_hrs     ~ "Sleep hours",
      age           ~ "Age (years)",
      sex           ~ "Sex",
      bmi           ~ "BMI",
      exercise      ~ "Exercise in past 30 days",
      income_cat    ~ "Income category (1-8)",
      smoker        ~ "Smoking status"
    )
  ) |>
  add_overall() |>
  add_p() |>
  bold_labels()
Characteristic Overall
N = 5,0001		 No
N = 4,2431		 Yes
N = 7571		 p-value
Physical unhealthy days 4 (9) 3 (8) 10 (13)
Sleep hours 7.00 (1.48) 7.09 (1.40) 6.51 (1.83)
Age (years) 56 (16) 57 (16) 50 (16)
Sex



    Male 2,701 (54%) 2,378 (56%) 323 (43%)
    Female 2,299 (46%) 1,865 (44%) 434 (57%)
BMI 28.5 (6.3) 28.4 (6.2) 29.3 (7.0)
Exercise in past 30 days 3,673 (73%) 3,192 (75%) 481 (64%)
Income category (1-8) 5.85 (2.11) 6.00 (2.04) 5.05 (2.29)
Smoking status



    Former/Never 3,280 (66%) 2,886 (68%) 394 (52%)
    Current 1,720 (34%) 1,357 (32%) 363 (48%)
1 Mean (SD); n (%)

Task 1: Build a Multiple Logistic Regression Model (15 points)

1a. (5 pts) Fit a multiple logistic regression model predicting fmd from at least 5 predictors of your choice.

mod_fmd <- glm(fmd ~ exercise + smoker + sex + sleep_hrs + income_cat,data = brfss_logistic,family = binomial)

1b. (5 pts) Report the adjusted ORs with 95% CIs in a publication-quality table using tbl_regression().

mod_fmd <- glm(fmd ~ exercise + smoker + sex + sleep_hrs + income_cat,data = brfss_logistic,family = binomial) 

mod_fmd |>
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      exercise   ~ "Exercise",
      smoker     ~ "Smoking Status",
      sex        ~ "Sex",
      sleep_hrs  ~ "Sleep Hours",
      income_cat ~ "Income category (1-8)"
      
    )
  ) |>
  bold_labels() |>
  bold_p()
Characteristic OR 95% CI p-value
Exercise


    No — —
    Yes 0.69 0.58, 0.82 <0.001
Smoking Status


    Former/Never — —
    Current 1.55 1.31, 1.82 <0.001
Sex


    Male — —
    Female 1.61 1.37, 1.90 <0.001
Sleep Hours 0.79 0.75, 0.83 <0.001
Income category (1-8) 0.86 0.83, 0.90 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

1c. (5 pts) Interpret the adjusted OR for two predictors of your choice in 1-2 sentences each. Make sure to mention what the OR represents (per unit change for continuous; reference category for categorical).

The adjusted odds ratio for sleep hours is 0.79 (95% CI: 0.75,0.83). This means that for each additional hour of sleep, the odds of frequent mental distress decreases by 21%, holding all other variables constant.

The adjusted odds ratio for smoking status is 1.55(95% CI:1.31,1.82), comparing current smokers to former/never smokers(reference group). This shows that current smokers have 55% higher odds of frequent mental distress compared to former/never smokers.

Task 2: Wald and Likelihood Ratio Tests (15 points)

2a. (5 pts) Identify the Wald p-value for each predictor in your model from the tidy() or summary() output.

mod_fmd <- glm(fmd ~ exercise + smoker + sex + sleep_hrs + income_cat,data = brfss_logistic,family = binomial) 


tidy(mod_fmd, conf.int = TRUE, exponentiate = TRUE) |>
  mutate(across(c(estimate, std.error, statistic, conf.low, conf.high),
                \(x) round(x, 3)),
         p.value = format.pval(p.value, digits = 3)) |>
  kable(caption = "Wald Tests for Each Coefficient (Exponentiated)") |>
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Wald Tests for Each Coefficient (Exponentiated)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 1.705 0.224 2.379 0.0173 1.098 2.645
exerciseYes 0.694 0.088 -4.174 2.99e-05 0.585 0.824
smokerCurrent 1.546 0.084 5.185 2.16e-07 1.311 1.823
sexFemale 1.613 0.082 5.810 6.25e-09 1.373 1.896
sleep_hrs 0.790 0.027 -8.805 < 2e-16 0.749 0.832
income_cat 0.864 0.019 -7.693 1.44e-14 0.833 0.897

2b. (5 pts) Fit a reduced model that drops one predictor of your choice. Perform a likelihood ratio test comparing the full and reduced models using anova(reduced, full, test = "LRT").

mod_reduced <- glm(fmd ~ exercise + sex + sleep_hrs + income_cat,data = brfss_logistic,family = binomial)

anova(mod_reduced, mod_fmd, test = "LRT") |>
  kable(digits = 3,
        caption = "LR Test: Does removing smoker improve the model?") |>
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
LR Test: Does removing smoker improve the model?
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
4995 3989.548 NA NA NA
4994 3962.991 1 26.557 0

2c. (5 pts) Compare the conclusions from the Wald test (for the dropped predictor) and the LR test. Do they agree? In 2-3 sentences, explain when the two tests might disagree.

The Wald test for smoker shows a highly statistical significant association with frequent mental distress(p=2.16e^07), showing that smoking status is an important predictor in this model. The likelihood ratio test compares the reduced and full model which yields a highly significant result, showing that removing smoker variable significant worsens model fit. The two tests agree in that smoking status should be retained in the model. The tests might disagree when sample sizes are small.

Task 3: Test an Interaction (20 points)

3a. (5 pts) Fit a model that includes an interaction between two predictors of your choice (e.g., exercise * sex or smoker * age).

mod_int <- glm(fmd ~ exercise * sex,data = brfss_logistic,family = binomial)

3b. (5 pts) Perform a likelihood ratio test comparing the model with the interaction to the model without it.

mod_fmd <- glm(fmd ~ exercise + smoker + sex + sleep_hrs + income_cat,data = brfss_logistic,family = binomial) 

mod_int <- glm(fmd ~ exercise* sex,data = brfss_logistic,family = binomial) 

anova(mod_int, mod_fmd, test = "LRT") |>
  kable(digits = 3,
        caption = "LR Test: Does exercise*sex improve the model") |>
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
LR Test: Does exercise*sex improve the model
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
4996 4165.132 NA NA NA
4994 3962.991 2 202.141 0

3c. (5 pts) Create a visualization of the interaction using ggpredict() and plot().

mod_int <- glm(fmd ~ exercise* sex,data = brfss_logistic,family = binomial) 

ggpredict(mod_int, terms = c("exercise", "sex")) |>
  plot() +
  labs(title = "Predicted Probability of FMD by Exercise and Sex",
       x = "Exercise", y = "Predicted Probability of FMD",
       color = "Sex") +
  theme_minimal()

3d. (5 pts) In 3-4 sentences, interpret the interaction. Does the effect of one predictor differ across levels of the other? If statistically significant, report the stratum-specific odds ratios.

The interaction plot shows that the effect of exercise on frequent mental distress differs somewhat by sex. For both male and females who exercised has a lower predicted probability of frequent mental distress compared to those who don’t exercise. However, females consistently have higher predicted probabilities of FMD than males at both levels of exercise. The likelihood ratio test comparing the model with and without the interaction term is statistically significant(x^2=202.14, p<0.001), showing us that including the interaction improves model fit.

Task 4: Goodness-of-Fit and Discrimination (25 points)

4a. (5 pts) Compute McFadden’s pseudo-R² for your full model using performance::r2_mcfadden().

performance::r2_mcfadden(mod_fmd)
## # R2 for Generalized Linear Regression
##        R2: 0.068
##   adj. R2: 0.067

4b. (5 pts) Perform the Hosmer-Lemeshow goodness-of-fit test using ResourceSelection::hoslem.test(). Report the test statistic and p-value. Comment on the interpretation given your sample size.

The hosmer-lemeshow goodness of fit test showed a statistic of X^2=14.80 with 8 degrees of freedom and a p-value of 0.063. The p-value is greater than 0.05, showing that it is not statistically significant which indicates adequate fit. Given the large sample size, the hosmer-lemeshow test can be over powered and detect trivial misfit.

hl_test <- hoslem.test(
  x = as.numeric(brfss_logistic$fmd) - 1,
  y = fitted(mod_fmd),
  g = 10
)

hl_test
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  as.numeric(brfss_logistic$fmd) - 1, fitted(mod_fmd)
## X-squared = 14.801, df = 8, p-value = 0.06314

4c. (5 pts) Create a calibration plot showing observed vs. predicted probabilities by decile. Comment on whether the model appears well calibrated.

From the calibration plot it shows that observed proportions of frequent mental distress generally tracks closely to the predicted probabilities. Most of the points lie near the dashed line indicating that their is a good agreement between observed and predicted values. There are some minor deviations in the lower probability ranges, where the model slightly over or under predicts the outcome. These differences are relatively small and not systematic.

brfss_logistic |>
  mutate(pred_prob = fitted(mod_fmd),
         obs_outcome = as.numeric(fmd) - 1,
         decile = ntile(pred_prob, 10)) |>
  group_by(decile) |>
  summarise(
    mean_pred = mean(pred_prob),
    mean_obs  = mean(obs_outcome),
    .groups = "drop"
  ) |>
  ggplot(aes(x = mean_pred, y = mean_obs)) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  geom_point(size = 3, color = "steelblue") +
  geom_line(color = "steelblue") +
  labs(title = "Calibration Plot: Observed vs. Predicted Probability of FMD",
       subtitle = "Points should fall on the dashed line for perfect calibration",
       x = "Mean Predicted Probability (within decile)",
       y = "Observed Proportion (within decile)") +
  theme_minimal()

4d. (10 pts) Compute and plot the ROC curve using pROC::roc(). Report the AUC. Based on the AUC value, how would you describe the model’s discrimination ability (poor, acceptable, excellent, outstanding)?

The AUC for the model is 0.689. This shows that the model has poor acceptable discrimination, meaning that it is limited in its ability to distinguish between individuals with and without FMD.

roc_obj <- roc(
  response = brfss_logistic$fmd,
  predictor = fitted(mod_fmd),
  levels = c("No", "Yes"),
  direction = "<"
)

auc_value <- auc(roc_obj)

ggroc(roc_obj, color = "steelblue", linewidth = 1.2) +
  geom_abline(slope = 1, intercept = 1, linetype = "dashed", color = "red") +
  labs(title = "ROC Curve for Frequent Mental Distress Model",
       subtitle = paste0("AUC = ", round(auc_value, 3)),
       x = "Specificity", y = "Sensitivity") +
  theme_minimal()

End of Lab Activity