DM FIM Program Outcomes
2026-09-24
============================================================
SECTION 1: DATA LOAD & COLUMN CONFIG
Edit column names here; nowhere else.
============================================================
df_raw <- read_excel("DMFIM Dataframe and Roster PHI Redacted.xlsx",
col_types = "text", sheet = "DF") %>%
clean_names()
stopifnot("enrolled" %in% names(df_raw))
stopifnot("cohort" %in% names(df_raw))
# The "DF" sheet carries Excel dropdown/data-validation formatting far beyond
# the real data (read_excel returns ~1,048,487 rows, almost all blank).
# Keep only rows with a meaningful amount of data actually entered.
n_before <- nrow(df_raw)
df <- df_raw[rowSums(!is.na(df_raw)) > 5, ]
cat("Loaded", n_before, "raw rows from the sheet;", nrow(df), "are real patient rows.\n")## Loaded 1048487 raw rows from the sheet; 1406 are real patient rows.# ---- Chart-review outcome columns (pre-built change scores; validated against
# absolute values below and used only for cross-checks, not directly) ----
a1c_t1_col <- "a1c_change_t1"
ed_t1_col <- "ed_utilization_t1"
ip_t1_col <- "ip_utilization_t1"
a1c_t2_col <- "a1c_change_t2"
ed_t2_col <- "ed_utilization_t2"
hosp_t2_col <- "hospital_utilization_t2"
# ---- HbA1c value columns (absolute values — source of truth for A1C change) ----
# Baseline is patient-relative (chart-pulled from each patient's own 3 months
# pre-enrollment) and is the only one of these anchored to the individual
# patient's timeline. T1/T2/9mo/12mo are chart-review batch pulls, not
# calendar-verified follow-up windows - treat them as sequential checkpoints.
col_a1c_all <- "hb_a1c_value_last_3_months" # baseline
col_a1c_t1 <- "hemoglobin_a1c_11_20_2_20_t1" # T1
col_a1c_t2 <- "hemoglobin_a1c_t2" # T2
col_a1c_9mo <- "hemoglobin_a1c_9months_post" # 9 months post
col_a1c_12mo <- "hemoglobin_a1c_12_months_post" # 12 months post
# ---- Utilization columns (ED and IP, same baseline/T1/T2/9mo/12mo checkpoints
# as A1C above - each wave's utilization is chart-pulled alongside that
# wave's A1C draw) ----
col_ed_base <- "ed_vis_last_90_days"
col_ip_base <- "hosp_last_90_days"
col_ed_9mo <- "ed_utilization_9months_post"
col_ip_9mo <- "hospital_utilization_9months_post"
col_ed_12mo <- "ed_utilization_12_months_post"
col_ip_12mo <- "hospital_utilization_12_months_post"
# ---- Clinical risk ----
col_lace <- "lace_readmission_score_score_column"
# ---- Food insecurity (Hunger Vital Sign, 3-item; positive = any "Yes") ----
col_ps_worry_t0 <- "ps_in_the_past_12_months_were_you_ever_worried_your_food_would_run_out_before_you_had_money_to_buy_more"
col_ps_lastmonth_t0 <- "ps_in_the_past_12_months_did_the_food_you_could_afford_not_last_until_the_end_of_the_month_and_you_couldn_t_get_more"
col_ps_skip_t0 <- "ps_in_the_last_month_did_anyone_in_your_household_have_to_skip_meals_due_to_lack_of_food"
col_worry_t2 <- "t2_in_the_past_6_months_were_you_ever_worried_your_food_would_run_out_before_you_had_money_to_buy_more"
col_lastmonth_t2 <- "t2_in_the_past_6_months_did_the_food_you_could_afford_not_last_until_the_end_of_the_month_and_you_couldn_t_get_more"
col_skip_t2 <- "t2_in_the_last_month_did_anyone_in_your_household_have_to_skip_meals_due_to_lack_of_food"
# ---- Patient experience survey columns ----
col_help_choices <- "t1_did_the_nutrition_material_you_got_from_white_plains_hospital_help_you_make_healthier_food_choices"
col_help_choices_t2 <- "t2_did_the_nutrition_material_you_got_from_white_plains_hospital_help_you_make_healthier_food_choices"
col_help_control <- "t1_how_much_did_this_program_help_you_feel_more_in_control_of_your_diabetes"
col_help_control_t2 <- "t2_how_much_did_this_program_help_you_feel_more_in_control_of_your_diabetes"
col_ease_t1 <- "t1_in_the_past_week_how_easy_was_it_to_use_the_food_in_your_meals"
col_ease_t2 <- "t2_in_the_past_week_how_easy_was_it_to_use_the_food_in_your_meals"
col_meals_t1 <- "t1_in_the_past_week_how_many_meals_did_you_make_using_food_from_the_program"
col_meals_t2 <- "t2_in_the_past_week_how_many_meals_did_you_make_using_food_from_the_program"
col_waste_t1 <- "t1_did_you_throw_away_any_of_the_food_from_your_most_recent_delivery"
col_waste_t2 <- "t2_did_you_throw_away_any_of_the_food_from_your_most_recent_delivery"
# ---- Confidence survey columns (T1/T2, 3 items each; T0 only has 2 of the
# 3 - meal prep confidence was never asked at baseline) ----
col_conf_food_t0 <- "t0_are_you_confident_in_making_food_choices_that_help_control_your_blood_sugar"
col_conf_labels_t0 <- "t0_are_you_confident_in_understanding_how_to_read_food_labels_and_nutrition_information"
col_conf_food_t1 <- "t1_how_confident_are_you_in_making_food_choices_that_help_control_your_blood_sugar"
col_conf_labels_t1 <- "t1_how_confident_are_you_in_understanding_food_labels_or_nutrition_information"
col_conf_prepare_t1 <- "t1_how_confident_are_you_in_your_ability_to_prepare_healthy_meals_with_the_food_you_typically_have_at_home"
col_conf_food_t2 <- "t2_how_confident_are_you_in_making_food_choices_that_help_control_your_blood_sugar"
col_conf_labels_t2 <- "t2_how_confident_are_you_in_understanding_food_labels_or_nutrition_information"
col_conf_prepare_t2 <- "t2_how_confident_are_you_in_your_ability_to_prepare_healthy_meals_with_the_foods_you_typically_have_at_home"
# ---- Demographic / socioeconomic columns ----
race_col <- "race"
eth_col <- "ethnicity"
sex_col <- "legal_sex"
ins_col <- "primary_cvg"
# ---- Enrollment funnel columns (eligibility screen -> food-insecurity
# prescreen -> enrollment) ----
col_zip <- "zip_code"
col_t0_contact <- "t0_contact" # outreach/contact outcome prior to enrollment
# Westchester County ZIP codes (delivery area after PDSA 2 expanded from a
# 15-mile radius to the full county) - used only for the eligibility funnel
# below, not for any outcome analysis.
westchester_zips <- c(
"10501","10502","10503","10504","10505","10506","10507","10508",
"10509","10510","10511","10512","10514","10516","10517","10518",
"10519","10520","10521","10522","10523","10524","10526","10527",
"10528","10530","10532","10533","10535","10536","10537","10538",
"10540","10541","10542","10543","10545","10546","10547","10548",
"10549","10550","10551","10552","10553","10560","10562","10564",
"10566","10567","10570","10571","10572","10573","10576","10577",
"10578","10579","10580","10583","10587","10588","10589","10590",
"10591","10594","10595","10596","10597","10598","10601","10602",
"10603","10604","10605","10606","10607","10701","10702","10703",
"10704","10705","10706","10707","10708","10709","10710","10801",
"10802","10803","10804","10805","10901","10920","10954","10960",
"10965","10970","10976","10977","10980","10982","10983","10984",
"10986","10987","10989","10992","10993","10994","10996"
)
# ---- Race code map ----
RACE_MAP <- c(
R1 = "American Indian / Alaska Native",
R2 = "Asian",
R3 = "Black or African American",
R4 = "Native Hawaiian / Pacific Islander",
R5 = "White",
R9 = "Other"
)
# ---- Color palette ----
COL_IMPROVED <- "#2E7D32"
COL_NOCHANGE <- "#757575"
COL_WORSENED <- "#C62828"
COL_ENROLLED <- "#243B5A"
COL_ALL <- "#CFE1F2"
COL_TEAL <- "#20A39E"
COL_BLUE <- "#1F4E79"
COL_MID_BLUE <- "#2E86C1"============================================================
SECTION 2: HELPER FUNCTIONS (defined once)
============================================================
is_nonblank <- function(x) {
x_std <- str_to_upper(str_squish(as.character(x)))
!(is.na(x_std) | x_std == "" | x_std %in% c("NA", "N/A", "NULL"))
}
parse_count <- function(x) {
s <- str_squish(as.character(x))
s[s == ""] <- NA_character_
s[str_to_upper(s) %in% c("NA", "N/A", "NULL")] <- NA_character_
suppressWarnings(as.numeric(str_replace_all(s, "[^0-9.-]", "")))
}
clean_a1c <- function(x) {
x <- str_squish(as.character(x))
x[x %in% c("", "NA", "N/A", "na", "n/a", "NULL", "null")] <- NA_character_
x[x == "<4.0"] <- "4.0"
x[x == ">14.0"] <- "14.0"
x[x == ">20.0"] <- "20.0"
x <- gsub("%", "", x)
suppressWarnings(as.numeric(x))
}
clean_yesno <- function(x) {
s <- str_to_upper(str_squish(as.character(x)))
case_when(
s %in% c("YES", "Y") ~ TRUE,
s %in% c("NO", "N") ~ FALSE,
TRUE ~ NA
)
}
clean_impact <- function(x) {
s <- str_to_upper(str_squish(as.character(x)))
case_when(
is.na(s) ~ NA_character_,
s %in% c("", "NA", "N/A", "NULL") ~ NA_character_,
s %in% c("5") ~ "A lot",
s %in% c("4") ~ "Somewhat",
s %in% c("3", "2") ~ "A little",
s %in% c("1") ~ "Not at all",
str_detect(s, "^A LOT") ~ "A lot",
str_detect(s, "^SOMEWHAT") ~ "Somewhat",
str_detect(s, "A LITTLE") ~ "A little",
str_detect(s, "NOT AT ALL") ~ "Not at all",
s %in% c("YES", "Y") ~ "Yes",
s %in% c("NO", "N") ~ "No",
TRUE ~ "Other / unclear"
)
}
clean_ease <- function(x) {
s <- str_to_upper(str_squish(as.character(x)))
case_when(
is.na(s) ~ NA_character_,
s %in% c("", "NA", "N/A", "NULL") ~ NA_character_,
s == "1" ~ "Very easy",
s == "2" ~ "Easy",
s == "3" ~ "Neutral",
s == "4" ~ "Difficult",
str_detect(s, "VERY EASY") ~ "Very easy",
# "Somehwat"/"Somewhat" difficult - catches the observed typo directly,
# must run before the general DIFFICULT check below
str_detect(s, "SOME.HWAT DIFFICULT|SOMEWHAT DIFFICULT") ~ "Somewhat difficult",
str_detect(s, "^EASY") ~ "Easy",
str_detect(s, "NEUT.?L|MODERATELY") ~ "Neutral", # covers the "Neutal" typo
str_detect(s, "DIFFICULT") ~ "Difficult",
TRUE ~ NA_character_
)
}
# Ordinal score for clean_ease(), for paired before/after comparison.
# Originally only Neutral/Easy/Very easy were recognized at all - Difficult
# and Somewhat difficult responses (incl. typos like "Somehwat"/"Neutal")
# were silently dropped as NA, inflating the apparent % easy. Fixed above;
# scored here on a full 5-point ordinal scale.
ease_score <- function(x) {
v <- clean_ease(x)
case_when(
v == "Difficult" ~ 1, v == "Somewhat difficult" ~ 2, v == "Neutral" ~ 3,
v == "Easy" ~ 4, v == "Very easy" ~ 5, TRUE ~ NA_real_
)
}
# "In the past week, how many meals did you make using food from the
# program?" - free text, not a clean number: digits ("3"), spelled-out
# ("three"), ranges ("four to six", "2-3x"), "N+"/"N or more", and "Nx"
# (N times). Ambiguous entries (different units like "3 a day", "unsure of
# #", "all") are excluded rather than guessed.
clean_meals_numeric <- function(x) {
s <- str_to_lower(str_squish(as.character(x)))
s[s %in% c("", "na", "n/a", "null")] <- NA
s[str_detect(s, "unsure|unk\\b")] <- NA
s[str_detect(s, "^all$")] <- NA
s[str_detect(s, "a day")] <- NA # different unit (per day, not per week)
s <- ifelse(str_detect(s, "^o\\b"), str_replace(s, "^o\\b", "0"), s) # "O" typo for "0"
s <- str_replace(s, "to(\\d)", "to \\1") # "TO1" -> "to 1"
word_map <- c(zero = "0", one = "1", two = "2", three = "3", four = "4",
five = "5", six = "6", seven = "7", eight = "8", nine = "9", ten = "10")
for (w in names(word_map)) s <- str_replace_all(s, paste0("\\b", w, "\\b"), word_map[[w]])
s[str_detect(s, "ran out")] <- "0"
rng <- str_match(s, "(\\d+)\\s*(?:to|-)\\s*(\\d+)")
is_range <- !is.na(rng[, 1])
out <- suppressWarnings(as.numeric(s))
out[is_range] <- (as.numeric(rng[is_range, 2]) + as.numeric(rng[is_range, 3])) / 2
plus <- str_match(s, "(\\d+)\\s*(?:\\+|or more|more)")
is_plus <- !is.na(plus[, 1]) & is.na(out)
out[is_plus] <- as.numeric(plus[is_plus, 2])
xsuffix <- str_match(s, "^(\\d+)x$")
is_x <- !is.na(xsuffix[, 1]) & is.na(out)
out[is_x] <- as.numeric(xsuffix[is_x, 2])
out
}
# 1 (not confident) - 4 (very confident)
clean_conf_score <- function(x) {
s <- str_to_upper(str_squish(as.character(x)))
case_when(
is.na(s) ~ NA_real_,
s %in% c("", "NA", "N/A", "NULL", "NOT SURVEYED", "NOT ASKED") ~ NA_real_,
s == "1" ~ 1, s == "2" ~ 2, s == "3" ~ 3, s %in% c("4", "5") ~ 4,
str_detect(s, "VERY") & str_detect(s, "CONF") ~ 4,
str_detect(s, "SOMEWHAT") & str_detect(s, "NOT") & str_detect(s, "CONF") ~ 2,
str_detect(s, "NOT") & str_detect(s, "CONF") ~ 1,
str_detect(s, "SOMEWHAT") & str_detect(s, "CONF") ~ 3,
str_detect(s, "CONF") ~ 3,
TRUE ~ NA_real_
)
}
# High/Low confidence bucket from a clean 1-4 score (>=3 = Confident/Very
# confident = "High"; <3 = Not/Somewhat not confident = "Low"). A one-step
# move within a bucket (e.g. Very confident -> Confident) is NOT a bucket
# change - only crossing this line counts as improved/worsened.
conf_bucket <- function(score) if_else(score >= 3, "High", "Low")
# T0's confidence responses use a messier, inconsistent format (plain 0-4
# numbers, Yes/No, free text) that doesn't line up cleanly with T1/T2's 1-4
# scale. Bucketed directly rather than scored numerically: 0-1 = Low,
# 2 = Unsure (excluded from paired comparisons, not forced into a
# direction), 3-4 = High.
conf_bucket_t0 <- function(x) {
s <- str_to_upper(str_squish(as.character(x)))
case_when(
is.na(s) | s %in% c("", "NA", "N/A", "NULL", "MAYBE", "M") ~ NA_character_,
s %in% c("0", "1") ~ "Low",
s == "2" ~ "Unsure",
s %in% c("3", "4") ~ "High",
str_detect(s, "NOT AT ALL") ~ "Low",
str_detect(s, "SOMEWHAT") & str_detect(s, "NOT") ~ "Low",
str_detect(s, "NEUTRAL") ~ "Unsure",
str_detect(s, "NOT") & str_detect(s, "CONF") ~ "Low",
str_detect(s, "VERY") & str_detect(s, "CONF") ~ "High",
str_detect(s, "MOSTLY") & str_detect(s, "CONF") ~ "High",
str_detect(s, "SOMEWHAT") & str_detect(s, "CONF") ~ "High",
s %in% c("YES", "Y") ~ "High",
s %in% c("NO", "N") ~ "Low",
str_detect(s, "CONF") ~ "High",
TRUE ~ NA_character_
)
}
collapse_insurance <- function(x) {
s <- str_to_upper(str_squish(as.character(x)))
s[s %in% c("", "NA", "N/A", "NULL")] <- NA_character_
case_when(
is.na(s) ~ NA_character_,
str_detect(s, "SELF PAY|UNINSURED|WORKERS COMP|NO FAULT") ~ "Uninsured / Self Pay",
str_detect(s, "DUAL|SNP|ARCHCARE SENIOR LIFE|PACE") ~ "Dual Eligible",
str_detect(s, "MEDICARE") ~ "Medicare",
str_detect(s, "MEDICAID|HARP|EMERGENCY MEDICAID|OUT OF STATE MEDICAID|MEDICAID PENDING") ~ "Medicaid",
TRUE ~ "Commercial"
)
}
# ---- Table formatters ----
fmt_mean_sd <- function(x, digits = 2) {
x <- x[!is.na(x)]
if (length(x) == 0) return("–")
sd_str <- if (length(x) < 2) "–" else format(round(sd(x), digits), nsmall = digits)
paste0(format(round(mean(x), digits), nsmall = digits), " (", sd_str, ")")
}
fmt_np <- function(n, denom) {
if (is.na(n) || denom == 0) return("–")
paste0(n, " (", sprintf("%.1f", 100 * n / denom), "%)")
}
# Same as fmt_np but shows the denominator (n/N), used in the outcomes table
fmt_nN <- function(n, denom) {
if (is.na(n) || denom == 0) return("–")
paste0(n, "/", denom, " (", sprintf("%.1f", 100 * n / denom), "%)")
}
fmt_p <- function(p) {
if (is.na(p)) return("–")
if (p < 0.001) return("<0.001")
format(round(p, 3), nsmall = 3)
}
# p-value formatted with a trailing asterisk when significant at .05
fmt_p_star <- function(p) {
s <- fmt_p(p)
if (!is.na(p) && p < 0.05) paste0(s, "*") else s
}
# "mean (SD), n=X" - same as fmt_mean_sd but with the sample size shown inline
fmt_mean_sd_n <- function(x, digits = 2) {
x <- x[!is.na(x)]
if (length(x) == 0) return("–")
sd_str <- if (length(x) < 2) "–" else format(round(sd(x), digits), nsmall = digits)
paste0(format(round(mean(x), digits), nsmall = digits), " (", sd_str, "), n=", length(x))
}
# Categorical cell for Table 1: "n (pct%)" of non-missing in that column
cat_stat <- function(x, level) {
denom <- sum(!is.na(x))
if (denom == 0) return("–")
n <- sum(x == level, na.rm = TRUE)
fmt_np(n, denom)
}
# p-values for a continuous variable: one-way ANOVA across Cohorts 1-5 and
# Not-enrolled, followed by Dunnett's test (each cohort vs. the not-enrolled
# reference, single-step multiplicity-adjusted). This is the standard method
# for "which of several groups differs from one common reference group" -
# more appropriate here than 5 separate unadjusted pairwise tests.
COHORT_LEVELS <- c("Cohort 1", "Cohort 2", "Cohort 3", "Cohort 4", "Cohort 5")
dunnett_pvals <- function(x, group) {
out <- setNames(rep(NA_real_, length(COHORT_LEVELS)), COHORT_LEVELS)
d <- data.frame(x = x, group = group, stringsAsFactors = FALSE)
d <- d[!is.na(d$x) & !is.na(d$group), ]
d$group <- factor(d$group)
if (!"Not enrolled" %in% levels(d$group) || nlevels(droplevels(d$group)) < 2) return(out)
d$group <- stats::relevel(d$group, ref = "Not enrolled")
result <- tryCatch({
mod <- stats::aov(x ~ group, data = d)
dt <- multcomp::glht(mod, linfct = multcomp::mcp(group = "Dunnett"))
summary(dt)
}, error = function(e) NULL)
if (is.null(result)) return(out)
nm <- names(result$test$coefficients)
pv <- as.numeric(result$test$pvalues)
for (i in seq_along(nm)) {
cohort_name <- trimws(strsplit(nm[i], " - ")[[1]][1])
if (cohort_name %in% names(out)) out[cohort_name] <- pv[i]
}
out
}
# p-value for a categorical variable: chi-square, falling back to simulated
# Fisher's exact test when any expected cell count is small
p_categorical <- function(vals_a, vals_b) {
tab <- table(
c(rep("A", length(vals_a)), rep("B", length(vals_b))),
c(vals_a, vals_b)
)
if (nrow(tab) < 2 || ncol(tab) < 2) return(NA_real_)
chi <- suppressWarnings(chisq.test(tab))
if (any(chi$expected < 5)) {
set.seed(1)
fisher.test(tab, simulate.p.value = TRUE, B = 5000)$p.value
} else {
chi$p.value
}
}
# Each cohort vs. not-enrolled, chi-square/Fisher, Bonferroni-adjusted for
# the N simultaneous comparisons (there's no categorical equivalent of
# Dunnett's test in base R, so a manual multiplicity correction is used
# instead)
cat_pvals_bonf <- function(col, groups_list, ref_data) {
raw <- vapply(COHORT_LEVELS, function(cl) {
p_categorical(groups_list[[cl]][[col]], ref_data[[col]])
}, numeric(1))
setNames(pmin(1, raw * length(COHORT_LEVELS)), COHORT_LEVELS) # pmin() drops names - restore them
}
# ---- Minimal raw-HTML table builder ----
# Used instead of kable() %>% kable_styling() for tables that embed raw HTML
# inside cells (e.g. a bold/asterisked p-value on its own line within a
# value's own cell): kableExtra's kable_styling()/row_spec()/pack_rows() all
# re-escape already-inserted raw HTML in cell content, corrupting it. This
# builder never round-trips through kableExtra, so nothing gets re-escaped.
render_html_table <- function(df, caption, bold_rows = integer(0),
sections = list(), footnote = "") {
nc <- ncol(df)
th <- paste0(
"<th style='background:#1A2B4A;color:white;padding:6px 10px;text-align:",
ifelse(seq_len(nc) == 1, "left", "center"), ";'>", names(df), "</th>",
collapse = "")
body <- character(0)
for (i in seq_len(nrow(df))) {
for (sname in names(sections)) {
if (isTRUE(sections[[sname]] == i)) {
body <- c(body, paste0(
"<tr><td colspan='", nc, "' style='background:#F0F4F8;color:#1A2B4A;",
"font-weight:bold;padding:6px 10px;'>", sname, "</td></tr>"))
}
}
row_bold <- if (i %in% bold_rows) "font-weight:bold;" else ""
cells <- vapply(seq_len(nc), function(j) {
align <- if (j == 1) "left" else "center"
paste0("<td style='padding:4px 10px;text-align:", align, ";", row_bold, "'>",
df[i, j], "</td>")
}, character(1))
body <- c(body, paste0("<tr>", paste(cells, collapse = ""), "</tr>"))
}
htmltools::HTML(paste0(
"<div style='overflow-x:auto;'>",
"<div style='font-weight:bold;font-size:1.15em;margin-bottom:6px;'>", caption, "</div>",
"<table style='border-collapse:collapse;width:100%;font-size:13px;'>",
"<thead><tr>", th, "</tr></thead>",
"<tbody>", paste(body, collapse = ""), "</tbody>",
"</table>",
"<div style='font-size:11px;color:#555;margin-top:6px;'><em>Notes:</em> ", footnote, "</div>",
"</div>"
))
}============================================================
SECTION 3: DATA PREPARATION
============================================================
# ---- Disenrollment signal detection (T1/T2 chart-review notes) ----
# Text-pattern screen over t1_notes/t2_notes. This flags CANDIDATES for
# disenrollment/attrition, not a confirmed clinical determination - the
# project's own convention (see prior case-by-case review) is that these
# need human review before being treated as final. "Did not receive
# educational materials" is explicitly NOT treated as a disenrollment signal.
disenroll_pattern <- paste0(
"DECEAS|DIED|PASSED|HOSPICE|",
"CANCEL.*DELIVER|HOLD ON THE DELIVERIES|",
"WITHDREW|WITHDRAW|REQUESTED.*REMOVED|REQUESTED TO BE REMOVED|",
"NEVER RECEIV|NO DELIVERIES|DID NOT RECEIV|HAS NOT RECEIVED|",
"NEVER RECEIVED|NEVER GOTTEN|NEVER GOT|",
"ONLY RECEIVED FIRST|PAST THREE MONTHS|PAST 3M|LAST 3M|",
"NO LONGER IN SERVICE|OUT OF SERVICE|NUMBER.*NO LONGER|",
"NEW ADDRESS|WRONG ADDRESS|SHE MOVED|HE MOVED|PT MOVED|DIDNT TELL|",
"NOT RECEPTIVE|DOES NOT RECALL ENROLLING|DOESN.T RECALL ENROLLING|",
"HOSPITALIZED|NURSING HOME|\\bNH\\b|REHAB SINCE"
)
disenroll_exclude_pattern <- "DID NOT RECEIVE EDUCATIONAL|DID NOT RECIEVE EDUCATIONAL|EDUCATIONAL MATERIALS"
classify_disenroll_signal <- function(t1_raw, t2_raw) {
combined <- str_to_upper(paste(t1_raw, t2_raw))
case_when(
str_detect(combined, "DECEAS|DIED|PASSED|HOSPICE") ~ "Deceased",
str_detect(combined, "REQUESTED.*REMOVED|REQUESTED TO BE REMOVED") ~ "Requested removal",
str_detect(combined, "NURSING HOME|\\bNH\\b|REHAB") ~ "Hospitalized / facility",
str_detect(combined, "CANCEL.*DELIVER|HOLD ON THE DELIVER") ~ "Delivery cancelled / hold",
str_detect(combined, "NEW ADDRESS|WRONG ADDRESS|MOVED|DIDNT TELL") ~ "Address issue",
str_detect(combined, "NOT RECEPTIVE|DOES NOT RECALL|DOESN.T RECALL") ~ "Unresponsive / disengaged",
str_detect(combined, "NEVER RECEIV|NO DELIVERIES|DID NOT RECEIV|HAS NOT RECEIVED|NEVER RECEIVED|NEVER GOTTEN|NEVER GOT|ONLY RECEIVED FIRST|PAST THREE MONTHS|PAST 3M|LAST 3M") ~ "Never/stopped receiving deliveries",
TRUE ~ "Other"
)
}df_status <- df %>%
mutate(
enrolled_clean = str_to_upper(str_squish(as.character(enrolled))),
cohort_num = suppressWarnings(as.integer(str_extract(str_squish(as.character(cohort)), "\\d+"))),
enrolled_group = if_else(enrolled_clean == "Y", "Enrolled", "Not enrolled"),
cohort_final = case_when(
enrolled_group == "Enrolled" & cohort_num == 1 ~ "Cohort 1",
enrolled_group == "Enrolled" & cohort_num == 2 ~ "Cohort 2",
enrolled_group == "Enrolled" & cohort_num == 3 ~ "Cohort 3",
enrolled_group == "Enrolled" & cohort_num == 4 ~ "Cohort 4",
enrolled_group == "Enrolled" & cohort_num == 5 ~ "Cohort 5",
TRUE ~ "Not enrolled"
)
)
# ---- Demographic / socioeconomic cleaning ----
df_status <- df_status %>%
mutate(
age_num = suppressWarnings(as.numeric(age)),
sex_group = case_when(
str_to_upper(str_squish(.data[[sex_col]])) == "MALE" ~ "Male",
str_to_upper(str_squish(.data[[sex_col]])) == "FEMALE" ~ "Female",
TRUE ~ NA_character_
),
race_clean = {
# NOTE: do not str_squish() before splitting - squish collapses the
# embedded \r\n that separates multi-race entries into a single space,
# which silently drops every code after the first. Trim edges only,
# split on the real separators, then clean each item.
raw <- as.character(.data[[race_col]])
sapply(raw, function(s) {
if (is.na(s)) return(NA_character_)
s <- str_trim(s)
if (s == "") return(NA_character_)
items <- str_trim(unlist(str_split(s, "\\r?\\n|\\s*[;,]\\s*")))
items <- items[items != ""]
codes <- str_extract(items, "^[A-Za-z]\\d+")
recoded <- ifelse(!is.na(codes) & codes %in% names(RACE_MAP), RACE_MAP[codes], items)
paste(sort(unique(recoded)), collapse = "; ")
}, USE.NAMES = FALSE)
},
race_group = case_when(
str_detect(race_clean, "White") ~ "White",
str_detect(race_clean, "Black") ~ "Black",
!is.na(race_clean) ~ "Other",
TRUE ~ NA_character_
),
ethnicity_clean = {
raw <- str_to_upper(str_squish(as.character(.data[[eth_col]])))
case_when(
str_detect(raw, "E1") ~ "Spanish / Hispanic / Latino",
str_detect(raw, "E2") ~ "Not Spanish / Hispanic / Latino",
str_detect(raw, "S1") ~ "Patient Declined",
str_detect(raw, "S2") ~ "Patient Unavailable",
TRUE ~ NA_character_
)
},
# Declined and Unavailable are collapsed for display/testing - both are
# tiny (n=48 and n=3) and represent the same "unknown" status
ethnicity_display = if_else(ethnicity_clean %in% c("Patient Declined", "Patient Unavailable"),
"Declined / Unavailable", ethnicity_clean),
insurance_group = collapse_insurance(.data[[ins_col]]),
# ---- Clinical: A1C at every checkpoint ----
a1c_base = clean_a1c(.data[[col_a1c_all]]),
a1c_t1 = clean_a1c(.data[[col_a1c_t1]]),
a1c_t2 = clean_a1c(.data[[col_a1c_t2]]),
a1c_9mo = clean_a1c(.data[[col_a1c_9mo]]),
a1c_12mo = clean_a1c(.data[[col_a1c_12mo]]),
lace_num = parse_count(.data[[col_lace]]),
# ---- Utilization at every checkpoint (ED and IP, same waves as A1C) ----
ed_base = parse_count(.data[[col_ed_base]]),
ip_base = parse_count(.data[[col_ip_base]]),
ed_t1 = parse_count(.data[[ed_t1_col]]),
ip_t1 = parse_count(.data[[ip_t1_col]]),
ed_t2 = parse_count(.data[[ed_t2_col]]),
ip_t2 = parse_count(.data[[hosp_t2_col]]),
ed_9mo = parse_count(.data[[col_ed_9mo]]),
ip_9mo = parse_count(.data[[col_ip_9mo]]),
ed_12mo = parse_count(.data[[col_ed_12mo]]),
ip_12mo = parse_count(.data[[col_ip_12mo]]),
# ---- A1C change: computed from absolute values, not the pre-built delta
# columns. QA below shows why. ----
a1c_change_t1_calc = a1c_t1 - a1c_base,
a1c_change_t2_calc = a1c_t2 - a1c_base,
# ---- Food insecurity (Hunger Vital Sign, positive = any "Yes") ----
fi_t0_pos = clean_yesno(.data[[col_ps_worry_t0]]) | clean_yesno(.data[[col_ps_lastmonth_t0]]) | clean_yesno(.data[[col_ps_skip_t0]]),
fi_t0_any = is_nonblank(.data[[col_ps_worry_t0]]) | is_nonblank(.data[[col_ps_lastmonth_t0]]) | is_nonblank(.data[[col_ps_skip_t0]]),
fi_t2_pos = clean_yesno(.data[[col_worry_t2]]) | clean_yesno(.data[[col_lastmonth_t2]]) | clean_yesno(.data[[col_skip_t2]]),
fi_t2_any = is_nonblank(.data[[col_worry_t2]]) | is_nonblank(.data[[col_lastmonth_t2]]) | is_nonblank(.data[[col_skip_t2]]),
# ---- Food ease, T1/T2 (ordinal 1-5, for paired change) ----
ease_t1_score = ease_score(.data[[col_ease_t1]]),
ease_t2_score = ease_score(.data[[col_ease_t2]]),
# ---- Meals prepared from program food per week, T1/T2 ----
meals_t1 = clean_meals_numeric(.data[[col_meals_t1]]),
meals_t2 = clean_meals_numeric(.data[[col_meals_t2]]),
# ---- Patient confidence: food choices + labels have a baseline (T0)
# item, so their composite is compared baseline -> T2. Meal prep was
# never asked at T0, so it stays T1 -> T2 on its own. ----
conf_food_t0_bucket = conf_bucket_t0(.data[[col_conf_food_t0]]),
conf_labels_t0_bucket = conf_bucket_t0(.data[[col_conf_labels_t0]]),
conf_food_t1_s = clean_conf_score(.data[[col_conf_food_t1]]),
conf_labels_t1_s = clean_conf_score(.data[[col_conf_labels_t1]]),
conf_prepare_t1_s = clean_conf_score(.data[[col_conf_prepare_t1]]),
conf_food_t2_s = clean_conf_score(.data[[col_conf_food_t2]]),
conf_labels_t2_s = clean_conf_score(.data[[col_conf_labels_t2]]),
conf_prepare_t2_s = clean_conf_score(.data[[col_conf_prepare_t2]]),
conf_food_t1_bucket = conf_bucket(conf_food_t1_s),
conf_labels_t1_bucket = conf_bucket(conf_labels_t1_s),
conf_food_t2_bucket = conf_bucket(conf_food_t2_s),
conf_labels_t2_bucket = conf_bucket(conf_labels_t2_s),
conf_prepare_t1_bucket = conf_bucket(conf_prepare_t1_s),
conf_prepare_t2_bucket = conf_bucket(conf_prepare_t2_s),
# Composite (food choices + labels) bucket: "High"/"Low" only when both
# items agree; a mixed pair (one High, one Low/Unsure) is excluded as
# ambiguous rather than averaged across mismatched scales
conf_composite_t0_bucket = case_when(
conf_food_t0_bucket == "High" & conf_labels_t0_bucket == "High" ~ "High",
conf_food_t0_bucket == "Low" & conf_labels_t0_bucket == "Low" ~ "Low",
TRUE ~ NA_character_
),
conf_composite_t2_bucket = case_when(
conf_food_t2_bucket == "High" & conf_labels_t2_bucket == "High" ~ "High",
conf_food_t2_bucket == "Low" & conf_labels_t2_bucket == "Low" ~ "Low",
TRUE ~ NA_character_
),
# ---- Food waste, T1/T2 (TRUE = threw food away) ----
waste_t1 = clean_yesno(.data[[col_waste_t1]]),
waste_t2 = clean_yesno(.data[[col_waste_t2]]),
# ---- Attrition / disenrollment signal (T1 + T2 chart-review notes) ----
t1_notes_raw = str_squish(as.character(t1_notes)),
t2_notes_raw = str_squish(as.character(t2_notes)),
disenroll_hit = (is_nonblank(t1_notes_raw) & str_detect(str_to_upper(t1_notes_raw), disenroll_pattern)) |
(is_nonblank(t2_notes_raw) & str_detect(str_to_upper(t2_notes_raw), disenroll_pattern)),
disenroll_edu_only = str_detect(str_to_upper(coalesce(t1_notes_raw, "")), disenroll_exclude_pattern) &
str_detect(str_to_upper(coalesce(t2_notes_raw, "")), "^NA$|^N/A$|^$") &
!str_detect(str_to_upper(coalesce(t1_notes_raw, "")), disenroll_pattern),
disenroll_flag = disenroll_hit & !disenroll_edu_only,
disenroll_signal = if_else(disenroll_flag, classify_disenroll_signal(t1_notes_raw, t2_notes_raw), NA_character_)
)
# ---- QA: why A1C change is computed from absolute values ----
# The spreadsheet's pre-built a1c_change_t2 column disagrees with
# (T2 value - baseline value) for most of Cohort 2, and for Cohort 3 it
# contains values for patients who have no T2 result at all (turned out to
# be -1 x their T1 value - a copy/formula artifact). Recomputing directly
# from the absolute HbA1c columns avoids both problems.
qa_a1c_t2 <- df_status %>%
filter(cohort_final %in% c("Cohort 1","Cohort 2","Cohort 3","Cohort 4")) %>%
mutate(
delta_precomputed = parse_count(.data[[a1c_t2_col]]),
has_t2_value = !is.na(a1c_t2),
agrees = !is.na(delta_precomputed) & !is.na(a1c_change_t2_calc) &
abs(delta_precomputed - a1c_change_t2_calc) < 0.05
) %>%
group_by(cohort_final) %>%
summarise(
n_with_t2_value = sum(has_t2_value),
n_precomputed_delta = sum(!is.na(delta_precomputed)),
n_precomputed_but_no_t2 = sum(!is.na(delta_precomputed) & !has_t2_value),
n_precomputed_agrees = sum(agrees),
.groups = "drop"
)
cat("QA - pre-built a1c_change_t2 vs recomputed (T2 - baseline):\n")## QA - pre-built a1c_change_t2 vs recomputed (T2 - baseline):print(qa_a1c_t2)## # A tibble: 4 × 5
## cohort_final n_with_t2_value n_precomputed_delta n_precomputed_but_no_t2
## <chr> <int> <int> <int>
## 1 Cohort 1 30 30 1
## 2 Cohort 2 19 17 0
## 3 Cohort 3 0 15 15
## 4 Cohort 4 0 0 0
## # ℹ 1 more variable: n_precomputed_agrees <int># ---- Cohort counts ----
cat("\nCohort counts:\n")##
## Cohort counts:df_status %>% count(cohort_final, name = "N") %>% arrange(cohort_final) %>% print()## # A tibble: 6 × 2
## cohort_final N
## <chr> <int>
## 1 Cohort 1 60
## 2 Cohort 2 40
## 3 Cohort 3 42
## 4 Cohort 4 36
## 5 Cohort 5 9
## 6 Not enrolled 1219# ---- QA FLAG: Cohort 3 has essentially zero T2 data despite being well past
# its graduation window - this looks like a data collection/entry gap, not a
# "not yet due" situation. Every T2 column (chart review AND every T2 survey
# question) is blank for all 42 patients.
qa_c3_t2 <- df_status %>%
filter(cohort_final == "Cohort 3") %>%
mutate(disch = as.Date(suppressWarnings(as.numeric(disch_date_time)), origin = "1899-12-30"),
months_since_disch = round(as.numeric(as.Date("2026-09-08") - disch) / 30.44, 1)) %>%
summarise(
n = n(),
min_months_since_discharge = min(months_since_disch, na.rm = TRUE),
max_months_since_discharge = max(months_since_disch, na.rm = TRUE),
n_with_any_t2_survey_response = sum(is_nonblank(t2_notes) | is_nonblank(t2_contact_y_n) |
is_nonblank(.data[[col_ease_t2]]) | is_nonblank(.data[[col_waste_t2]])),
n_with_t2_a1c = sum(!is.na(a1c_t2)),
n_with_t2_ed_ip = sum(!is.na(ed_t2) | !is.na(ip_t2))
)
cat("\nQA FLAG - Cohort 3 T2 completeness (all patients are well past a ~5-6 month graduation window,\n",
"matching Cohorts 1-2's program length, yet T2 data is essentially absent):\n", sep = "")##
## QA FLAG - Cohort 3 T2 completeness (all patients are well past a ~5-6 month graduation window,
## matching Cohorts 1-2's program length, yet T2 data is essentially absent):print(qa_c3_t2)## # A tibble: 1 × 6
## n min_months_since_discharge max_months_since_dis…¹ n_with_any_t2_survey…²
## <int> <dbl> <dbl> <int>
## 1 42 6.3 11.2 1
## # ℹ abbreviated names: ¹max_months_since_discharge,
## # ²n_with_any_t2_survey_response
## # ℹ 2 more variables: n_with_t2_a1c <int>, n_with_t2_ed_ip <int># ---- Analysis subsets ----
df_c1 <- df_status %>% filter(cohort_final == "Cohort 1")
df_c2 <- df_status %>% filter(cohort_final == "Cohort 2")
df_c3 <- df_status %>% filter(cohort_final == "Cohort 3")
df_c4 <- df_status %>% filter(cohort_final == "Cohort 4")
df_c5 <- df_status %>% filter(cohort_final == "Cohort 5")
df_not <- df_status %>% filter(cohort_final == "Not enrolled")
df_enrolled <- bind_rows(df_c1, df_c2, df_c3, df_c4, df_c5)
df_analysis <- bind_rows(df_enrolled, df_not)
# Cohorts 4 and 5 are still in their active intervention window - neither has
# any legitimate 9-month/12-month post-intervention data yet. Excluded
# explicitly (rather than relying on those cells happening to be blank) from
# every post-intervention checkpoint below.
df_post <- df_enrolled %>% filter(!cohort_final %in% c("Cohort 4", "Cohort 5"))
# ---- QA: disenrollment signal counts by cohort and type ----
cat("\nDisenrollment signal candidates by cohort (pending case-by-case review):\n")##
## Disenrollment signal candidates by cohort (pending case-by-case review):df_enrolled %>% count(cohort_final, disenroll_flag) %>% filter(disenroll_flag) %>% select(-disenroll_flag) %>% print()## # A tibble: 3 × 2
## cohort_final n
## <chr> <int>
## 1 Cohort 1 10
## 2 Cohort 2 8
## 3 Cohort 3 4cat("\nBy signal type:\n")##
## By signal type:df_enrolled %>% filter(disenroll_flag) %>% count(disenroll_signal, sort = TRUE) %>% print()## # A tibble: 7 × 2
## disenroll_signal n
## <chr> <int>
## 1 Never/stopped receiving deliveries 10
## 2 Address issue 3
## 3 Deceased 2
## 4 Delivery cancelled / hold 2
## 5 Hospitalized / facility 2
## 6 Unresponsive / disengaged 2
## 7 Requested removal 1# ---- QA FLAG: T0 confidence questions were NOT administered on the same
# scale as T1/T2. T1/T2 use a genuine 4-point Likert scale ("Not at all
# confident" -> "Very confident"). At T0, most respondents instead answered
# a plain Yes/No/Maybe question - a much lower bar for a positive answer -
# with only a minority giving a 0-4 numeric or free-text rating. Because of
# this, a T0 -> T2 CHANGE comparison for food-choice/food-label confidence is
# confounded by instrument, not just by time: an apparent decline (or gain)
# may simply reflect the change in question format rather than a true change
# in patient confidence. The composite confidence rows in Table 2
# (conf_composite_t0_bucket / conf_composite_t2_bucket, and
# row_conf_summary) inherit this limitation and should be read with that
# caveat - see manuscript limitations.
qa_conf_t0_format <- bind_rows(
df_enrolled %>% transmute(item = "Food choices", raw = str_squish(as.character(.data[[col_conf_food_t0]]))),
df_enrolled %>% transmute(item = "Food labels", raw = str_squish(as.character(.data[[col_conf_labels_t0]])))
) %>%
filter(is_nonblank(raw)) %>%
mutate(format_type = case_when(
str_to_upper(raw) %in% c("YES", "Y", "NO", "N", "MAYBE", "M") ~ "Yes/No/Maybe",
str_detect(raw, "^[0-4]$") ~ "Numeric 0-4",
TRUE ~ "Free text"
)) %>%
count(item, format_type) %>%
group_by(item) %>%
mutate(pct = round(100 * n / sum(n), 1)) %>%
ungroup()
cat("\nQA FLAG - T0 confidence question response formats (not the T1/T2 4-point scale):\n")##
## QA FLAG - T0 confidence question response formats (not the T1/T2 4-point scale):print(qa_conf_t0_format)## # A tibble: 6 × 4
## item format_type n pct
## <chr> <chr> <int> <dbl>
## 1 Food choices Free text 14 7.8
## 2 Food choices Numeric 0-4 42 23.5
## 3 Food choices Yes/No/Maybe 123 68.7
## 4 Food labels Free text 14 7.8
## 5 Food labels Numeric 0-4 43 24
## 6 Food labels Yes/No/Maybe 122 68.2cat("Implication: a Yes/No/Maybe baseline question sets a much lower bar for a positive\n",
"answer than T1/T2's 4-point scale, so baseline->T2 comparisons for these two items\n",
"(food choices, food labels) are confounded by INSTRUMENT, not just by time. Treat the\n",
"composite confidence rows as descriptive baseline prevalence, not a validated paired\n",
"change metric, until baseline is re-collected on the same 4-point scale.\n", sep = "")## Implication: a Yes/No/Maybe baseline question sets a much lower bar for a positive
## answer than T1/T2's 4-point scale, so baseline->T2 comparisons for these two items
## (food choices, food labels) are confounded by INSTRUMENT, not just by time. Treat the
## composite confidence rows as descriptive baseline prevalence, not a validated paired
## change metric, until baseline is re-collected on the same 4-point scale.============================================================
SECTION 3B: ENROLLMENT FUNNEL & PROCESS MEASURES
Reproduces the eligible -> screened -> positive -> enrolled pipeline
reported in the manuscript's "Enrollment and Population" paragraph and
the PDSA deck's Slide 1 funnel, plus T1/T2 survey completion by cohort.
Previously computed ad hoc in scratch scripts outside this file.
============================================================
funnel <- df %>%
mutate(
zip_5 = str_extract(.data[[col_zip]], "^\\d{5}"),
a1c_numeric = clean_a1c(.data[[col_a1c_all]]),
in_westchester = zip_5 %in% westchester_zips,
meets_a1c = !is.na(a1c_numeric) & a1c_numeric >= 8.0,
eligible = in_westchester & meets_a1c,
t0_contact_clean = case_when(
str_to_upper(str_squish(as.character(.data[[col_t0_contact]]))) == "CONTACTED" ~ "Contacted",
str_to_upper(str_squish(as.character(.data[[col_t0_contact]]))) == "DECLINED" ~ "Declined",
str_to_upper(str_squish(as.character(.data[[col_t0_contact]]))) == "DID NOT MEET CRITERIA" ~ "Did Not Meet Criteria",
str_to_upper(str_squish(as.character(.data[[col_t0_contact]]))) == "OUTSIDE OF DELIVERY RANGE" ~ "Outside of Delivery Range",
str_to_upper(str_squish(as.character(.data[[col_t0_contact]]))) == "LEFT MESSAGE" ~ "Left Message",
str_to_upper(str_squish(as.character(.data[[col_t0_contact]]))) %in% c("UNABLE TO REACH (NO MESSAGE)", "UNABLE TO REACH") ~ "Unable to Reach",
str_to_upper(str_squish(as.character(.data[[col_t0_contact]]))) %in% c("REQUEST FOR CALL BACK", "REQUEST FOR A CALL BACK") ~ "Requested Call Back",
TRUE ~ NA_character_
),
ps_any_yes = case_when(
str_to_upper(str_squish(as.character(.data[[col_ps_worry_t0]]))) %in% c("Y","YES") ~ TRUE,
str_to_upper(str_squish(as.character(.data[[col_ps_lastmonth_t0]]))) %in% c("Y","YES") ~ TRUE,
str_to_upper(str_squish(as.character(.data[[col_ps_skip_t0]]))) %in% c("Y","YES") ~ TRUE,
TRUE ~ FALSE
),
ps_any_filled = is_nonblank(.data[[col_ps_worry_t0]]) |
is_nonblank(.data[[col_ps_lastmonth_t0]]) |
is_nonblank(.data[[col_ps_skip_t0]])
)
eligible_pool <- funnel %>% filter(eligible)
screened <- eligible_pool %>% filter(ps_any_filled)
not_screened <- eligible_pool %>% filter(!ps_any_filled)
pos_pool <- screened %>% filter(ps_any_yes)
enrolled_n <- sum(df_status$cohort_final != "Not enrolled")
funnel_summary <- tibble(
Stage = c("Total in system", "Eligible (Westchester zip + A1C>=8.0)",
"Screened for food insecurity", "Not screened",
"Screened positive", "Screened negative", "Enrolled"),
n = c(nrow(df), nrow(eligible_pool), nrow(screened), nrow(not_screened),
nrow(pos_pool), nrow(screened) - nrow(pos_pool), enrolled_n)
)
cat("Enrollment funnel (reproduces PDSA deck Slide 1 / manuscript Enrollment paragraph):\n")## Enrollment funnel (reproduces PDSA deck Slide 1 / manuscript Enrollment paragraph):print(funnel_summary)## # A tibble: 7 × 2
## Stage n
## <chr> <int>
## 1 Total in system 1406
## 2 Eligible (Westchester zip + A1C>=8.0) 988
## 3 Screened for food insecurity 472
## 4 Not screened 516
## 5 Screened positive 209
## 6 Screened negative 263
## 7 Enrolled 187cat("\nEnrolled as % of screened-positive pool:",
round(100 * enrolled_n / nrow(pos_pool), 1), "%\n")##
## Enrolled as % of screened-positive pool: 89.5 %not_screened_reasons <- not_screened %>% count(t0_contact_clean, sort = TRUE)
cat("\nReasons for not being screened:\n")##
## Reasons for not being screened:print(not_screened_reasons)## # A tibble: 8 × 2
## t0_contact_clean n
## <chr> <int>
## 1 Unable to Reach 179
## 2 Left Message 121
## 3 Declined 77
## 4 Did Not Meet Criteria 63
## 5 Outside of Delivery Range 44
## 6 <NA> 16
## 7 Requested Call Back 12
## 8 Contacted 4unable_pct <- round(100 * sum(not_screened$t0_contact_clean == "Unable to Reach", na.rm = TRUE) / nrow(not_screened), 1)
cat("\n% of not-screened patients coded 'Unable to Reach':", unable_pct, "%\n")##
## % of not-screened patients coded 'Unable to Reach': 34.7 %# ---- T1/T2 survey completion by cohort (contact-flag based) ----
completion <- df_status %>%
filter(cohort_final != "Not enrolled") %>%
mutate(
t1_contact_clean = str_to_upper(str_squish(as.character(t1_contact_y_n))),
t2_contact_clean = str_to_upper(str_squish(as.character(t2_contact_y_n)))
) %>%
group_by(cohort_final) %>%
summarise(
n = n(),
t1_completed = sum(t1_contact_clean == "Y", na.rm = TRUE),
t1_pct = round(100 * t1_completed / n, 1),
t2_completed = sum(t2_contact_clean == "Y", na.rm = TRUE),
t2_pct = round(100 * t2_completed / n, 1),
.groups = "drop"
)
cat("\nT1/T2 survey completion by cohort:\n")##
## T1/T2 survey completion by cohort:print(completion)## # A tibble: 5 × 6
## cohort_final n t1_completed t1_pct t2_completed t2_pct
## <chr> <int> <int> <dbl> <int> <dbl>
## 1 Cohort 1 60 46 76.7 48 80
## 2 Cohort 2 40 28 70 26 65
## 3 Cohort 3 42 19 45.2 0 0
## 4 Cohort 4 36 1 2.8 0 0
## 5 Cohort 5 9 0 0 0 0# ---- Item-level nonresponse: "did nutrition materials help?" ----
# row_help_t1/row_help_t2 in Table 2 use sum(!is.na(v)) as their denominator,
# i.e. "answered this specific question" (Yes/No/etc., anything non-blank).
# That's different from "completed the T1/T2 survey wave" - a patient can be
# actively engaged in that wave (answering several OTHER T1/T2 questions)
# and still have this one specific item blank. This chunk separates the two:
# for each cohort, how many patients were "active" in the wave (answered at
# least one of the other core T1/T2 questions) vs. how many of those active
# patients specifically answered (or left blank) the nutrition-material item.
t1_other_cols <- c(col_ease_t1, col_waste_t1, col_meals_t1, col_conf_food_t1, col_conf_labels_t1, col_conf_prepare_t1, col_help_control)
t2_other_cols <- c(col_ease_t2, col_waste_t2, col_meals_t2, col_conf_food_t2, col_conf_labels_t2, col_conf_prepare_t2, col_help_control_t2)
# Local group list (mirrors groups2 in Section 5, which isn't defined yet at
# this point in the document - built here from the same analysis subsets).
groups_nr <- list(
`Overall Enrolled` = df_enrolled,
`Cohort 1` = df_c1,
`Cohort 2` = df_c2,
`Cohort 3` = df_c3,
`Cohort 4` = df_c4,
`Cohort 5` = df_c5
)
nonresponse_check <- function(d, other_cols, nutrition_col) {
active <- Reduce(`|`, lapply(other_cols, function(cn) is_nonblank(d[[cn]])))
nutrition_answered <- is_nonblank(d[[nutrition_col]])
tibble(
n_active_in_wave = sum(active),
n_answered_nutrition = sum(active & nutrition_answered),
n_missing_nutrition = sum(active & !nutrition_answered),
pct_missing_of_active = if (sum(active) > 0) round(100 * sum(active & !nutrition_answered) / sum(active), 1) else NA_real_
)
}
nonresponse_t1 <- bind_rows(lapply(names(groups_nr), function(gn) {
cbind(Group = gn, nonresponse_check(groups_nr[[gn]], t1_other_cols, col_help_choices))
}))
nonresponse_t2 <- bind_rows(lapply(names(groups_nr), function(gn) {
cbind(Group = gn, nonresponse_check(groups_nr[[gn]], t2_other_cols, col_help_choices_t2))
}))
cat("Nutrition-material item ('did it help?'), T1 - active-in-wave vs. answered vs. missing:\n")## Nutrition-material item ('did it help?'), T1 - active-in-wave vs. answered vs. missing:print(nonresponse_t1)## Group n_active_in_wave n_answered_nutrition n_missing_nutrition
## 1 Overall Enrolled 85 75 10
## 2 Cohort 1 40 35 5
## 3 Cohort 2 26 26 0
## 4 Cohort 3 19 14 5
## 5 Cohort 4 0 0 0
## 6 Cohort 5 0 0 0
## pct_missing_of_active
## 1 11.8
## 2 12.5
## 3 0.0
## 4 26.3
## 5 NA
## 6 NAcat("\nNutrition-material item ('did it help?'), T2 - active-in-wave vs. answered vs. missing:\n")##
## Nutrition-material item ('did it help?'), T2 - active-in-wave vs. answered vs. missing:print(nonresponse_t2)## Group n_active_in_wave n_answered_nutrition n_missing_nutrition
## 1 Overall Enrolled 66 61 5
## 2 Cohort 1 43 43 0
## 3 Cohort 2 23 18 5
## 4 Cohort 3 0 0 0
## 5 Cohort 4 0 0 0
## 6 Cohort 5 0 0 0
## pct_missing_of_active
## 1 7.6
## 2 0.0
## 3 21.7
## 4 NA
## 5 NA
## 6 NAcat("\n'Active in wave' = answered at least one other T1/T2 survey question (ease, waste, meals,\n",
"confidence x3, program-control) that timepoint. 'Missing' = active but this specific item was left\n",
"blank - i.e. genuine item-level nonresponse, not absence from the survey wave altogether.\n", sep = "")##
## 'Active in wave' = answered at least one other T1/T2 survey question (ease, waste, meals,
## confidence x3, program-control) that timepoint. 'Missing' = active but this specific item was left
## blank - i.e. genuine item-level nonresponse, not absence from the survey wave altogether.============================================================
SECTION 4: TABLE 1 - BASELINE CHARACTERISTICS
============================================================
groups1 <- list(
Total = df_analysis,
`Cohort 1` = df_c1,
`Cohort 2` = df_c2,
`Cohort 3` = df_c3,
`Cohort 4` = df_c4,
`Cohort 5` = df_c5,
`Not enrolled` = df_not
)
make_row1 <- function(label, cells) {
row <- as.data.frame(matrix(cells, nrow = 1), stringsAsFactors = FALSE)
names(row) <- names(cells)
cbind(Characteristic = label, row, stringsAsFactors = FALSE)
}
# Embed a cohort's p-value (vs. not-enrolled) directly into its own cell,
# under the value it belongs to - bold + asterisked when significant.
# NOTE: built as a plain HTML string, not via kableExtra::cell_spec() -
# kable_styling()/row_spec()/pack_rows() re-escape already-inserted raw HTML
# in cells (confirmed independently of this document), so Table 1 is
# rendered with a small custom HTML builder below instead of the usual
# kable() %>% kable_styling() chain.
value_p_cell <- function(value_str, p) {
sig <- !is.na(p) & p < 0.05
p_str <- paste0("p = ", fmt_p(p), if (isTRUE(sig)) "*" else "")
weight <- if (isTRUE(sig)) "font-weight:bold;" else ""
paste0("<span style='", weight, "'>", value_str, "<br><span style='font-size:85%;color:#555;", weight, "'>", p_str, "</span></span>")
}
# For a categorical header row: the test applies to the whole variable, not
# one value, so the cohort's cell holds just the p-value
p_only_cell <- function(p) {
sig <- !is.na(p) & p < 0.05
txt <- paste0("p = ", fmt_p(p), if (isTRUE(sig)) "*" else "")
weight <- if (isTRUE(sig)) "font-weight:bold;" else ""
paste0("<span style='", weight, "'>", txt, "</span>")
}
# ---- Significance tests reused by Table 2 ----
# Paired Wilcoxon signed-rank test on an already-computed change/difference
# vector: is the median change significantly different from zero?
p_change_from_zero <- function(x) {
x <- x[!is.na(x)]
if (length(x) < 3 || length(unique(x)) < 2) return(NA_real_)
suppressWarnings(wilcox.test(x, mu = 0)$p.value)
}
# Exact sign test on Improved vs. Worsened counts ("No change" pairs
# excluded) - the standard test for "did more patients improve than worsen"
p_sign_test <- function(n_improved, n_worsened) {
n_total <- n_improved + n_worsened
if (n_total < 1) return(NA_real_)
suppressWarnings(binom.test(n_improved, n_total, p = 0.5)$p.value)
}
# Continuous: one-way ANOVA + Dunnett's test (each cohort vs. not-enrolled,
# multiplicity-adjusted in one step - see dunnett_pvals in Section 2)
cont_row1 <- function(label, col, digits = 2) {
raw_p <- dunnett_pvals(df_analysis[[col]], df_analysis$cohort_final)
cohort_cells <- setNames(
vapply(COHORT_LEVELS, function(cl) value_p_cell(fmt_mean_sd(groups1[[cl]][[col]], digits), raw_p[[cl]]), character(1)),
COHORT_LEVELS
)
cells <- c(Total = fmt_mean_sd(df_analysis[[col]], digits), cohort_cells,
`Not enrolled` = fmt_mean_sd(df_not[[col]], digits))
make_row1(label, cells)
}
# Categorical: chi-square/Fisher per cohort vs. not-enrolled, Bonferroni-
# adjusted for the simultaneous comparisons (one per cohort)
cat_block1 <- function(label, col, levels_vec) {
raw_p <- cat_pvals_bonf(col, groups1, df_not)
cohort_p_cells <- setNames(
vapply(COHORT_LEVELS, function(cl) p_only_cell(raw_p[[cl]]), character(1)),
COHORT_LEVELS
)
header <- make_row1(label, c(Total = "", cohort_p_cells, `Not enrolled` = ""))
rows <- do.call(rbind, lapply(levels_vec, function(lv) {
cells <- vapply(groups1, function(d) cat_stat(d[[col]], lv), character(1))
make_row1(lv, cells)
}))
rbind(header, rows)
}
demo_block <- rbind(
cont_row1("Age, years", "age_num"),
cat_block1("Sex", "sex_group", c("Female", "Male")),
cat_block1("Race", "race_group", c("White", "Black", "Other")),
cat_block1("Ethnicity", "ethnicity_display", c(
"Not Spanish / Hispanic / Latino", "Spanish / Hispanic / Latino", "Declined / Unavailable")),
cat_block1("Insurance", "insurance_group", c(
"Commercial", "Medicare", "Medicaid", "Dual Eligible", "Uninsured / Self Pay"))
)
clinical_block <- rbind(
cont_row1("HbA1c Value (Last 3 Months)", "a1c_base"),
cont_row1("LACE+ Readmission Score", "lace_num")
)
util_block <- rbind(
cont_row1("Hospitalizations (pre, 90 days)", "ip_base"),
cont_row1("ED visits (pre, 90 days)", "ed_base")
)
# ---- Attrition / disenrollment (enrolled patients only - no not-enrolled
# comparison exists for this, so plain counts, no p-value) ----
disenroll_cells <- function(match_fn) {
cohort_cells <- setNames(
vapply(COHORT_LEVELS, function(cl) fmt_np(sum(match_fn(groups1[[cl]])), nrow(groups1[[cl]])), character(1)),
COHORT_LEVELS
)
c(Total = fmt_np(sum(match_fn(df_enrolled)), nrow(df_enrolled)), cohort_cells, `Not enrolled` = "–")
}
row_disenroll <- make_row1("Disenrollment signal identified (T1/T2 chart-review notes)",
disenroll_cells(function(d) d$disenroll_flag))
disenroll_types <- df_enrolled %>% filter(disenroll_flag) %>%
count(disenroll_signal, sort = TRUE) %>% pull(disenroll_signal)
disenroll_type_rows <- do.call(rbind, lapply(disenroll_types, function(sig) {
make_row1(paste0(" ", sig),
disenroll_cells(function(d) !is.na(d$disenroll_signal) & d$disenroll_signal == sig))
}))
attrition_block <- rbind(row_disenroll, disenroll_type_rows)
table1_df <- rbind(demo_block, clinical_block, util_block, attrition_block)
header_rows <- which(table1_df$Characteristic %in% c("Sex", "Race", "Ethnicity", "Insurance",
"Disenrollment signal identified (T1/T2 chart-review notes)"))
n_demo <- nrow(demo_block); n_clin <- nrow(clinical_block); n_util <- nrow(util_block)
render_html_table(table1_df, "Table 1. Baseline Characteristics",
bold_rows = header_rows,
sections = list(
"Demographics & Socioeconomic Characteristics" = 1,
"Clinical Characteristics" = n_demo + 1,
"Baseline Healthcare Utilization (90 days pre-index)" = n_demo + n_clin + 1,
"Attrition / Disenrollment (Enrolled Only)" = n_demo + n_clin + n_util + 1
),
footnote = paste0(
"Age, HbA1c value (last 3 months), LACE+ readmission score, hospitalizations in the past 90 days, and ED ",
"visits in the past 90 days are continuous variables summarized as mean (SD), compared with a one-way ANOVA ",
"across Cohorts 1-5 and Not-enrolled followed by Dunnett's test (each cohort vs. the not-enrolled reference; ",
"p-values shown are already multiplicity-adjusted for the 5 simultaneous comparisons). Sex, race, ethnicity, ",
"and insurance are categorical variables summarized as n (%); each cohort was compared to not-enrolled with a ",
"chi-square test (Fisher's exact with simulated p-values when expected cell counts were small), Bonferroni-",
"adjusted for the 5 comparisons and shown on the variable's header row. p-values are each cohort vs. ",
"not-enrolled; bold with * = p < .05. Race was collapsed into White, Black, and Other for analysis. ",
"Ethnicity's Declined and Unavailable responses were combined (n=48 and n=3) given their small size. Insurance ",
"types were grouped into Commercial, Medicare, Medicaid, Dual Eligible, and Uninsured / Self Pay. ",
"Disenrollment is a text-pattern screen over T1/T2 chart-review notes (deceased, requested removal, ",
"hospitalized/facility, address issues, unresponsive, stopped receiving deliveries, etc.) - these are ",
"candidates flagged for case-by-case clinical review, not a confirmed disenrollment count; \"did not receive ",
"educational materials\" alone is explicitly excluded as a signal. Not applicable to the not-enrolled group. ",
"Cohort 5 (n=", nrow(df_c5), ") is newly enrolled with only baseline data so far - included throughout, but ",
"expect \"-\" for any T1/T2/post-intervention cell."
))Table 1. Baseline Characteristics
| Characteristic | Total | Cohort 1 | Cohort 2 | Cohort 3 | Cohort 4 | Cohort 5 | Not enrolled |
|---|---|---|---|---|---|---|---|
| Demographics & Socioeconomic Characteristics | |||||||
| Age, years | 63.86 (15.81) | 59.10 (14.02) p = 0.036* | 55.27 (14.52) p = 0.001* | 58.45 (16.02) p = 0.057 | 60.86 (14.25) p = 0.558 | 60.22 (17.03) p = 0.920 | 64.67 (15.82) |
| Sex | p = 0.510 | p = 1.000 | p = 1.000 | p = 1.000 | p = 1.000 | ||
| Female | 625 (44.5%) | 33 (55.0%) | 20 (50.0%) | 21 (50.0%) | 20 (55.6%) | 2 (22.2%) | 529 (43.4%) |
| Male | 781 (55.5%) | 27 (45.0%) | 20 (50.0%) | 21 (50.0%) | 16 (44.4%) | 7 (77.8%) | 690 (56.6%) |
| Race | p = 0.005* | p = 0.121 | p = 0.009* | p = 0.290 | p = 0.300 | ||
| White | 576 (41.0%) | 14 (23.3%) | 9 (22.5%) | 8 (19.0%) | 9 (25.0%) | 1 (11.1%) | 535 (43.9%) |
| Black | 356 (25.3%) | 25 (41.7%) | 14 (35.0%) | 18 (42.9%) | 13 (36.1%) | 2 (22.2%) | 284 (23.3%) |
| Other | 474 (33.7%) | 21 (35.0%) | 17 (42.5%) | 16 (38.1%) | 14 (38.9%) | 6 (66.7%) | 400 (32.8%) |
| Ethnicity | p = 0.401 | p = 0.490 | p = 1.000 | p = 1.000 | p = 0.057 | ||
| Not Spanish / Hispanic / Latino | 921 (65.6%) | 36 (60.0%) | 21 (52.5%) | 27 (64.3%) | 21 (58.3%) | 2 (22.2%) | 814 (66.9%) |
| Spanish / Hispanic / Latino | 432 (30.8%) | 24 (40.0%) | 18 (45.0%) | 15 (35.7%) | 14 (38.9%) | 7 (77.8%) | 354 (29.1%) |
| Declined / Unavailable | 51 (3.6%) | 0 (0.0%) | 1 (2.5%) | 0 (0.0%) | 1 (2.8%) | 0 (0.0%) | 49 (4.0%) |
| Insurance | p = 1.000 | p = 0.153 | p = 0.041* | p = 1.000 | p = 0.146 | ||
| Commercial | 366 (26.0%) | 18 (30.0%) | 10 (25.0%) | 11 (26.2%) | 9 (25.0%) | 0 (0.0%) | 318 (26.1%) |
| Medicare | 657 (46.7%) | 23 (38.3%) | 12 (30.0%) | 12 (28.6%) | 14 (38.9%) | 3 (33.3%) | 593 (48.6%) |
| Medicaid | 321 (22.8%) | 17 (28.3%) | 17 (42.5%) | 15 (35.7%) | 12 (33.3%) | 6 (66.7%) | 254 (20.8%) |
| Dual Eligible | 50 (3.6%) | 2 (3.3%) | 1 (2.5%) | 2 (4.8%) | 1 (2.8%) | 0 (0.0%) | 44 (3.6%) |
| Uninsured / Self Pay | 12 (0.9%) | 0 (0.0%) | 0 (0.0%) | 2 (4.8%) | 0 (0.0%) | 0 (0.0%) | 10 (0.8%) |
| Clinical Characteristics | |||||||
| HbA1c Value (Last 3 Months) | 9.95 (2.15) | 10.20 (2.10) p = 0.821 | 10.23 (2.08) p = 0.879 | 9.95 (2.00) p = 1.000 | 10.50 (2.68) p = 0.417 | 11.70 (3.27) p = 0.061 | 9.90 (2.13) |
| LACE+ Readmission Score | 58.38 (16.06) | 59.33 (14.58) p = 0.998 | 54.42 (17.68) p = 0.452 | 56.45 (15.95) p = 0.931 | 59.36 (14.67) p = 0.999 | 56.11 (19.58) p = 0.995 | 58.51 (16.10) |
| Baseline Healthcare Utilization (90 days pre-index) | |||||||
| Hospitalizations (pre, 90 days) | 1.52 (1.26) | 0.67 (0.80) p = <0.001* | 1.60 (0.98) p = 1.000 | 1.36 (0.88) p = 0.796 | 1.25 (0.77) p = 0.495 | 1.33 (0.71) p = 0.985 | 1.57 (1.30) |
| ED visits (pre, 90 days) | 1.69 (1.75) | 1.13 (2.17) p = 0.041* | 1.60 (1.10) p = 0.991 | 1.40 (1.06) p = 0.707 | 1.31 (0.95) p = 0.526 | 1.89 (1.05) p = 1.000 | 1.74 (1.78) |
| Attrition / Disenrollment (Enrolled Only) | |||||||
| Disenrollment signal identified (T1/T2 chart-review notes) | 22 (11.8%) | 10 (16.7%) | 8 (20.0%) | 4 (9.5%) | 0 (0.0%) | 0 (0.0%) | – |
| Never/stopped receiving deliveries | 10 (5.3%) | 4 (6.7%) | 4 (10.0%) | 2 (4.8%) | 0 (0.0%) | 0 (0.0%) | – |
| Address issue | 3 (1.6%) | 2 (3.3%) | 1 (2.5%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) | – |
| Deceased | 2 (1.1%) | 0 (0.0%) | 1 (2.5%) | 1 (2.4%) | 0 (0.0%) | 0 (0.0%) | – |
| Delivery cancelled / hold | 2 (1.1%) | 2 (3.3%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) | – |
| Hospitalized / facility | 2 (1.1%) | 2 (3.3%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) | – |
| Unresponsive / disengaged | 2 (1.1%) | 0 (0.0%) | 2 (5.0%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) | – |
| Requested removal | 1 (0.5%) | 0 (0.0%) | 0 (0.0%) | 1 (2.4%) | 0 (0.0%) | 0 (0.0%) | – |
Notes: Age, HbA1c value (last 3 months), LACE+ readmission score, hospitalizations in the past 90 days, and ED visits in the past 90 days are continuous variables summarized as mean (SD), compared with a one-way ANOVA across Cohorts 1-5 and Not-enrolled followed by Dunnett's test (each cohort vs. the not-enrolled reference; p-values shown are already multiplicity-adjusted for the 5 simultaneous comparisons). Sex, race, ethnicity, and insurance are categorical variables summarized as n (%); each cohort was compared to not-enrolled with a chi-square test (Fisher's exact with simulated p-values when expected cell counts were small), Bonferroni-adjusted for the 5 comparisons and shown on the variable's header row. p-values are each cohort vs. not-enrolled; bold with * = p < .05. Race was collapsed into White, Black, and Other for analysis. Ethnicity's Declined and Unavailable responses were combined (n=48 and n=3) given their small size. Insurance types were grouped into Commercial, Medicare, Medicaid, Dual Eligible, and Uninsured / Self Pay. Disenrollment is a text-pattern screen over T1/T2 chart-review notes (deceased, requested removal, hospitalized/facility, address issues, unresponsive, stopped receiving deliveries, etc.) - these are candidates flagged for case-by-case clinical review, not a confirmed disenrollment count; "did not receive educational materials" alone is explicitly excluded as a signal. Not applicable to the not-enrolled group. Cohort 5 (n=9) is newly enrolled with only baseline data so far - included throughout, but expect "-" for any T1/T2/post-intervention cell.
============================================================
SECTION 5: TABLE 2 - PRIMARY AND SECONDARY OUTCOMES
============================================================
groups2 <- list(
`Overall Enrolled` = df_enrolled,
`Cohort 1` = df_c1,
`Cohort 2` = df_c2,
`Cohort 3` = df_c3,
`Cohort 4` = df_c4,
`Cohort 5` = df_c5
)
make_row2 <- function(label, cells) {
row <- as.data.frame(matrix(cells, nrow = 1), stringsAsFactors = FALSE)
names(row) <- names(cells)
cbind(Outcome = label, row, stringsAsFactors = FALSE)
}
# Header row (n paired) + Improved/No change/Worsened sub-rows.
# `label` should include its own timeframe (e.g. ", paired T1->T2") since
# different measures here span different windows. Default classify_fn is a
# plain numeric comparison (higher = better); pass one of the bucket_*
# functions below for measures where a one-step move should count as
# "No change" rather than a real improvement/worsening.
default_change_classify <- function(b, a) case_when(a > b ~ "Improved", a == b ~ "No change", TRUE ~ "Worsened")
bucket_change_classify <- function(b, a) case_when(
b == "Low" & a == "High" ~ "Improved",
b == "High" & a == "Low" ~ "Worsened",
TRUE ~ "No change" # includes staying in the same bucket
)
paired_change_block2 <- function(label, before_col, after_col, classify_fn = default_change_classify) {
classify <- function(d) {
b <- d[[before_col]]; a <- d[[after_col]]
out <- rep(NA_character_, length(b))
ok <- !is.na(b) & !is.na(a)
out[ok] <- classify_fn(b[ok], a[ok])
out
}
n_paired <- vapply(groups2, function(d) sum(!is.na(d[[before_col]]) & !is.na(d[[after_col]])), integer(1))
# Significance: exact sign test on Improved vs. Worsened (ties excluded) -
# is improvement more common than worsening among those who changed at all?
header_cells <- vapply(names(groups2), function(gn) {
cls <- classify(groups2[[gn]])
p <- p_sign_test(sum(cls == "Improved", na.rm = TRUE), sum(cls == "Worsened", na.rm = TRUE))
value_p_cell(paste0("n=", n_paired[[gn]]), p)
}, character(1))
header <- make_row2(label, setNames(header_cells, names(groups2)))
rows <- do.call(rbind, lapply(c("Improved", "No change", "Worsened"), function(lv) {
cells <- vapply(groups2, function(d) {
cls <- classify(d)
fmt_nN(sum(cls == lv, na.rm = TRUE), sum(!is.na(cls)))
}, character(1))
make_row2(paste0(" ", lv), cells)
}))
rbind(header, rows)
}
# One-line version of the block above, for Table 2's headline summary - full
# Improved/No change/Worsened detail (and the before/after prevalence) moves
# to the Patient Experience Detail section instead of bulking up Table 2
paired_change_summary_row2 <- function(label, before_col, after_col, classify_fn = default_change_classify) {
classify <- function(d) {
b <- d[[before_col]]; a <- d[[after_col]]
out <- rep(NA_character_, length(b))
ok <- !is.na(b) & !is.na(a)
out[ok] <- classify_fn(b[ok], a[ok])
out
}
cells <- vapply(names(groups2), function(gn) {
cls <- classify(groups2[[gn]])
n <- sum(!is.na(cls))
if (n == 0) return("–")
n_imp <- sum(cls == "Improved", na.rm = TRUE); n_wor <- sum(cls == "Worsened", na.rm = TRUE)
p <- p_sign_test(n_imp, n_wor)
val <- paste0(round(100 * n_imp / n, 1), "% improved, ", round(100 * n_wor / n, 1), "% worsened (n=", n, ")")
value_p_cell(val, p)
}, character(1))
make_row2(label, setNames(cells, names(groups2)))
}
# Continuous change row with a per-column paired Wilcoxon signed-rank test
# (H0: median change = 0), embedded inline like Table 1's p-values
change_row_with_p <- function(label, value_fn) {
cells <- vapply(names(groups2), function(gn) {
x <- value_fn(groups2[[gn]])
value_p_cell(fmt_mean_sd_n(x), p_change_from_zero(x))
}, character(1))
make_row2(label, setNames(cells, names(groups2)))
}
# ---- Primary outcome ----
# T2 is the endpoint/graduation, so this change score already represents the
# full intervention effect (baseline -> T2). See the trajectory tables in
# Section 6 for the midpoint (T1) breakdown.
row_a1c_change <- change_row_with_p("HbA1c change, baseline → T2 (full intervention), mean (SD)",
function(d) d$a1c_change_t2_calc)
row_a1c_reduction <- make_row2("HbA1c reduction ≥1.0% at T2 (full intervention)",
vapply(groups2, function(d) {
x <- d$a1c_change_t2_calc[!is.na(d$a1c_change_t2_calc)]
fmt_nN(sum(x <= -1.0), length(x))
}, character(1)))
# ---- Secondary: healthcare utilization ----
row_ed_base <- make_row2("ED visits at baseline (past 90 days), mean (SD)",
vapply(groups2, function(d) fmt_mean_sd_n(d$ed_base), character(1)))
row_ip_base <- make_row2("Hospitalizations at baseline (past 90 days), mean (SD)",
vapply(groups2, function(d) fmt_mean_sd_n(d$ip_base), character(1)))
row_ed_change <- change_row_with_p("Change in ED visits, baseline → T2 (full intervention), mean (SD)",
function(d) { ok <- !is.na(d$ed_base) & !is.na(d$ed_t2); d$ed_t2[ok] - d$ed_base[ok] })
row_ip_change <- change_row_with_p("Change in hospitalizations, baseline → T2 (full intervention), mean (SD)",
function(d) { ok <- !is.na(d$ip_base) & !is.na(d$ip_t2); d$ip_t2[ok] - d$ip_base[ok] })
# ---- Secondary: food insecurity ----
row_fi_t0 <- make_row2("Food insecurity positive at baseline (T0)",
vapply(groups2, function(d) fmt_nN(sum(d$fi_t0_pos, na.rm = TRUE), sum(d$fi_t0_any)), character(1)))
row_fi_t2 <- make_row2("Food insecurity positive at T2",
vapply(groups2, function(d) fmt_nN(sum(d$fi_t2_pos, na.rm = TRUE), sum(d$fi_t2_any)), character(1)))
# Paired significance: exact sign test (McNemar-style) on Resolved
# (positive->negative) vs. Newly positive (negative->positive), among
# patients with a definitive Yes/No at both T0 and T2 ("No change" pairs
# excluded from the test, same convention as ease/confidence below)
fi_classify <- function(b, a) case_when(b == TRUE & a == FALSE ~ "Improved",
b == FALSE & a == TRUE ~ "Worsened", TRUE ~ "No change")
row_fi_summary <- paired_change_summary_row2("Food insecurity resolved, paired baseline (T0)→T2",
"fi_t0_pos", "fi_t2_pos", classify_fn = fi_classify)
# ---- Secondary: patient experience ----
# Each question's T1 and T2 versions are asked and answered independently
# (different, non-overlapping respondents in general) - shown as separate
# rows with their own n so they're never mistaken for a pooled or paired
# figure, unlike the waste/meals rows below which explicitly pool T1+T2.
row_help_t1 <- make_row2("Nutrition material helped healthier food choices (Yes) - T1",
vapply(groups2, function(d) {
v <- clean_impact(d[[col_help_choices]])
fmt_nN(sum(v %in% c("A lot", "Somewhat", "Yes"), na.rm = TRUE), sum(!is.na(v)))
}, character(1)))
row_help_t2 <- make_row2("Nutrition material helped healthier food choices (Yes) - T2",
vapply(groups2, function(d) {
v <- clean_impact(d[[col_help_choices_t2]])
fmt_nN(sum(v %in% c("A lot", "Somewhat", "Yes"), na.rm = TRUE), sum(!is.na(v)))
}, character(1)))
# ---- Item-level nonresponse for the two rows above: how many patients were
# "active" that wave (answered at least one other T1/T2 question - ease,
# waste, meals, confidence x3, program-control) but left THIS question
# blank, vs. how many actually answered it. Distinguishes a genuine skipped
# item from simple absence from the survey wave (see Section 3B for the
# full per-cohort breakdown this summarizes). ----
t1_other_cols <- c(col_ease_t1, col_waste_t1, col_meals_t1, col_conf_food_t1, col_conf_labels_t1, col_conf_prepare_t1, col_help_control)
t2_other_cols <- c(col_ease_t2, col_waste_t2, col_meals_t2, col_conf_food_t2, col_conf_labels_t2, col_conf_prepare_t2, col_help_control_t2)
item_nonresponse_row <- function(label, other_cols, nutrition_col) {
make_row2(label, vapply(groups2, function(d) {
active <- Reduce(`|`, lapply(other_cols, function(cn) is_nonblank(d[[cn]])))
answered <- is_nonblank(d[[nutrition_col]])
fmt_nN(sum(active & !answered), sum(active))
}, character(1)))
}
row_help_t1_missing <- item_nonresponse_row(
" Nutrition-material item left blank despite active T1 survey wave", t1_other_cols, col_help_choices)
row_help_t2_missing <- item_nonresponse_row(
" Nutrition-material item left blank despite active T2 survey wave", t2_other_cols, col_help_choices_t2)
row_control_t1 <- make_row2("Program helped patient feel more in control (A lot / Somewhat) - T1",
vapply(groups2, function(d) {
v <- clean_impact(d[[col_help_control]])
fmt_nN(sum(v %in% c("A lot", "Somewhat"), na.rm = TRUE), sum(!is.na(v)))
}, character(1)))
row_control_t2 <- make_row2("Program helped patient feel more in control (A lot / Somewhat) - T2",
vapply(groups2, function(d) {
v <- clean_impact(d[[col_help_control_t2]])
fmt_nN(sum(v %in% c("A lot", "Somewhat"), na.rm = TRUE), sum(!is.na(v)))
}, character(1)))
row_waste <- make_row2("Reported no food waste (T1 or T2, pooled)",
vapply(groups2, function(d) {
resp <- c(d$waste_t1, d$waste_t2)
resp <- resp[!is.na(resp)]
fmt_nN(sum(!resp), length(resp))
}, character(1)))
row_meals <- make_row2("Meals prepared from program food per week, mean (SD) (T1 or T2, pooled)",
vapply(groups2, function(d) fmt_mean_sd_n(c(d$meals_t1, d$meals_t2)), character(1)))
# ---- Prevalence rows: the actual % reporting each answer at each
# timepoint, both the positive response AND its opposite, so the paired
# Improved/No change/Worsened blocks below aren't the only view - a reader
# can see directly whether "easy" (or "confident") became more or less
# common over the intervention, not just how individuals moved.
ease_prevalence_row <- function(timepoint_label, col) {
make_row2(paste0("Food ease at ", timepoint_label, ": Easy/Very easy vs. everything else"),
vapply(groups2, function(d) {
v <- clean_ease(d[[col]]); v <- v[!is.na(v)]; n <- length(v)
if (n == 0) return("–")
pos <- sum(v %in% c("Easy", "Very easy")); neg <- n - pos
paste0(round(100 * pos / n, 1), "% easy/very easy, ", round(100 * neg / n, 1),
"% neutral/difficult (n=", n, ")")
}, character(1)))
}
row_ease_t1_prev <- ease_prevalence_row("T1", col_ease_t1)
row_ease_t2_prev <- ease_prevalence_row("T2", col_ease_t2)
# Ease has no baseline (T0) equivalent - patients can't rate ease of using
# delivered food before any food has been delivered - so T1->T2 is the only
# valid window, unlike confidence below. Any one-step move (Neutral<->Easy<->
# Very easy) is treated as a real, ordered change since ease has no natural
# "confident bucket"-style grouping the way the 4-point confidence scale does.
ease_block <- paired_change_block2("Food ease, paired T1→T2", "ease_t1_score", "ease_t2_score")
conf_prevalence_row <- function(label, bucket_col) {
make_row2(label, vapply(groups2, function(d) {
b <- d[[bucket_col]]; b <- b[!is.na(b)]; n <- length(b)
if (n == 0) return("–")
pos <- sum(b == "High"); neg <- sum(b == "Low")
paste0(round(100 * pos / n, 1), "% high, ", round(100 * neg / n, 1), "% low (n=", n, ")")
}, character(1)))
}
row_conf_composite_t0_prev <- conf_prevalence_row("Confidence (food choices + labels) at baseline (T0): High vs. Low", "conf_composite_t0_bucket")
row_conf_composite_t2_prev <- conf_prevalence_row("Confidence (food choices + labels) at T2: High vs. Low", "conf_composite_t2_bucket")
# Food choices + labels have a genuine baseline (T0) item, so compared
# baseline -> T2 (bucket-based, per the lateral-movement fix above).
conf_block <- paired_change_block2("Patient confidence (food choices + labels), paired baseline→T2",
"conf_composite_t0_bucket", "conf_composite_t2_bucket", classify_fn = bucket_change_classify)
row_conf_prepare_t1_prev <- conf_prevalence_row("Confidence in meal preparation at T1: High vs. Low", "conf_prepare_t1_bucket")
row_conf_prepare_t2_prev <- conf_prevalence_row("Confidence in meal preparation at T2: High vs. Low", "conf_prepare_t2_bucket")
# Meal prep confidence was never asked at T0, so it's kept as its own T1->T2
# comparison rather than folded into the baseline-anchored composite above.
conf_prepare_block <- paired_change_block2("Confidence in meal preparation, paired T1→T2 (no T0 baseline collected)",
"conf_prepare_t1_bucket", "conf_prepare_t2_bucket", classify_fn = bucket_change_classify)
row_ease_summary <- paired_change_summary_row2("Food ease improved, paired T1→T2",
"ease_t1_score", "ease_t2_score")
row_conf_summary <- paired_change_summary_row2("Patient confidence (food choices + labels) improved, paired baseline→T2",
"conf_composite_t0_bucket", "conf_composite_t2_bucket", classify_fn = bucket_change_classify)
row_conf_prepare_summary <- paired_change_summary_row2("Confidence in meal preparation improved, paired T1→T2",
"conf_prepare_t1_bucket", "conf_prepare_t2_bucket", classify_fn = bucket_change_classify)
primary_block <- rbind(row_a1c_change, row_a1c_reduction)
util_block2 <- rbind(row_ed_base, row_ip_base, row_ed_change, row_ip_change)
fi_block2 <- rbind(row_fi_t0, row_fi_t2, row_fi_summary)
experience_block <- rbind(row_help_t1, row_help_t1_missing, row_help_t2, row_help_t2_missing,
row_control_t1, row_control_t2, row_waste, row_meals,
row_ease_summary, row_conf_summary, row_conf_prepare_summary)
table2_df <- rbind(primary_block, util_block2, fi_block2, experience_block)
render_html_table(table2_df, "Table 2. Primary and Secondary Outcomes",
sections = list(
"Primary Outcome" = 1,
"Secondary Outcomes: Healthcare Utilization" = nrow(primary_block) + 1,
"Secondary Outcomes: Food Insecurity" = nrow(primary_block) + nrow(util_block2) + 1,
"Secondary Outcomes: Patient Experience" = nrow(primary_block) + nrow(util_block2) + nrow(fi_block2) + 1
),
footnote = paste0(
"Outcomes are summarized among enrolled participants (Cohorts 1-5). T2 is the endpoint/graduation, so ",
"\"baseline → T2\" change rows represent the full-intervention effect (see Sections 7-9 for the baseline → T1 ",
"midpoint breakdown and the 9/12-month post-intervention trajectory; see Section 6 below for the full before/",
"after prevalence and paired Improved/No change/Worsened detail behind the food ease and confidence \"improved\" ",
"summary rows here). Continuous change rows (HbA1c, ED, IP) are tested with a paired Wilcoxon signed-rank test ",
"(H0: median change = 0); the \"improved\"/\"resolved\" summary rows (food insecurity, food ease, confidence) ",
"are tested with an exact sign test on Improved vs. Worsened counts (\"No change\" pairs excluded) - is ",
"improvement more common than worsening among those who changed at all. Food insecurity's test compares ",
"Resolved (positive at T0, negative at T2) vs. Newly positive (negative at T0, positive at T2), among patients ",
"with a definitive answer at both timepoints. Bold with * = p < .05 in both cases. Rows reporting a single ",
"proportion (reduction ≥1%, food insecurity prevalence, nutrition/control impact, food waste) have no ",
"comparable reference group and are not significance-tested. Nutrition-material and program-control questions ",
"are shown separately by T1 and T2 (each with its own respondent count) rather than pooled, since - unlike the ",
"food waste/meals rows below, which explicitly pool T1+T2 - these are largely non-overlapping respondents and ",
"a pooled percentage would not be directly comparable to either timepoint alone. The indented \"left blank ",
"despite active survey wave\" rows show item-level nonresponse specifically: among patients who answered at ",
"least one OTHER T1/T2 question that wave (ease, waste, meals, confidence x3, program-control), how many left ",
"the nutrition-material question itself blank - i.e. this one item was skipped, not that the patient was ",
"absent from the wave altogether (see Section 3B for the full per-cohort breakdown). For ",
"HbA1c, negative change reflects improvement; computed as the checkpoint's absolute value minus baseline value ",
"(the spreadsheet's own a1c_change_t2 column was found to disagree with the underlying lab values for most of ",
"Cohort 2 and to contain artifacts for Cohort 3 - see the QA check in the Data Preparation section; Cohort 3 ",
"also has essentially no T2 data at all despite being well past its graduation window - see the Cohort 3 QA ",
"flag in that same section). The primary outcome is defined as a reduction in HbA1c of 1.0% or more from ",
"baseline to T2. Food insecurity was defined as a positive response to one or more of the three Hunger Vital ",
"Sign items at that timepoint. Baseline utilization reflects the past 90 days; change was calculated as T2 ",
"minus baseline. Food waste pools the T1 and T2 delivery-waste responses together. Meals/week pools T1 and T2 ",
"responses (free text - digits, spelled-out numbers, and ranges are parsed to a number; ambiguous entries like ",
"\"3 a day\" or \"unsure\" are excluded). Food ease is scored Difficult/Somewhat difficult/Neutral/Easy/Very ",
"easy (no baseline item exists - patients can't rate ease of using delivered food before ",
"enrollment). Confidence is bucketed as High (Confident/Very confident) vs. Low (Not/Somewhat not confident); ",
"a one-step move within a bucket (e.g. Very confident to Confident) counts as no change, not a real ",
"improvement or worsening. Food choices and food labels have a genuine baseline (T0) item and are compared as ",
"a composite (High/Low only when both items agree); meal preparation confidence was never asked at T0, so it ",
"is T1-to-T2 instead."
))Table 2. Primary and Secondary Outcomes
| Outcome | Overall Enrolled | Cohort 1 | Cohort 2 | Cohort 3 | Cohort 4 | Cohort 5 |
|---|---|---|---|---|---|---|
| Primary Outcome | ||||||
| HbA1c change, baseline → T2 (full intervention), mean (SD) | -1.78 (2.43), n=49 p = <0.001* | -2.05 (2.36), n=30 p = <0.001* | -1.35 (2.55), n=19 p = 0.052 | – p = – | – p = – | – p = – |
| HbA1c reduction ≥1.0% at T2 (full intervention) | 27/49 (55.1%) | 18/30 (60.0%) | 9/19 (47.4%) | – | – | – |
| Secondary Outcomes: Healthcare Utilization | ||||||
| ED visits at baseline (past 90 days), mean (SD) | 1.36 (1.50), n=187 | 1.13 (2.17), n=60 | 1.60 (1.10), n=40 | 1.40 (1.06), n=42 | 1.31 (0.95), n=36 | 1.89 (1.05), n=9 |
| Hospitalizations at baseline (past 90 days), mean (SD) | 1.17 (0.92), n=187 | 0.67 (0.80), n=60 | 1.60 (0.98), n=40 | 1.36 (0.88), n=42 | 1.25 (0.77), n=36 | 1.33 (0.71), n=9 |
| Change in ED visits, baseline → T2 (full intervention), mean (SD) | -0.80 (1.45), n=81 p = <0.001* | -0.64 (1.46), n=58 p = <0.001* | -1.22 (1.38), n=23 p = <0.001* | – p = – | – p = – | – p = – |
| Change in hospitalizations, baseline → T2 (full intervention), mean (SD) | -0.84 (1.04), n=82 p = <0.001* | -0.56 (0.77), n=59 p = <0.001* | -1.57 (1.27), n=23 p = <0.001* | – p = – | – p = – | – p = – |
| Secondary Outcomes: Food Insecurity | ||||||
| Food insecurity positive at baseline (T0) | 184/185 (99.5%) | 59/60 (98.3%) | 40/40 (100.0%) | 42/42 (100.0%) | 34/34 (100.0%) | 9/9 (100.0%) |
| Food insecurity positive at T2 | 43/67 (64.2%) | 28/43 (65.1%) | 15/24 (62.5%) | – | – | – |
| Food insecurity resolved, paired baseline (T0)→T2 | 34.3% improved, 0% worsened (n=67) p = <0.001* | 32.6% improved, 0% worsened (n=43) p = <0.001* | 37.5% improved, 0% worsened (n=24) p = 0.004* | – | – | – |
| Secondary Outcomes: Patient Experience | ||||||
| Nutrition material helped healthier food choices (Yes) - T1 | 67/75 (89.3%) | 32/35 (91.4%) | 25/26 (96.2%) | 10/14 (71.4%) | – | – |
| Nutrition-material item left blank despite active T1 survey wave | 10/85 (11.8%) | 5/40 (12.5%) | 0/26 (0.0%) | 5/19 (26.3%) | – | – |
| Nutrition material helped healthier food choices (Yes) - T2 | 48/61 (78.7%) | 36/43 (83.7%) | 12/18 (66.7%) | – | – | – |
| Nutrition-material item left blank despite active T2 survey wave | 5/66 (7.6%) | 0/43 (0.0%) | 5/23 (21.7%) | – | – | – |
| Program helped patient feel more in control (A lot / Somewhat) - T1 | 65/83 (78.3%) | 33/38 (86.8%) | 19/26 (73.1%) | 13/19 (68.4%) | – | – |
| Program helped patient feel more in control (A lot / Somewhat) - T2 | 56/63 (88.9%) | 41/43 (95.3%) | 15/20 (75.0%) | – | – | – |
| Reported no food waste (T1 or T2, pooled) | 121/143 (84.6%) | 67/77 (87.0%) | 39/47 (83.0%) | 15/19 (78.9%) | – | – |
| Meals prepared from program food per week, mean (SD) (T1 or T2, pooled) | 3.89 (2.06), n=145 | 3.33 (1.80), n=79 | 4.14 (2.08), n=48 | 5.69 (2.00), n=18 | – | – |
| Food ease improved, paired T1→T2 | 21.4% improved, 21.4% worsened (n=42) p = 1.000 | 11.5% improved, 30.8% worsened (n=26) p = 0.227 | 37.5% improved, 6.2% worsened (n=16) p = 0.125 | – | – | – |
| Patient confidence (food choices + labels) improved, paired baseline→T2 | 10% improved, 0% worsened (n=30) p = 0.250 | 0% improved, 0% worsened (n=15) p = – | 20% improved, 0% worsened (n=15) p = 0.250 | – | – | – |
| Confidence in meal preparation improved, paired T1→T2 | 6.2% improved, 18.8% worsened (n=32) p = 0.289 | 9.5% improved, 23.8% worsened (n=21) p = 0.453 | 0% improved, 9.1% worsened (n=11) p = 1.000 | – | – | – |
Notes: Outcomes are summarized among enrolled participants (Cohorts 1-5). T2 is the endpoint/graduation, so "baseline → T2" change rows represent the full-intervention effect (see Sections 7-9 for the baseline → T1 midpoint breakdown and the 9/12-month post-intervention trajectory; see Section 6 below for the full before/after prevalence and paired Improved/No change/Worsened detail behind the food ease and confidence "improved" summary rows here). Continuous change rows (HbA1c, ED, IP) are tested with a paired Wilcoxon signed-rank test (H0: median change = 0); the "improved"/"resolved" summary rows (food insecurity, food ease, confidence) are tested with an exact sign test on Improved vs. Worsened counts ("No change" pairs excluded) - is improvement more common than worsening among those who changed at all. Food insecurity's test compares Resolved (positive at T0, negative at T2) vs. Newly positive (negative at T0, positive at T2), among patients with a definitive answer at both timepoints. Bold with * = p < .05 in both cases. Rows reporting a single proportion (reduction ≥1%, food insecurity prevalence, nutrition/control impact, food waste) have no comparable reference group and are not significance-tested. Nutrition-material and program-control questions are shown separately by T1 and T2 (each with its own respondent count) rather than pooled, since - unlike the food waste/meals rows below, which explicitly pool T1+T2 - these are largely non-overlapping respondents and a pooled percentage would not be directly comparable to either timepoint alone. The indented "left blank despite active survey wave" rows show item-level nonresponse specifically: among patients who answered at least one OTHER T1/T2 question that wave (ease, waste, meals, confidence x3, program-control), how many left the nutrition-material question itself blank - i.e. this one item was skipped, not that the patient was absent from the wave altogether (see Section 3B for the full per-cohort breakdown). For HbA1c, negative change reflects improvement; computed as the checkpoint's absolute value minus baseline value (the spreadsheet's own a1c_change_t2 column was found to disagree with the underlying lab values for most of Cohort 2 and to contain artifacts for Cohort 3 - see the QA check in the Data Preparation section; Cohort 3 also has essentially no T2 data at all despite being well past its graduation window - see the Cohort 3 QA flag in that same section). The primary outcome is defined as a reduction in HbA1c of 1.0% or more from baseline to T2. Food insecurity was defined as a positive response to one or more of the three Hunger Vital Sign items at that timepoint. Baseline utilization reflects the past 90 days; change was calculated as T2 minus baseline. Food waste pools the T1 and T2 delivery-waste responses together. Meals/week pools T1 and T2 responses (free text - digits, spelled-out numbers, and ranges are parsed to a number; ambiguous entries like "3 a day" or "unsure" are excluded). Food ease is scored Difficult/Somewhat difficult/Neutral/Easy/Very easy (no baseline item exists - patients can't rate ease of using delivered food before enrollment). Confidence is bucketed as High (Confident/Very confident) vs. Low (Not/Somewhat not confident); a one-step move within a bucket (e.g. Very confident to Confident) counts as no change, not a real improvement or worsening. Food choices and food labels have a genuine baseline (T0) item and are compared as a composite (High/Low only when both items agree); meal preparation confidence was never asked at T0, so it is T1-to-T2 instead.
============================================================
SECTION 6: PATIENT EXPERIENCE DETAIL (SUPPLEMENTARY)
============================================================
# Full detail behind Table 2's condensed "improved"/"resolved" rows: the
# before/after prevalence (both the positive response and its opposite) and
# the complete paired Improved/No change/Worsened breakdown with
# significance, for food insecurity plus the 3 patient-experience measures.
detail_fi <- paired_change_block2("Food insecurity resolved, paired baseline (T0)→T2",
"fi_t0_pos", "fi_t2_pos", classify_fn = fi_classify)
detail_ease <- rbind(row_ease_t1_prev, row_ease_t2_prev, ease_block)
detail_conf <- rbind(row_conf_composite_t0_prev, row_conf_composite_t2_prev, conf_block)
detail_prepare <- rbind(row_conf_prepare_t1_prev, row_conf_prepare_t2_prev, conf_prepare_block)
# ---- Confidence (food choices, food labels), T1->T2 same-scale alternative ----
# The composite above is anchored to baseline (T0), which used a mostly
# yes/no/maybe format rather than T1/T2's 4-point scale (see qa_conf_t0_format
# in Section 3B) - a real instrument-mismatch limitation, not just noise. T1
# and T2 use the identical scale for every item, so this block reports each
# item's T1->T2 change on its own (no T0 involved) as a methodologically
# clean alternative - the same framing used for the corrected PDSA-deck
# confidence chart.
row_conf_food_t1_prev <- conf_prevalence_row("Food choices confidence at T1: High vs. Low", "conf_food_t1_bucket")
row_conf_food_t2_prev <- conf_prevalence_row("Food choices confidence at T2: High vs. Low", "conf_food_t2_bucket")
conf_food_t1t2_block <- paired_change_block2("Food choices confidence improved, paired T1→T2",
"conf_food_t1_bucket", "conf_food_t2_bucket", classify_fn = bucket_change_classify)
row_conf_labels_t1_prev <- conf_prevalence_row("Food labels confidence at T1: High vs. Low", "conf_labels_t1_bucket")
row_conf_labels_t2_prev <- conf_prevalence_row("Food labels confidence at T2: High vs. Low", "conf_labels_t2_bucket")
conf_labels_t1t2_block <- paired_change_block2("Food labels confidence improved, paired T1→T2",
"conf_labels_t1_bucket", "conf_labels_t2_bucket", classify_fn = bucket_change_classify)
detail_conf_t1t2 <- rbind(row_conf_food_t1_prev, row_conf_food_t2_prev, conf_food_t1t2_block,
row_conf_labels_t1_prev, row_conf_labels_t2_prev, conf_labels_t1t2_block)
detail_df <- rbind(detail_fi, detail_ease, detail_conf, detail_conf_t1t2, detail_prepare)
detail_header_rows <- which(detail_df$Outcome %in% c(
"Food insecurity resolved, paired baseline (T0)→T2",
"Food ease, paired T1→T2",
"Patient confidence (food choices + labels), paired baseline→T2",
"Food choices confidence improved, paired T1→T2",
"Food labels confidence improved, paired T1→T2",
"Confidence in meal preparation, paired T1→T2 (no T0 baseline collected)"
))
render_html_table(detail_df, "Patient Experience Detail: Food Insecurity, Food Ease, and Confidence",
bold_rows = detail_header_rows,
sections = list(
"Food Insecurity" = 1,
"Food Ease" = nrow(detail_fi) + 1,
"Patient Confidence (Food Choices + Labels), Baseline→T2" = nrow(detail_fi) + nrow(detail_ease) + 1,
"Patient Confidence (Food Choices, Food Labels), T1→T2 Same-Scale Alternative" =
nrow(detail_fi) + nrow(detail_ease) + nrow(detail_conf) + 1,
"Confidence in Meal Preparation" =
nrow(detail_fi) + nrow(detail_ease) + nrow(detail_conf) + nrow(detail_conf_t1t2) + 1
),
footnote = paste0(
"Supplementary detail for the \"improved\"/\"resolved\" summary rows in Table 2. Food insecurity's Improved = ",
"resolved (positive at T0, negative at T2); Worsened = newly positive (negative at T0, positive at T2), among ",
"patients with a definitive Yes/No at both timepoints. \"At [timepoint]: ... vs. ...\" rows for food ease and ",
"confidence are simple prevalence snapshots (% reporting the positive response and % reporting the opposite, ",
"independent of pairing), shown before and after so the shift in how common each answer is over the ",
"intervention is visible directly. The paired Improved/No change/Worsened block below each one instead tracks ",
"the same individuals (exact sign test on Improved vs. Worsened, ties excluded; bold with * = p < .05). ",
"Food ease has no baseline (T0) item (patients can't rate ease of using delivered food before enrollment), so ",
"it is T1 to T2 only. Confidence is bucketed High (Confident/Very confident) vs. Low (Not/Somewhat not ",
"confident) - a one-step move within a bucket (e.g. Very confident to Confident) counts as no change. Food ",
"choices and food labels have a genuine baseline (T0) item and are compared as a composite (High/Low only ",
"when both items agree; a mixed pair is excluded as ambiguous); T0's response format was inconsistent with ",
"T1/T2 (plain 0-4 values, Yes/No, free text) - 0-1 was treated as Low, 2 as Unsure (excluded from the paired ",
"comparison), 3-4 as High. The \"T1->T2 Same-Scale Alternative\" block reports food choices and food labels ",
"confidence again, this time comparing T1 to T2 only (both on the identical 4-point scale) instead of baseline ",
"(T0) to T2 - the same framing used for the corrected PDSA-deck confidence chart, free of the T0 instrument ",
"mismatch noted above. Meal preparation confidence was never asked at T0, so it is shown as its own ",
"T1-to-T2 comparison."
))Patient Experience Detail: Food Insecurity, Food Ease, and Confidence
| Outcome | Overall Enrolled | Cohort 1 | Cohort 2 | Cohort 3 | Cohort 4 | Cohort 5 |
|---|---|---|---|---|---|---|
| Food Insecurity | ||||||
| Food insecurity resolved, paired baseline (T0)→T2 | n=67 p = <0.001* | n=43 p = <0.001* | n=24 p = 0.004* | n=0 p = – | n=0 p = – | n=0 p = – |
| Improved | 23/67 (34.3%) | 14/43 (32.6%) | 9/24 (37.5%) | – | – | – |
| No change | 44/67 (65.7%) | 29/43 (67.4%) | 15/24 (62.5%) | – | – | – |
| Worsened | 0/67 (0.0%) | 0/43 (0.0%) | 0/24 (0.0%) | – | – | – |
| Food Ease | ||||||
| Food ease at T1: Easy/Very easy vs. everything else | 94.8% easy/very easy, 5.2% neutral/difficult (n=77) | 100% easy/very easy, 0% neutral/difficult (n=35) | 91.7% easy/very easy, 8.3% neutral/difficult (n=24) | 88.9% easy/very easy, 11.1% neutral/difficult (n=18) | – | – |
| Food ease at T2: Easy/Very easy vs. everything else | 96.8% easy/very easy, 3.2% neutral/difficult (n=62) | 95% easy/very easy, 5% neutral/difficult (n=40) | 100% easy/very easy, 0% neutral/difficult (n=22) | – | – | – |
| Food ease, paired T1→T2 | n=42 p = 1.000 | n=26 p = 0.227 | n=16 p = 0.125 | n=0 p = – | n=0 p = – | n=0 p = – |
| Improved | 9/42 (21.4%) | 3/26 (11.5%) | 6/16 (37.5%) | – | – | – |
| No change | 24/42 (57.1%) | 15/26 (57.7%) | 9/16 (56.2%) | – | – | – |
| Worsened | 9/42 (21.4%) | 8/26 (30.8%) | 1/16 (6.2%) | – | – | – |
| Patient Confidence (Food Choices + Labels), Baseline→T2 | ||||||
| Confidence (food choices + labels) at baseline (T0): High vs. Low | 91.3% high, 8.7% low (n=127) | 94.3% high, 5.7% low (n=35) | 90.6% high, 9.4% low (n=32) | 89.3% high, 10.7% low (n=28) | 89.3% high, 10.7% low (n=28) | 100% high, 0% low (n=4) |
| Confidence (food choices + labels) at T2: High vs. Low | 92.9% high, 7.1% low (n=42) | 88% high, 12% low (n=25) | 100% high, 0% low (n=17) | – | – | – |
| Patient confidence (food choices + labels), paired baseline→T2 | n=30 p = 0.250 | n=15 p = – | n=15 p = 0.250 | n=0 p = – | n=0 p = – | n=0 p = – |
| Improved | 3/30 (10.0%) | 0/15 (0.0%) | 3/15 (20.0%) | – | – | – |
| No change | 27/30 (90.0%) | 15/15 (100.0%) | 12/15 (80.0%) | – | – | – |
| Worsened | 0/30 (0.0%) | 0/15 (0.0%) | 0/15 (0.0%) | – | – | – |
| Patient Confidence (Food Choices, Food Labels), T1→T2 Same-Scale Alternative | ||||||
| Food choices confidence at T1: High vs. Low | 90.9% high, 9.1% low (n=77) | 86.5% high, 13.5% low (n=37) | 100% high, 0% low (n=21) | 89.5% high, 10.5% low (n=19) | – | – |
| Food choices confidence at T2: High vs. Low | 88.7% high, 11.3% low (n=53) | 84.8% high, 15.2% low (n=33) | 95% high, 5% low (n=20) | – | – | – |
| Food choices confidence improved, paired T1→T2 | n=36 p = 1.000 | n=24 p = 1.000 | n=12 p = – | n=0 p = – | n=0 p = – | n=0 p = – |
| Improved | 3/36 (8.3%) | 3/24 (12.5%) | 0/12 (0.0%) | – | – | – |
| No change | 29/36 (80.6%) | 17/24 (70.8%) | 12/12 (100.0%) | – | – | – |
| Worsened | 4/36 (11.1%) | 4/24 (16.7%) | 0/12 (0.0%) | – | – | – |
| Food labels confidence at T1: High vs. Low | 91.1% high, 8.9% low (n=79) | 83.8% high, 16.2% low (n=37) | 95.7% high, 4.3% low (n=23) | 100% high, 0% low (n=19) | – | – |
| Food labels confidence at T2: High vs. Low | 80.4% high, 19.6% low (n=56) | 75% high, 25% low (n=36) | 90% high, 10% low (n=20) | – | – | – |
| Food labels confidence improved, paired T1→T2 | n=36 p = 0.508 | n=24 p = 0.508 | n=12 p = – | n=0 p = – | n=0 p = – | n=0 p = – |
| Improved | 3/36 (8.3%) | 3/24 (12.5%) | 0/12 (0.0%) | – | – | – |
| No change | 27/36 (75.0%) | 15/24 (62.5%) | 12/12 (100.0%) | – | – | – |
| Worsened | 6/36 (16.7%) | 6/24 (25.0%) | 0/12 (0.0%) | – | – | – |
| Confidence in Meal Preparation | ||||||
| Confidence in meal preparation at T1: High vs. Low | 94.6% high, 5.4% low (n=74) | 90.9% high, 9.1% low (n=33) | 100% high, 0% low (n=22) | 94.7% high, 5.3% low (n=19) | – | – |
| Confidence in meal preparation at T2: High vs. Low | 88.7% high, 11.3% low (n=53) | 85.3% high, 14.7% low (n=34) | 94.7% high, 5.3% low (n=19) | – | – | – |
| Confidence in meal preparation, paired T1→T2 (no T0 baseline collected) | n=32 p = 0.289 | n=21 p = 0.453 | n=11 p = 1.000 | n=0 p = – | n=0 p = – | n=0 p = – |
| Improved | 2/32 (6.2%) | 2/21 (9.5%) | 0/11 (0.0%) | – | – | – |
| No change | 24/32 (75.0%) | 14/21 (66.7%) | 10/11 (90.9%) | – | – | – |
| Worsened | 6/32 (18.8%) | 5/21 (23.8%) | 1/11 (9.1%) | – | – | – |
Notes: Supplementary detail for the "improved"/"resolved" summary rows in Table 2. Food insecurity's Improved = resolved (positive at T0, negative at T2); Worsened = newly positive (negative at T0, positive at T2), among patients with a definitive Yes/No at both timepoints. "At [timepoint]: ... vs. ..." rows for food ease and confidence are simple prevalence snapshots (% reporting the positive response and % reporting the opposite, independent of pairing), shown before and after so the shift in how common each answer is over the intervention is visible directly. The paired Improved/No change/Worsened block below each one instead tracks the same individuals (exact sign test on Improved vs. Worsened, ties excluded; bold with * = p < .05). Food ease has no baseline (T0) item (patients can't rate ease of using delivered food before enrollment), so it is T1 to T2 only. Confidence is bucketed High (Confident/Very confident) vs. Low (Not/Somewhat not confident) - a one-step move within a bucket (e.g. Very confident to Confident) counts as no change. Food choices and food labels have a genuine baseline (T0) item and are compared as a composite (High/Low only when both items agree; a mixed pair is excluded as ambiguous); T0's response format was inconsistent with T1/T2 (plain 0-4 values, Yes/No, free text) - 0-1 was treated as Low, 2 as Unsure (excluded from the paired comparison), 3-4 as High. The "T1->T2 Same-Scale Alternative" block reports food choices and food labels confidence again, this time comparing T1 to T2 only (both on the identical 4-point scale) instead of baseline (T0) to T2 - the same framing used for the corrected PDSA-deck confidence chart, free of the T0 instrument mismatch noted above. Meal preparation confidence was never asked at T0, so it is shown as its own T1-to-T2 comparison.
============================================================
SECTION 7: A1C PROGRESSION - BASELINE -> T1 -> T2 -> 9MO -> 12MO
============================================================
# Snapshot stats at one checkpoint (not necessarily paired to any other)
traj_stat <- function(x, label) {
x <- x[!is.na(x)]
tibble(timepoint = label, n = length(x),
mean = if (length(x) > 0) round(mean(x), 2) else NA_real_,
sd = if (length(x) > 1) round(sd(x), 2) else NA_real_)
}
# Paired change between two checkpoints (same patients only)
paired_stat <- function(before, after, label) {
ok <- !is.na(before) & !is.na(after)
b <- before[ok]; a <- after[ok]; n <- length(b)
if (n < 3) {
return(tibble(Comparison = label, n = n, `Mean change (SD)` = "–",
`Median change` = "–", `p-value` = "–"))
}
d <- a - b
wx <- suppressWarnings(wilcox.test(a, b, paired = TRUE))
tibble(Comparison = label, n = n,
`Mean change (SD)` = fmt_mean_sd(d),
`Median change` = as.character(round(median(d), 2)),
`p-value` = fmt_p(wx$p.value))
}
TP_LEVELS <- c("Baseline", "T1", "T2", "9 Months Post", "12 Months Post")# ---- QA: why ED/IP n at 9mo/12mo is much smaller than A1C n at 9mo/12mo ----
# This is a real pattern in the source data, not a parsing bug: at the
# 9-month and 12-month waves, the chart-review team recorded an A1C value for
# many patients but explicitly entered "N/A" text (not a blank cell) in the
# matching ED/IP utilization fields - i.e., utilization simply was not
# re-pulled at those waves for most patients. Verify directly against the
# spreadsheet: ed_utilization_9months_post, hospital_utilization_9months_post,
# ed_utilization_12_months_post, hospital_utilization_12_months_post.
completeness_qa <- tibble(
Checkpoint = c("9 Months Post", "12 Months Post"),
`A1C values available` = c(sum(!is.na(df_post$a1c_9mo)), sum(!is.na(df_post$a1c_12mo))),
`ED utilization available` = c(sum(!is.na(df_post$ed_9mo)), sum(!is.na(df_post$ed_12mo))),
`IP utilization available` = c(sum(!is.na(df_post$ip_9mo)), sum(!is.na(df_post$ip_12mo)))
)
cat("Checkpoint completeness (Cohorts 1-3 only - Cohort 4 excluded, still active):\n")## Checkpoint completeness (Cohorts 1-3 only - Cohort 4 excluded, still active):print(completeness_qa)## # A tibble: 2 × 4
## Checkpoint `A1C values available` `ED utilization available`
## <chr> <int> <int>
## 1 9 Months Post 28 1
## 2 12 Months Post 2 0
## # ℹ 1 more variable: `IP utilization available` <int>a1c_traj <- bind_rows(
traj_stat(df_enrolled$a1c_base, "Baseline"),
traj_stat(df_enrolled$a1c_t1, "T1"),
traj_stat(df_enrolled$a1c_t2, "T2"),
traj_stat(df_post$a1c_9mo, "9 Months Post"),
traj_stat(df_post$a1c_12mo, "12 Months Post")
) %>%
mutate(timepoint = factor(timepoint, levels = TP_LEVELS),
label = paste0(mean, "%\n(n=", n, ")"))
kable(a1c_traj %>% transmute(Timepoint = timepoint, n,
`Mean HbA1c (%)` = ifelse(is.na(mean), "–", as.character(mean)),
SD = ifelse(is.na(sd), "–", as.character(sd))),
caption = "A1C by Checkpoint (Enrolled Patients, Cohorts 1-5; 9/12-month excludes Cohorts 4-5 - still active)") %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = FALSE, font_size = 13) %>%
row_spec(0, bold = TRUE, background = "#1A2B4A", color = "white")| Timepoint | n | Mean HbA1c (%) | SD |
|---|---|---|---|
| Baseline | 187 | 10.28 | 2.26 |
| T1 | 68 | 8.16 | 1.72 |
| T2 | 49 | 8.58 | 2.16 |
| 9 Months Post | 28 | 9.15 | 3 |
| 12 Months Post | 2 | 9.8 | 0.85 |
ggplot(a1c_traj %>% filter(n > 0), aes(x = timepoint, y = mean, group = 1)) +
geom_line(color = "#0D7A6E", linewidth = 1.6) +
geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.12, linewidth = 0.9, color = "#0D7A6E") +
geom_point(color = "#0D7A6E", size = 4.5) +
geom_text(aes(label = label, y = mean + sd), vjust = -0.4, size = 3.6, fontface = "bold", color = "#1A2F4A") +
scale_y_continuous(expand = expansion(mult = c(0.08, 0.22))) +
labs(title = "Mean HbA1c by Checkpoint (All Enrolled Patients)",
subtitle = "Error bars = ± SD · 9/12-month checkpoints exclude Cohorts 4-5 (still active)",
x = NULL, y = "Mean HbA1c (%)") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold", color = "#1A2F4A", size = 14),
plot.subtitle = element_text(color = "#64748B", size = 9),
panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(),
axis.text.x = element_text(face = "bold", color = "#1A2F4A", size = 11, angle = 20, hjust = 1))
a1c_pairwise <- bind_rows(
paired_stat(df_enrolled$a1c_base, df_enrolled$a1c_t1, "Baseline → T1 (midpoint)"),
paired_stat(df_enrolled$a1c_t1, df_enrolled$a1c_t2, "T1 → T2 (endpoint)"),
paired_stat(df_enrolled$a1c_base, df_enrolled$a1c_t2, "Baseline → T2 (full intervention)"),
paired_stat(df_post$a1c_t2, df_post$a1c_9mo, "T2 → 9 Months Post"),
paired_stat(df_post$a1c_9mo, df_post$a1c_12mo, "9 → 12 Months Post"),
paired_stat(df_post$a1c_base, df_post$a1c_9mo, "Baseline → 9 Months Post"),
paired_stat(df_post$a1c_base, df_post$a1c_12mo, "Baseline → 12 Months Post")
)
kable(a1c_pairwise, caption = "A1C Change Between Checkpoints (paired patients only)",
align = c("l", "c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = TRUE, font_size = 13) %>%
row_spec(0, bold = TRUE, background = "#1A2B4A", color = "white") %>%
row_spec(3, bold = TRUE, background = "#E4F5F3") %>%
footnote(general = paste0(
"Negative change = improvement. Paired Wilcoxon signed-rank test. Baseline is patient-relative ",
"(chart-pulled from each patient's own 3 months pre-enrollment). T1/T2/9-month/12-month are chart-review ",
"checkpoints, not calendar-verified follow-up windows. Rows involving 9- or 12-month checkpoints exclude ",
"Cohorts 4-5 (still in their active intervention window)."),
general_title = "Notes: ", escape = FALSE)| Comparison | n | Mean change (SD) | Median change | p-value |
|---|---|---|---|---|
| Baseline → T1 (midpoint) | 68 | -2.11 (2.19) | -1.75 | <0.001 |
| T1 → T2 (endpoint) | 35 | -0.05 (1.03) | 0 | 1.000 |
| Baseline → T2 (full intervention) | 49 | -1.78 (2.43) | -1.5 | <0.001 |
| T2 → 9 Months Post | 16 | 0.18 (2.21) | 0.3 | 0.816 |
| 9 → 12 Months Post | 1 | – | – | – |
| Baseline → 9 Months Post | 28 | -1.11 (3.17) | -0.95 | 0.077 |
| Baseline → 12 Months Post | 2 | – | – | – |
| Notes: | ||||
| Negative change = improvement. Paired Wilcoxon signed-rank test. Baseline is patient-relative (chart-pulled from each patient's own 3 months pre-enrollment). T1/T2/9-month/12-month are chart-review checkpoints, not calendar-verified follow-up windows. Rows involving 9- or 12-month checkpoints exclude Cohorts 4-5 (still in their active intervention window). |
============================================================
SECTION 8: ED UTILIZATION PROGRESSION - BASELINE -> T1 -> T2 -> 9MO -> 12MO
============================================================
ed_traj <- bind_rows(
traj_stat(df_enrolled$ed_base, "Baseline"),
traj_stat(df_enrolled$ed_t1, "T1"),
traj_stat(df_enrolled$ed_t2, "T2"),
traj_stat(df_post$ed_9mo, "9 Months Post"),
traj_stat(df_post$ed_12mo, "12 Months Post")
) %>%
mutate(timepoint = factor(timepoint, levels = TP_LEVELS),
label = paste0(mean, "\n(n=", n, ")"))
kable(ed_traj %>% transmute(Timepoint = timepoint, n,
`Mean ED visits` = ifelse(is.na(mean), "–", as.character(mean)),
SD = ifelse(is.na(sd), "–", as.character(sd))),
caption = "ED Visits by Checkpoint (Enrolled Patients, Cohorts 1-5; 9/12-month excludes Cohorts 4-5 - still active)") %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = FALSE, font_size = 13) %>%
row_spec(0, bold = TRUE, background = "#1A2B4A", color = "white")| Timepoint | n | Mean ED visits | SD |
|---|---|---|---|
| Baseline | 187 | 1.36 | 1.5 |
| T1 | 142 | 0.4 | 0.83 |
| T2 | 81 | 0.53 | 1.23 |
| 9 Months Post | 1 | 1 | – |
| 12 Months Post | 0 | – | – |
ggplot(ed_traj %>% filter(n > 0), aes(x = timepoint, y = mean, group = 1)) +
geom_line(color = "#1F4E79", linewidth = 1.6) +
geom_errorbar(aes(ymin = pmax(0, mean - sd), ymax = mean + sd), width = 0.12, linewidth = 0.9, color = "#1F4E79") +
geom_point(color = "#1F4E79", size = 4.5) +
geom_text(aes(label = label, y = mean + sd), vjust = -0.4, size = 3.6, fontface = "bold", color = "#1A2F4A") +
scale_y_continuous(expand = expansion(mult = c(0.08, 0.22))) +
labs(title = "Mean ED Visits per Patient by Checkpoint",
subtitle = "Error bars = ± SD · 9/12-month checkpoints exclude Cohorts 4-5 (still active)",
x = NULL, y = "Mean ED visits") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold", color = "#1A2F4A", size = 14),
plot.subtitle = element_text(color = "#64748B", size = 9),
panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(),
axis.text.x = element_text(face = "bold", color = "#1A2F4A", size = 11, angle = 20, hjust = 1))
ed_pairwise <- bind_rows(
paired_stat(df_enrolled$ed_base, df_enrolled$ed_t1, "Baseline → T1 (midpoint)"),
paired_stat(df_enrolled$ed_t1, df_enrolled$ed_t2, "T1 → T2 (endpoint)"),
paired_stat(df_enrolled$ed_base, df_enrolled$ed_t2, "Baseline → T2 (full intervention)"),
paired_stat(df_post$ed_t2, df_post$ed_9mo, "T2 → 9 Months Post"),
paired_stat(df_post$ed_9mo, df_post$ed_12mo, "9 → 12 Months Post"),
paired_stat(df_post$ed_base, df_post$ed_9mo, "Baseline → 9 Months Post"),
paired_stat(df_post$ed_base, df_post$ed_12mo, "Baseline → 12 Months Post")
)
kable(ed_pairwise, caption = "ED Visit Change Between Checkpoints (paired patients only)",
align = c("l", "c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = TRUE, font_size = 13) %>%
row_spec(0, bold = TRUE, background = "#1A2B4A", color = "white") %>%
row_spec(3, bold = TRUE, background = "#E4F5F3") %>%
footnote(general = paste0(
"Negative change = fewer ED visits. Paired Wilcoxon signed-rank test. Rows involving 9- or 12-month ",
"checkpoints exclude Cohorts 4-5 (still in their active intervention window)."),
general_title = "Notes: ", escape = FALSE)| Comparison | n | Mean change (SD) | Median change | p-value |
|---|---|---|---|---|
| Baseline → T1 (midpoint) | 142 | -0.94 (1.41) | -1 | <0.001 |
| T1 → T2 (endpoint) | 81 | 0.06 (0.98) | 0 | 0.678 |
| Baseline → T2 (full intervention) | 81 | -0.80 (1.45) | -1 | <0.001 |
| T2 → 9 Months Post | 1 | – | – | – |
| 9 → 12 Months Post | 0 | – | – | – |
| Baseline → 9 Months Post | 1 | – | – | – |
| Baseline → 12 Months Post | 0 | – | – | – |
| Notes: | ||||
| Negative change = fewer ED visits. Paired Wilcoxon signed-rank test. Rows involving 9- or 12-month checkpoints exclude Cohorts 4-5 (still in their active intervention window). |
============================================================
SECTION 9: IP UTILIZATION PROGRESSION - BASELINE -> T1 -> T2 -> 9MO -> 12MO
============================================================
ip_traj <- bind_rows(
traj_stat(df_enrolled$ip_base, "Baseline"),
traj_stat(df_enrolled$ip_t1, "T1"),
traj_stat(df_enrolled$ip_t2, "T2"),
traj_stat(df_post$ip_9mo, "9 Months Post"),
traj_stat(df_post$ip_12mo, "12 Months Post")
) %>%
mutate(timepoint = factor(timepoint, levels = TP_LEVELS),
label = paste0(mean, "\n(n=", n, ")"))
kable(ip_traj %>% transmute(Timepoint = timepoint, n,
`Mean IP stays` = ifelse(is.na(mean), "–", as.character(mean)),
SD = ifelse(is.na(sd), "–", as.character(sd))),
caption = "IP Stays by Checkpoint (Enrolled Patients, Cohorts 1-5; 9/12-month excludes Cohorts 4-5 - still active)") %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = FALSE, font_size = 13) %>%
row_spec(0, bold = TRUE, background = "#1A2B4A", color = "white")| Timepoint | n | Mean IP stays | SD |
|---|---|---|---|
| Baseline | 187 | 1.17 | 0.92 |
| T1 | 142 | 0.22 | 0.51 |
| T2 | 82 | 0.12 | 0.33 |
| 9 Months Post | 0 | – | – |
| 12 Months Post | 0 | – | – |
ggplot(ip_traj %>% filter(n > 0), aes(x = timepoint, y = mean, group = 1)) +
geom_line(color = "#7B4F9E", linewidth = 1.6) +
geom_errorbar(aes(ymin = pmax(0, mean - sd), ymax = mean + sd), width = 0.12, linewidth = 0.9, color = "#7B4F9E") +
geom_point(color = "#7B4F9E", size = 4.5) +
geom_text(aes(label = label, y = mean + sd), vjust = -0.4, size = 3.6, fontface = "bold", color = "#1A2F4A") +
scale_y_continuous(expand = expansion(mult = c(0.08, 0.22))) +
labs(title = "Mean IP Stays per Patient by Checkpoint",
subtitle = "Error bars = ± SD · 9/12-month checkpoints exclude Cohorts 4-5 (still active)",
x = NULL, y = "Mean IP stays") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold", color = "#1A2F4A", size = 14),
plot.subtitle = element_text(color = "#64748B", size = 9),
panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(),
axis.text.x = element_text(face = "bold", color = "#1A2F4A", size = 11, angle = 20, hjust = 1))
ip_pairwise <- bind_rows(
paired_stat(df_enrolled$ip_base, df_enrolled$ip_t1, "Baseline → T1 (midpoint)"),
paired_stat(df_enrolled$ip_t1, df_enrolled$ip_t2, "T1 → T2 (endpoint)"),
paired_stat(df_enrolled$ip_base, df_enrolled$ip_t2, "Baseline → T2 (full intervention)"),
paired_stat(df_post$ip_t2, df_post$ip_9mo, "T2 → 9 Months Post"),
paired_stat(df_post$ip_9mo, df_post$ip_12mo, "9 → 12 Months Post"),
paired_stat(df_post$ip_base, df_post$ip_9mo, "Baseline → 9 Months Post"),
paired_stat(df_post$ip_base, df_post$ip_12mo, "Baseline → 12 Months Post")
)
kable(ip_pairwise, caption = "IP Stay Change Between Checkpoints (paired patients only)",
align = c("l", "c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = TRUE, font_size = 13) %>%
row_spec(0, bold = TRUE, background = "#1A2B4A", color = "white") %>%
row_spec(3, bold = TRUE, background = "#E4F5F3") %>%
footnote(general = paste0(
"Negative change = fewer IP stays. Paired Wilcoxon signed-rank test. Rows involving 9- or 12-month ",
"checkpoints exclude Cohorts 4-5 (still in their active intervention window)."),
general_title = "Notes: ", escape = FALSE)| Comparison | n | Mean change (SD) | Median change | p-value |
|---|---|---|---|---|
| Baseline → T1 (midpoint) | 142 | -0.92 (1.05) | -1 | <0.001 |
| T1 → T2 (endpoint) | 82 | -0.13 (0.60) | 0 | 0.047 |
| Baseline → T2 (full intervention) | 82 | -0.84 (1.04) | -1 | <0.001 |
| T2 → 9 Months Post | 0 | – | – | – |
| 9 → 12 Months Post | 0 | – | – | – |
| Baseline → 9 Months Post | 0 | – | – | – |
| Baseline → 12 Months Post | 0 | – | – | – |
| Notes: | ||||
| Negative change = fewer IP stays. Paired Wilcoxon signed-rank test. Rows involving 9- or 12-month checkpoints exclude Cohorts 4-5 (still in their active intervention window). |