dx <- load_t2dm(data_path)

1. Condition lookup table

Hypertension was provided. Obesity, hyperlipidemia, CKD and heart failure are added in t2dm_helpers.R; the table is printed here so it appears in the knitted notebook.

condition_lookup
# Sanity check: which actual codes does each pattern capture?
dx %>% filter(!is.na(icd10)) %>% distinct(icd10, DX_NAME) %>%
  crossing(condition_lookup %>% select(condition, pattern)) %>%
  filter(str_detect(icd10, pattern)) %>%
  group_by(condition) %>%
  summarise(n_codes = n(), example_codes = paste(head(sort(icd10), 6), collapse = ", "))

2. Build the table and check the outcome

pt <- build_patient_table(dx)
glimpse(pt)
## Rows: 26,018
## Columns: 21
## $ sid.x            <chr> "S0000103", "S0000535", "S0001681", "S0001858", "S000…
## $ Patient_sex      <fct> Male, Female, Female, Female, Female, Female, Male, F…
## $ Patient_Race     <fct> Black or African American, White, White, Unreported /…
## $ insurance        <fct> Commercial Plan A, Medicaid, Self-pay / Uninsured, Se…
## $ dob_date         <date> 1961-12-26, 1959-10-31, 1977-04-07, 1978-04-03, 1965…
## $ index_dt         <date> 2019-01-09, 2016-03-14, 2018-06-11, 2017-01-31, 2019…
## $ age_at_index     <dbl> 57, 56, 41, 38, 54, 48, 57, 44, 40, 45, 45, 18, 58, 4…
## $ n_encounters     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ n_specialties    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ n_dx             <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ n_other_dx       <int> 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1,…
## $ any_telemed      <int> 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0,…
## $ portal_active    <int> 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0,…
## $ main_specialty   <fct> Primary Care, Primary Care, Other specialty, Primary …
## $ hypertension     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ obesity          <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ hyperlipidemia   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ ckd              <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ heart_failure    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ any_complication <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0,…
## $ outcome_f        <fct> no, no, no, no, no, no, no, no, no, no, no, no, yes, …
# One row per patient?
stopifnot(nrow(pt) == n_distinct(dx$sid.x), !any(duplicated(pt$sid.x)))
tibble(rows = nrow(pt), unique_patients = n_distinct(pt$sid.x))
# Which codes define the outcome?
dx %>% filter(str_detect(icd10, outcome_pattern)) %>%
  count(icd10, DX_NAME, sort = TRUE) %>% head(20)
# Outcome prevalence
pt %>% count(any_complication) %>% mutate(pct = round(100 * n / sum(n), 1))
# Cross-check with an independent calculation straight from the long data
n_check <- dx %>% filter(str_detect(icd10, outcome_pattern)) %>% distinct(sid.x) %>% nrow()
stopifnot(n_check == sum(pt$any_complication))
n_check
## [1] 3007

3. Leakage check

# 1. No feature column is derived from E11 codes.
#    Condition patterns must not match any E11 code:
e11_codes <- unique(dx$icd10[str_detect(dx$icd10, e11_pattern) & !is.na(dx$icd10)])
leaks <- sapply(condition_lookup$pattern, function(p) any(str_detect(e11_codes, p)))
stopifnot(!any(leaks))

# 2. n_other_dx excludes E11 codes; n_dx does not.
pt %>% summarise(mean_n_dx = mean(n_dx), mean_n_other_dx = mean(n_other_dx))
# n_dx is mechanically higher for outcome patients (their complication codes are
# counted in it). Compare the two counts by outcome:
pt %>% group_by(any_complication) %>%
  summarise(mean_n_dx = mean(n_dx), mean_n_other_dx = mean(n_other_dx))
# 3. Confirm the model feature list does not include n_dx or the outcome
stopifnot(!any(c("n_dx", "any_complication", "outcome_f") %in% model_features))
model_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"

Note for discussion: in real coding practice, diabetic CKD (E11.22) is coded with an additional N18 code, so the ckd flag can partly encode the outcome even though it is not an E11 code. Check how often they co-occur:

dx_flags <- dx %>% filter(!is.na(icd10)) %>% group_by(sid.x) %>%
  summarise(has_e1122 = any(str_detect(icd10, "^E11\\.22")),
            has_n18   = any(str_detect(icd10, "^N18")), .groups = "drop")
dx_flags %>% count(has_e1122, has_n18)

4. Optional stretch features

# n_specialties : patients seen by more specialties may have more complex disease.
# any_telemed   : telemedicine use may reflect access barriers (near-constant, so expect little signal).
# portal_active : portal use is a proxy for engagement with care.
pt %>% select(all_of(stretch_features)) %>% summary()
##  n_specialties    any_telemed     portal_active    
##  Min.   :1.000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:1.000   1st Qu.:0.0000   1st Qu.:0.00000  
##  Median :1.000   Median :1.0000   Median :0.00000  
##  Mean   :1.001   Mean   :0.5045   Mean   :0.08682  
##  3rd Qu.:1.000   3rd Qu.:1.0000   3rd Qu.:0.00000  
##  Max.   :2.000   Max.   :1.0000   Max.   :1.00000

5. Label and feature summary

num_vars <- c("age_at_index", "n_encounters", "n_other_dx", "n_specialties")
bin_vars <- c(condition_lookup$condition, "any_telemed", "portal_active")

numeric_summary <- pt %>%
  group_by(any_complication) %>%
  summarise(across(all_of(num_vars), list(mean = ~ mean(.x, na.rm = TRUE),
                                          median = ~ median(.x, na.rm = TRUE)))) %>%
  pivot_longer(-any_complication) %>%
  pivot_wider(names_from = any_complication, values_from = value, names_prefix = "outcome_")
numeric_summary
binary_summary <- pt %>%
  group_by(any_complication) %>%
  summarise(across(all_of(bin_vars), ~ round(100 * mean(.x), 1))) %>%
  pivot_longer(-any_complication, names_to = "feature", values_to = "pct_with_feature") %>%
  pivot_wider(names_from = any_complication, values_from = pct_with_feature, names_prefix = "outcome_")
binary_summary
categorical_summary <- bind_rows(
  pt %>% count(level = Patient_sex)    %>% mutate(variable = "Patient_sex"),
  pt %>% count(level = Patient_Race)   %>% mutate(variable = "Patient_Race"),
  pt %>% count(level = insurance)      %>% mutate(variable = "insurance"),
  pt %>% count(level = main_specialty) %>% mutate(variable = "main_specialty")
) %>%
  left_join(
    bind_rows(
      pt %>% group_by(level = Patient_sex)    %>% summarise(pct_outcome = 100 * mean(any_complication)) %>% mutate(variable = "Patient_sex"),
      pt %>% group_by(level = Patient_Race)   %>% summarise(pct_outcome = 100 * mean(any_complication)) %>% mutate(variable = "Patient_Race"),
      pt %>% group_by(level = insurance)      %>% summarise(pct_outcome = 100 * mean(any_complication)) %>% mutate(variable = "insurance"),
      pt %>% group_by(level = main_specialty) %>% summarise(pct_outcome = 100 * mean(any_complication)) %>% mutate(variable = "main_specialty")
    ), by = c("variable", "level")) %>%
  mutate(pct_outcome = round(pct_outcome, 1)) %>%
  select(variable, level, n, pct_outcome)
categorical_summary
# Missing values in the model features
pt %>% summarise(across(all_of(model_features), ~ sum(is.na(.x)))) %>%
  pivot_longer(everything(), names_to = "feature", values_to = "n_missing") %>%
  filter(n_missing > 0)
saveRDS(pt, "patient_table.rds")