Week 3: Model
Diabetes Complication Classifier
# Rebuild the patient table from the raw data so this notebook runs on its own
pt <- load_t2dm(data_path) %>% build_patient_table() %>% analytic_cohort()
# Features: required set. To include your Week 2 stretch features, use
# features <- c(model_features, stretch_features)
features <- model_features
features## [1] "age_at_index" "Patient_sex" "Patient_Race" "insurance"
## [5] "main_specialty" "n_encounters" "n_other_dx" "hypertension"
## [9] "obesity" "hyperlipidemia" "ckd" "heart_failure"1. Stratified 75/25 split
sp <- stratified_split(pt, prop = 0.75, seed = 719)
train <- sp$train
test <- sp$test
bind_rows(train = train %>% summarise(n = n(), pct_outcome = 100 * mean(any_complication)),
test = test %>% summarise(n = n(), pct_outcome = 100 * mean(any_complication)),
.id = "set")# Classification threshold = training prevalence (see helpers for rationale)
threshold <- mean(train$any_complication)
threshold## [1] 0.11557022. Logistic regression
logit_fit <- fit_logit(train, features)
summary(logit_fit)##
## Call:
## glm(formula = make_formula(features), family = binomial(), data = train)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -2.247e+00 9.224e-01 -2.436
## age_at_index -9.089e-04 2.180e-03 -0.417
## Patient_sexMale 2.642e-02 5.093e-02 0.519
## Patient_RaceAsian -5.180e-02 1.969e-01 -0.263
## Patient_RaceBlack or African American -4.584e-02 5.930e-02 -0.773
## Patient_RaceMore than one race 1.967e-01 1.305e-01 1.508
## Patient_RaceNative Hawaiian / Pacific Islander 8.192e-02 2.173e-01 0.377
## Patient_RaceUnknown 1.194e-01 1.735e-01 0.688
## Patient_RaceUnreported / Declined -2.536e-02 8.770e-02 -0.289
## Patient_RaceOther race 1.507e-01 4.362e-01 0.346
## insuranceCommercial Plan A 7.808e-02 8.344e-02 0.936
## insuranceCommercial Plan B -4.984e-02 8.612e-02 -0.579
## insuranceMedicaid 1.113e-01 8.062e-02 1.380
## insuranceMedicare Advantage -1.492e-01 1.434e-01 -1.040
## insuranceSelf-pay / Uninsured 1.147e-01 9.876e-02 1.162
## insuranceSliding Fee Scale 8.507e-02 1.572e-01 0.541
## insuranceOther payer 3.407e-02 3.685e-01 0.092
## main_specialtyBehavioral Health -8.799e-01 5.417e-01 -1.624
## main_specialtyBurn Surgery 7.125e-01 4.290e-01 1.661
## main_specialtyCardiology -3.118e-01 4.041e-01 -0.771
## main_specialtyDermatology 6.880e-01 1.311e+00 0.525
## main_specialtyGeriatric Primary Care -1.465e-01 1.994e-01 -0.735
## main_specialtyObstetrics and Gynecology 2.931e-01 2.552e-01 1.148
## main_specialtyPodiatry 7.959e-02 1.990e-01 0.400
## main_specialtyOther specialty 2.111e-01 1.083e-01 1.949
## n_encounters 1.362e+00 9.109e-01 1.495
## n_other_dx -7.697e+00 7.112e-01 -10.822
## hypertension -9.961e+00 2.282e+02 -0.044
## obesity -9.985e+00 3.051e+02 -0.033
## hyperlipidemia -9.973e+00 3.923e+02 -0.025
## ckd -1.018e+01 1.071e+03 -0.009
## heart_failure -9.946e+00 9.111e+02 -0.011
## Pr(>|z|)
## (Intercept) 0.0148 *
## age_at_index 0.6767
## Patient_sexMale 0.6039
## Patient_RaceAsian 0.7924
## Patient_RaceBlack or African American 0.4395
## Patient_RaceMore than one race 0.1317
## Patient_RaceNative Hawaiian / Pacific Islander 0.7061
## Patient_RaceUnknown 0.4914
## Patient_RaceUnreported / Declined 0.7725
## Patient_RaceOther race 0.7297
## insuranceCommercial Plan A 0.3494
## insuranceCommercial Plan B 0.5628
## insuranceMedicaid 0.1675
## insuranceMedicare Advantage 0.2982
## insuranceSelf-pay / Uninsured 0.2454
## insuranceSliding Fee Scale 0.5885
## insuranceOther payer 0.9263
## main_specialtyBehavioral Health 0.1043
## main_specialtyBurn Surgery 0.0968 .
## main_specialtyCardiology 0.4405
## main_specialtyDermatology 0.5997
## main_specialtyGeriatric Primary Care 0.4625
## main_specialtyObstetrics and Gynecology 0.2509
## main_specialtyPodiatry 0.6892
## main_specialtyOther specialty 0.0513 .
## n_encounters 0.1350
## n_other_dx <2e-16 ***
## hypertension 0.9652
## obesity 0.9739
## hyperlipidemia 0.9797
## ckd 0.9924
## heart_failure 0.9913
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13952.2 on 19485 degrees of freedom
## Residual deviance: 9312.9 on 19454 degrees of freedom
## AIC: 9376.9
##
## Number of Fisher Scoring iterations: 17# Odds ratios with Wald 95% CIs
or_tbl <- tibble(
term = names(coef(logit_fit)),
OR = exp(coef(logit_fit)),
lo = exp(confint.default(logit_fit)[, 1]),
hi = exp(confint.default(logit_fit)[, 2]),
p = summary(logit_fit)$coefficients[, 4]
) %>%
filter(term != "(Intercept)") %>%
mutate(across(c(OR, lo, hi), ~ round(.x, 2)), p = signif(p, 3))
or_tblggplot(or_tbl, aes(x = OR, y = reorder(term, OR))) +
geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.2, colour = "#3B6EA8") +
geom_point(colour = "#3B6EA8", size = 2) +
scale_x_log10() +
labs(x = "Odds ratio (log scale, 95% CI)", y = NULL,
title = "Logistic regression: odds of a documented complication") +
theme_minimal()
3. Random forest
rf_fit <- fit_rf(train, features, seed = 719, num.trees = 500)
rf_fit## Ranger result
##
## Call:
## ranger(formula = as.formula(paste("outcome_f ~", paste(features, collapse = " + "))), data = train, num.trees = num.trees, probability = TRUE, importance = "permutation", min.node.size = 25, seed = seed)
##
## Type: Probability estimation
## Number of trees: 500
## Sample size: 19486
## Number of independent variables: 12
## Mtry: 3
## Target node size: 25
## Variable importance mode: permutation
## Splitrule: gini
## OOB prediction error (Brier s.): 0.082026674. Compare on the test set
p_logit <- predict_prob(logit_fit, test)
p_rf <- predict_prob(rf_fit, test)
results <- bind_rows(
`Logistic regression` = classification_metrics(test$any_complication, p_logit, threshold),
`Random forest` = classification_metrics(test$any_complication, p_rf, threshold),
.id = "model"
) %>% mutate(across(AUC:precision, ~ round(.x, 3)))
results# AUC confidence intervals and a formal comparison (DeLong test)
roc_logit <- roc(test$any_complication, p_logit, quiet = TRUE, direction = "<")
roc_rf <- roc(test$any_complication, p_rf, quiet = TRUE, direction = "<")
ci.auc(roc_logit); ci.auc(roc_rf)## 95% CI: 0.8411-0.8601 (DeLong)## 95% CI: 0.8346-0.854 (DeLong)roc.test(roc_logit, roc_rf, method = "delong")##
## DeLong's test for two correlated ROC curves
##
## data: roc_logit and roc_rf
## Z = 1.5855, p-value = 0.1129
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## -0.00149838 0.01418654
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8506204 0.8442763roc_df <- bind_rows(
tibble(model = "Logistic regression", fpr = 1 - roc_logit$specificities, tpr = roc_logit$sensitivities),
tibble(model = "Random forest", fpr = 1 - roc_rf$specificities, tpr = roc_rf$sensitivities)
)
ggplot(roc_df, aes(fpr, tpr, colour = model)) +
geom_abline(linetype = "dashed", colour = "grey60") +
geom_path(linewidth = 0.9) +
coord_equal() +
scale_colour_manual(values = c("#3B6EA8", "#D98E04")) +
labs(x = "1 - specificity", y = "Sensitivity", colour = NULL,
title = "ROC curves on the test set") +
theme_minimal() + theme(legend.position = "bottom")
# How do sensitivity and specificity move with the threshold?
sapply(c(0.5 * threshold, threshold, 1.5 * threshold), function(t)
unlist(classification_metrics(test$any_complication, p_logit, t)[, c("sensitivity", "specificity", "precision")])) %>%
`colnames<-`(paste0("t=", round(c(0.5, 1, 1.5) * threshold, 3)))## t=0.058 t=0.116 t=0.173
## sensitivity 0.9986684 0.9986684 0.9973369
## specificity 0.6976501 0.6976501 0.6983464
## precision 0.3015682 0.3015682 0.3017728saveRDS(list(train = train, test = test, features = features, threshold = threshold,
logit_fit = logit_fit, rf_fit = rf_fit), "week3_models.rds")5. Model choice
Write your paragraph here. Weigh test-set performance (is the AUC difference meaningful given the CIs and DeLong test?) against interpretability, how easily clinicians could trust and audit the model, and how each model would be maintained.