DM FIM Program Outcomes
2026-09-09
============================================================
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_control <- "t1_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_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"
# ---- 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",
str_detect(s, "VERY EASY") ~ "Very easy",
s == "EASY" ~ "Easy",
str_detect(s, "NEUTRAL") ~ "Neutral",
TRUE ~ NA_character_
)
}
# Ordinal score for clean_ease(), for paired before/after comparison
ease_score <- function(x) {
v <- clean_ease(x)
case_when(v == "Neutral" ~ 1, v == "Easy" ~ 2, v == "Very easy" ~ 3, TRUE ~ NA_real_)
}
# 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-4 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 4 separate unadjusted pairwise tests.
COHORT_LEVELS <- c("Cohort 1", "Cohort 2", "Cohort 3", "Cohort 4")
dunnett_pvals <- function(x, group) {
out <- setNames(rep(NA_real_, 4), 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 4 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 * 4), 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 (excluded)",
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-3, for paired change) ----
ease_t1_score = ease_score(.data[[col_ease_t1]]),
ease_t2_score = ease_score(.data[[col_ease_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_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, incl. excluded Cohort 5 ----
cat("\nCohort counts (Cohort 5 exists in the data but is excluded from all analysis below - too new / not yet validated):\n")##
## Cohort counts (Cohort 5 exists in the data but is excluded from all analysis below - too new / not yet validated):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 (excluded) 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 (Cohort 5 excluded throughout) ----
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_not <- df_status %>% filter(cohort_final == "Not enrolled")
df_enrolled <- bind_rows(df_c1, df_c2, df_c3, df_c4)
df_analysis <- bind_rows(df_enrolled, df_not)
# Cohort 4 is still in its active intervention window - it has no 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 != "Cohort 4")
# ---- 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============================================================
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,
`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)
cells <- c(
Total = fmt_mean_sd(df_analysis[[col]], digits),
`Cohort 1` = value_p_cell(fmt_mean_sd(df_c1[[col]], digits), raw_p["Cohort 1"]),
`Cohort 2` = value_p_cell(fmt_mean_sd(df_c2[[col]], digits), raw_p["Cohort 2"]),
`Cohort 3` = value_p_cell(fmt_mean_sd(df_c3[[col]], digits), raw_p["Cohort 3"]),
`Cohort 4` = value_p_cell(fmt_mean_sd(df_c4[[col]], digits), raw_p["Cohort 4"]),
`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 4 simultaneous comparisons
cat_block1 <- function(label, col, levels_vec) {
raw_p <- cat_pvals_bonf(col, groups1, df_not)
header <- make_row1(label, c(
Total = "",
`Cohort 1` = p_only_cell(raw_p["Cohort 1"]),
`Cohort 2` = p_only_cell(raw_p["Cohort 2"]),
`Cohort 3` = p_only_cell(raw_p["Cohort 3"]),
`Cohort 4` = p_only_cell(raw_p["Cohort 4"]),
`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) {
c(
Total = fmt_np(sum(match_fn(df_enrolled)), nrow(df_enrolled)),
`Cohort 1` = fmt_np(sum(match_fn(df_c1)), nrow(df_c1)),
`Cohort 2` = fmt_np(sum(match_fn(df_c2)), nrow(df_c2)),
`Cohort 3` = fmt_np(sum(match_fn(df_c3)), nrow(df_c3)),
`Cohort 4` = fmt_np(sum(match_fn(df_c4)), nrow(df_c4)),
`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-4 and Not-enrolled followed by Dunnett's test (each cohort vs. the not-enrolled reference; ",
"p-values shown are already multiplicity-adjusted for the 4 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 4 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=", sum(df_status$cohort_final == "Cohort 5 (excluded)"),
") is present in the source data but excluded pending validation."
))Table 1. Baseline Characteristics
| Characteristic | Total | Cohort 1 | Cohort 2 | Cohort 3 | Cohort 4 | Not enrolled |
|---|---|---|---|---|---|---|
| Demographics & Socioeconomic Characteristics | ||||||
| Age, years | 63.88 (15.81) | 59.10 (14.02) p = 0.029* | 55.27 (14.52) p = <0.001* | 58.45 (16.02) p = 0.045* | 60.86 (14.25) p = 0.479 | 64.67 (15.82) |
| Sex | p = 0.408 | p = 1.000 | p = 1.000 | p = 0.804 | ||
| Female | 623 (44.6%) | 33 (55.0%) | 20 (50.0%) | 21 (50.0%) | 20 (55.6%) | 529 (43.4%) |
| Male | 774 (55.4%) | 27 (45.0%) | 20 (50.0%) | 21 (50.0%) | 16 (44.4%) | 690 (56.6%) |
| Race | p = 0.004* | p = 0.097 | p = 0.007* | p = 0.232 | ||
| White | 575 (41.2%) | 14 (23.3%) | 9 (22.5%) | 8 (19.0%) | 9 (25.0%) | 535 (43.9%) |
| Black | 354 (25.3%) | 25 (41.7%) | 14 (35.0%) | 18 (42.9%) | 13 (36.1%) | 284 (23.3%) |
| Other | 468 (33.5%) | 21 (35.0%) | 17 (42.5%) | 16 (38.1%) | 14 (38.9%) | 400 (32.8%) |
| Ethnicity | p = 0.321 | p = 0.392 | p = 1.000 | p = 1.000 | ||
| Not Spanish / Hispanic / Latino | 919 (65.9%) | 36 (60.0%) | 21 (52.5%) | 27 (64.3%) | 21 (58.3%) | 814 (66.9%) |
| Spanish / Hispanic / Latino | 425 (30.5%) | 24 (40.0%) | 18 (45.0%) | 15 (35.7%) | 14 (38.9%) | 354 (29.1%) |
| Declined / Unavailable | 51 (3.7%) | 0 (0.0%) | 1 (2.5%) | 0 (0.0%) | 1 (2.8%) | 49 (4.0%) |
| Insurance | p = 1.000 | p = 0.122 | p = 0.033* | p = 1.000 | ||
| Commercial | 366 (26.2%) | 18 (30.0%) | 10 (25.0%) | 11 (26.2%) | 9 (25.0%) | 318 (26.1%) |
| Medicare | 654 (46.8%) | 23 (38.3%) | 12 (30.0%) | 12 (28.6%) | 14 (38.9%) | 593 (48.6%) |
| Medicaid | 315 (22.5%) | 17 (28.3%) | 17 (42.5%) | 15 (35.7%) | 12 (33.3%) | 254 (20.8%) |
| Dual Eligible | 50 (3.6%) | 2 (3.3%) | 1 (2.5%) | 2 (4.8%) | 1 (2.8%) | 44 (3.6%) |
| Uninsured / Self Pay | 12 (0.9%) | 0 (0.0%) | 0 (0.0%) | 2 (4.8%) | 0 (0.0%) | 10 (0.8%) |
| Clinical Characteristics | ||||||
| HbA1c Value (Last 3 Months) | 9.94 (2.14) | 10.20 (2.10) p = 0.745 | 10.23 (2.08) p = 0.813 | 9.95 (2.00) p = 1.000 | 10.50 (2.68) p = 0.347 | 9.90 (2.13) |
| LACE+ Readmission Score | 58.39 (16.04) | 59.33 (14.58) p = 0.992 | 54.42 (17.68) p = 0.380 | 56.45 (15.95) p = 0.881 | 59.36 (14.67) p = 0.996 | 58.51 (16.10) |
| Baseline Healthcare Utilization (90 days pre-index) | ||||||
| Hospitalizations (pre, 90 days) | 1.52 (1.27) | 0.67 (0.80) p = <0.001* | 1.60 (0.98) p = 1.000 | 1.36 (0.88) p = 0.722 | 1.25 (0.77) p = 0.423 | 1.57 (1.30) |
| ED visits (pre, 90 days) | 1.69 (1.75) | 1.13 (2.17) p = 0.034* | 1.60 (1.10) p = 0.977 | 1.40 (1.06) p = 0.627 | 1.31 (0.95) p = 0.452 | 1.74 (1.78) |
| Attrition / Disenrollment (Enrolled Only) | ||||||
| Disenrollment signal identified (T1/T2 chart-review notes) | 22 (12.4%) | 10 (16.7%) | 8 (20.0%) | 4 (9.5%) | 0 (0.0%) | – |
| Never/stopped receiving deliveries | 10 (5.6%) | 4 (6.7%) | 4 (10.0%) | 2 (4.8%) | 0 (0.0%) | – |
| Address issue | 3 (1.7%) | 2 (3.3%) | 1 (2.5%) | 0 (0.0%) | 0 (0.0%) | – |
| Deceased | 2 (1.1%) | 0 (0.0%) | 1 (2.5%) | 1 (2.4%) | 0 (0.0%) | – |
| Delivery cancelled / hold | 2 (1.1%) | 2 (3.3%) | 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%) | – |
| Unresponsive / disengaged | 2 (1.1%) | 0 (0.0%) | 2 (5.0%) | 0 (0.0%) | 0 (0.0%) | – |
| Requested removal | 1 (0.6%) | 0 (0.0%) | 0 (0.0%) | 1 (2.4%) | 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-4 and Not-enrolled followed by Dunnett's test (each cohort vs. the not-enrolled reference; p-values shown are already multiplicity-adjusted for the 4 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 4 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 present in the source data but excluded pending validation.
============================================================
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
)
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)))
# ---- Secondary: patient experience ----
row_help <- make_row2("Nutrition material helped healthier food choices (Yes)",
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_control <- make_row2("Program helped patient feel more in control (A lot / Somewhat)",
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_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)))
# ---- 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. Neutral"),
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 <- sum(v == "Neutral")
paste0(round(100 * pos / n, 1), "% easy/very easy, ", round(100 * neg / n, 1), "% neutral (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)
experience_block <- rbind(row_help, row_control, row_waste, 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-4). 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 3 \"improved\" summary rows here). ",
"Continuous change rows (HbA1c, ED, IP) are tested with a paired Wilcoxon signed-rank test (H0: median change ",
"= 0); the 3 \"improved\" rows 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. Bold ",
"with * = p < .05 in both cases. Rows reporting a single proportion (reduction ≥1%, food insecurity, ",
"nutrition/control impact, food waste) have no comparable reference group and are not significance-tested. 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. Food ease is scored ",
"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 |
|---|---|---|---|---|---|
| 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 = – |
| 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.34 (1.52), n=178 | 1.13 (2.17), n=60 | 1.60 (1.10), n=40 | 1.40 (1.06), n=42 | 1.31 (0.95), n=36 |
| Hospitalizations at baseline (past 90 days), mean (SD) | 1.16 (0.93), n=178 | 0.67 (0.80), n=60 | 1.60 (0.98), n=40 | 1.36 (0.88), n=42 | 1.25 (0.77), n=36 |
| 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 = – |
| 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 = – |
| Secondary Outcomes: Food Insecurity | |||||
| Food insecurity positive at baseline (T0) | 175/176 (99.4%) | 59/60 (98.3%) | 40/40 (100.0%) | 42/42 (100.0%) | 34/34 (100.0%) |
| Food insecurity positive at T2 | 43/67 (64.2%) | 28/43 (65.1%) | 15/24 (62.5%) | – | – |
| Secondary Outcomes: Patient Experience | |||||
| Nutrition material helped healthier food choices (Yes) | 67/75 (89.3%) | 32/35 (91.4%) | 25/26 (96.2%) | 10/14 (71.4%) | – |
| Program helped patient feel more in control (A lot / Somewhat) | 65/83 (78.3%) | 33/38 (86.8%) | 19/26 (73.1%) | 13/19 (68.4%) | – |
| Reported no food waste (T1 or T2, pooled) | 121/143 (84.6%) | 67/77 (87.0%) | 39/47 (83.0%) | 15/19 (78.9%) | – |
| Food ease improved, paired T1→T2 | 26.7% improved, 20% worsened (n=45) p = 0.664 | 20.7% improved, 27.6% worsened (n=29) p = 0.791 | 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-4). 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 3 "improved" summary rows here). Continuous change rows (HbA1c, ED, IP) are tested with a paired Wilcoxon signed-rank test (H0: median change = 0); the 3 "improved" rows 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. Bold with * = p < .05 in both cases. Rows reporting a single proportion (reduction ≥1%, food insecurity, nutrition/control impact, food waste) have no comparable reference group and are not significance-tested. 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. Food ease is scored 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 3 condensed "improved" rows: the before/after
# prevalence (both the positive response and its opposite) and the complete
# paired Improved/No change/Worsened breakdown with significance, for each
# of the 3 patient-experience measures.
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)
detail_df <- rbind(detail_ease, detail_conf, detail_prepare)
detail_header_rows <- which(detail_df$Outcome %in% c(
"Food ease, paired T1→T2",
"Patient confidence (food choices + labels), paired baseline→T2",
"Confidence in meal preparation, paired T1→T2 (no T0 baseline collected)"
))
render_html_table(detail_df, "Patient Experience Detail: Food Ease and Confidence",
bold_rows = detail_header_rows,
sections = list(
"Food Ease" = 1,
"Patient Confidence (Food Choices + Labels)" = nrow(detail_ease) + 1,
"Confidence in Meal Preparation" = nrow(detail_ease) + nrow(detail_conf) + 1
),
footnote = paste0(
"Supplementary detail for the 3 \"improved\" summary rows in Table 2. \"At [timepoint]: ... vs. ...\" rows 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. Meal preparation confidence was never asked at T0, so it is shown as its own T1-to-T2 comparison."
))Patient Experience Detail: Food Ease and Confidence
| Outcome | Overall Enrolled | Cohort 1 | Cohort 2 | Cohort 3 | Cohort 4 |
|---|---|---|---|---|---|
| Food Ease | |||||
| Food ease at T1: Easy/Very easy vs. Neutral | 92.4% easy/very easy, 7.6% neutral (n=79) | 92.1% easy/very easy, 7.9% neutral (n=38) | 91.7% easy/very easy, 8.3% neutral (n=24) | 94.1% easy/very easy, 5.9% neutral (n=17) | – |
| Food ease at T2: Easy/Very easy vs. Neutral | 96.8% easy/very easy, 3.2% neutral (n=62) | 95% easy/very easy, 5% neutral (n=40) | 100% easy/very easy, 0% neutral (n=22) | – | – |
| Food ease, paired T1→T2 | n=45 p = 0.664 | n=29 p = 0.791 | n=16 p = 0.125 | n=0 p = – | n=0 p = – |
| Improved | 12/45 (26.7%) | 6/29 (20.7%) | 6/16 (37.5%) | – | – |
| No change | 24/45 (53.3%) | 15/29 (51.7%) | 9/16 (56.2%) | – | – |
| Worsened | 9/45 (20.0%) | 8/29 (27.6%) | 1/16 (6.2%) | – | – |
| Patient Confidence (Food Choices + Labels) | |||||
| Confidence (food choices + labels) at baseline (T0): High vs. Low | 91.1% high, 8.9% low (n=123) | 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) |
| 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 = – |
| 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%) | – | – |
| 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 = – |
| 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 3 "improved" summary rows in Table 2. "At [timepoint]: ... vs. ..." rows 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. 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-4; 9/12-month excludes Cohort 4 - 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 | 178 | 10.21 | 2.19 |
| 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 Cohort 4 (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 ",
"Cohort 4 (still in its 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 Cohort 4 (still in its 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-4; 9/12-month excludes Cohort 4 - 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 | 178 | 1.34 | 1.52 |
| 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 Cohort 4 (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 Cohort 4 (still in its 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 Cohort 4 (still in its 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-4; 9/12-month excludes Cohort 4 - 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 | 178 | 1.16 | 0.93 |
| 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 Cohort 4 (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 Cohort 4 (still in its 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 Cohort 4 (still in its active intervention window). |