Where-matrix imputation, weighted Tables 1a/1b, and protocol v0.3.1

.github/workflows/protocol-version.yml

@@ -0,0 +1,48 @@
+name: Protocol version check
+
+# The study protocol is a prespecified governance document: any change to
+# docs/protocol/ must be recorded in the version history, not edited silently.
+# This check fails a PR that touches protocol files without bumping
+# version-summary.version in full-protocol.qmd.
+
+on:
+  pull_request:
+    paths:
+      - "docs/protocol/**"
+
+jobs:
+  version-bumped:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+
+      - name: Require a version bump when protocol files change
+        env:
+          BASE_REF: ${{ github.base_ref }}
+        run: |
+          base="origin/${BASE_REF}"
+          git fetch origin "${BASE_REF}" --depth=1
+
+          changed=$(git diff --name-only "$base"...HEAD -- docs/protocol/)
+          echo "Protocol files changed:"
+          echo "$changed"
+
+          get_version() {
+            git show "$1:docs/protocol/full-protocol.qmd" 2>/dev/null \
+              | awk '/^version-summary:/{f=1; next} f && /version:/{gsub(/[" ]/, "", $2); print $2; exit}'
+          }
+          base_version=$(get_version "$base")
+          head_version=$(get_version HEAD)
+          echo "Base version: ${base_version:-<none>} | Head version: ${head_version:-<none>}"
+
+          if [ -z "$head_version" ]; then
+            echo "::error::Could not parse version-summary.version from full-protocol.qmd"
+            exit 1
+          fi
+          if [ "$base_version" = "$head_version" ]; then
+            echo "::error::docs/protocol/ changed but version-summary.version is still ${head_version}. Bump the version and add a version-history entry describing the change."
+            exit 1
+          fi
+          echo "Version bumped: ${base_version:-<none>} -> ${head_version}"

CLAUDE.md

@@ -133,10 +133,10 @@ Roles are comma-separated in `cshm-variables.csv`. A variable may carry multiple
 |------|-------|---------|
 | `design` | Survey design | Survey infrastructure (SurveyCycle, WTS_M) |
 | `intermediate` | Harmonization | Raw cchsflow input needed to derive a unified variable; not used directly by pipeline code |
-| `predictor` | Model | Covariate in the APC model or descriptive analysis |
+| `predictor` | Model | Covariate in the APC model |
 | `model-stratifier` | Model | Stratifies APC into separate fits (e.g. DHH_SEX) |
-| `table1` | Descriptive | Row in Table 1 descriptive statistics |
-| `table1-stratifier` | Descriptive | Stratifies Table 1 columns |
+| `table1` | Descriptive | Row in Table 1 (drives row selection in `get_cshm_desc_data()`) |
+| `table1-stratifier` | Descriptive | Reserved for cycle/extra stratification of descriptive tables (not yet consumed by code) |
 | `apc-numerator` | APC data prep | Defines the event indicator in Stage 7 |
 | `apc-denominator` | APC data prep | Constructs the at-risk person-year denominator in Stage 7 |
 | `imputation-predictor` | Imputation | Included in MICE imputation model |
@@ -148,6 +148,10 @@ Role vocabulary (single source of truth): [schemas/cshm-variables.yaml](schemas/
 
 `haven::tagged_na()` throughout: **NA(a)** = not applicable · **NA(b)** = don't know/refused · **NA(c)** = not asked this cycle
 
+## Protocol versioning
+
+The study protocol is prespecified: any change under `docs/protocol/` (including Appendix D) must bump `version-summary.version` in [docs/protocol/full-protocol.qmd](docs/protocol/full-protocol.qmd) and add a dated `version-history` entry describing the change — in the same commit. Enforced on PRs by `.github/workflows/protocol-version.yml`.
+
 ## Code style
 
 - Follow tidyverse design principles; snake_case for all function and variable names

R/create-descriptive-tables.R

@@ -5,30 +5,58 @@
 NA_c_label <- "Missing from survey"
 
 categorical_predictor_footnote <- paste0(
-  "For categorical variables the values are displayed as N (percent). ",
+  "For categorical variables the values are displayed as unweighted N ",
+  "(survey-weighted percent where weights were supplied; otherwise unweighted percent). ",
   "Percents may not sum to 100 due to missingness."
 )
 continuous_predictor_footnote <- paste0(
-  "For continuous predictors, the values are displayed as min - max, median (IQR)."
+  "For continuous predictors, the values are displayed as min - max, median (IQR); ",
+  "median and IQR are survey-weighted where weights were supplied."
 )
 
 # ---- Formatting helpers -----------------------------------------------------
+# When the descriptive data carry weighted columns (weight_var supplied to the
+# engine), cells show unweighted n with weighted percent / weighted median (IQR)
+# per protocol v0.3.0 §3.4.1. Without weights they fall back to unweighted stats.
 
 format_cat_descriptive_data <- function(descriptive_data_row) {
-  if (is.na(descriptive_data_row[1, "n"]) || descriptive_data_row[1, "n"] == 0) {
+  if (nrow(descriptive_data_row) == 0 || is.na(descriptive_data_row[1, "n"])) {
+    stop("Descriptive table lookup matched no engine row — worksheet/data ",
+         "category mismatch (wiring bug), not an empty stratum.", call. = FALSE)
+  }
+  if (descriptive_data_row[1, "n"] == 0) {
     return("No data")
   }
-  formatted_n <- format(descriptive_data_row[1, "n"], big.mark = ",")
-  paste0(formatted_n, "\n (", round(descriptive_data_row[1, "percent"] * 100, 1), ")")
+  pct <- if ("wtd_percent" %in% colnames(descriptive_data_row) &&
+             !is.na(descriptive_data_row[1, "wtd_percent"])) {
+    descriptive_data_row[1, "wtd_percent"]
+  } else {
+    descriptive_data_row[1, "percent"]
+  }
+  # MI-averaged n can be fractional; display as a rounded count
+  formatted_n <- format(round(descriptive_data_row[1, "n"]), big.mark = ",")
+  paste0(formatted_n, "\n (", round(pct * 100, 1), ")")
 }
 
 format_cont_descriptive_data <- function(descriptive_data_row) {
-  if (descriptive_data_row[1, "n"] == 0) return("No data")
+  if (nrow(descriptive_data_row) == 0 || is.na(descriptive_data_row[1, "n"])) {
+    stop("Descriptive table lookup matched no engine row — worksheet/data ",
+         "mismatch (wiring bug), not an empty stratum.", call. = FALSE)
+  }
+  if (descriptive_data_row[1, "n"] == 0) {
+    return("No data")
+  }
+  weighted <- "wtd_median" %in% colnames(descriptive_data_row) &&
+    !is.na(descriptive_data_row[1, "wtd_median"])
+  med <- if (weighted) descriptive_data_row[1, "wtd_median"]
+         else descriptive_data_row[1, "median"]
+  p25 <- if (weighted) descriptive_data_row[1, "wtd_percentile25"]
+         else descriptive_data_row[1, "percentile25"]
+  p75 <- if (weighted) descriptive_data_row[1, "wtd_percentile75"]
+         else descriptive_data_row[1, "percentile75"]
   paste0(
     descriptive_data_row[1, "min"], " - ", descriptive_data_row[1, "max"], ",\n",
-    descriptive_data_row[1, "median"],
-    " (", descriptive_data_row[1, "percentile25"],
-    " - ", descriptive_data_row[1, "percentile75"], ")"
+    med, " (", p25, " - ", p75, ")"
   )
 }
 
@@ -386,6 +414,7 @@ create_cycle_specific_descriptive_table <- function(
   cycle_col,
   cycle_labels,
   column_stratifier = NULL,
+  weight_var        = NULL,
   sections_order    = NULL,
   include_na        = TRUE
 ) {
@@ -409,7 +438,7 @@ create_cycle_specific_descriptive_table <- function(
 
     cycle_desc <- get_descriptive_data(
       cycle_data, variables_sheet, variable_details_sheet,
-      variables, stratify_config
+      variables, stratify_config, weight_var = weight_var
     )
     cycle_tables[[key]] <- .build_descriptive_table_data(
       cycle_desc, variables_sheet, variable_details_sheet,

R/descriptive-data.R

@@ -1,19 +1,30 @@
 # descriptive-data.R
 # CSHM-specific wrapper for calculating descriptive statistics.
 # Pipeline targets: table_1a_data, table_1b_data
+#
+# Presentation per protocol v0.3.0 §3.4.1: single table with unweighted n and
+# survey-weighted statistics (weighted % for categories and NA-type rows;
+# weighted median/IQR for continuous). Table 1b averages across the m
+# completed imputation datasets.
 
 #' Calculate descriptive statistics for the CSHM study population
 #'
-#' Computes statistics for all predictor variables stratified by model-stratifier.
-#' No row stratification is applied in the base table.
+#' Computes statistics for all table1-role variables stratified by
+#' model-stratifier. No row stratification is applied in the base table.
 #'
 #' @param data Cleaned or imputed study data frame
 #' @param variables_sheet Variables worksheet data frame
 #' @param variable_details_sheet Variable details worksheet data frame
+#' @param weight_var Survey weight column name (e.g. "WTS_M"). NULL = unweighted.
 #' @return Data frame of descriptive statistics (input to create_descriptive_table)
-get_cshm_desc_data <- function(data, variables_sheet, variable_details_sheet) {
-  predictor_vars <- select_vars_by_role("predictor", variables_sheet)
+get_cshm_desc_data <- function(data, variables_sheet, variable_details_sheet,
+                               weight_var = NULL) {
+  # Table rows are selected by the table1 role (the documented wiring);
+  # columns are stratified by model-stratifier (sex). The cycle-specific
+  # appendix table is wired from config (survey_var cycle/sex), not roles.
+  predictor_vars <- select_vars_by_role("table1", variables_sheet)
   sex_stratifier <- select_vars_by_role("model-stratifier", variables_sheet)[1]
+  stopifnot(!is.na(sex_stratifier))
 
   # Only describe variables that were actually harmonized into data
   # (some variables in the sheet may be absent if no variable_details rows matched)
@@ -32,6 +43,53 @@ get_cshm_desc_data <- function(data, variables_sheet, variable_details_sheet) {
     variables_sheet         = variables_sheet,
     variables_details_sheet = variable_details_sheet,
     variables               = available,
-    stratify_config         = stratify_config
+    stratify_config         = stratify_config,
+    weight_var              = weight_var
   )
 }
+
+#' Descriptive statistics averaged across multiple imputations (Table 1b)
+#'
+#' Computes the descriptive statistics within each completed dataset and
+#' averages the estimates across imputations (means for central statistics
+#' and proportions; min of minima, max of maxima). With m = 1 this reduces
+#' to a single get_cshm_desc_data() call.
+#'
+#' @param imputation_result Output of impute_data() (list with $datasets)
+#' @param variables_sheet Variables worksheet data frame
+#' @param variable_details_sheet Variable details worksheet data frame
+#' @param weight_var Survey weight column name. NULL = unweighted.
+#' @return Data frame of descriptive statistics, averaged across imputations
+get_cshm_desc_data_mi <- function(imputation_result, variables_sheet,
+                                  variable_details_sheet, weight_var = NULL) {
+  datasets <- imputation_result$datasets
+  stopifnot(length(datasets) >= 1)
+
+  per_imp <- lapply(datasets, function(d) {
+    get_cshm_desc_data(d, variables_sheet, variable_details_sheet, weight_var)
+  })
+  if (length(per_imp) == 1) return(per_imp[[1]])
+
+  stacked <- dplyr::bind_rows(per_imp, .id = ".imp")
+  # All imputations must produce identical row sets (worksheet-driven rows;
+  # factor levels preserved by the write-back) — averaging mismatched groups
+  # would be silent corruption.
+  stopifnot(nrow(stacked) == nrow(per_imp[[1]]) * length(per_imp))
+
+  mean_cols <- intersect(
+    c("median", "percentile25", "percentile75", "n", "percent",
+      "wtd_percentile25", "wtd_median", "wtd_percentile75", "wtd_percent"),
+    colnames(stacked)
+  )
+  key_cols <- setdiff(colnames(stacked), c(".imp", mean_cols, "min", "max"))
+
+  stacked |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(key_cols))) |>
+    dplyr::summarise(
+      dplyr::across(dplyr::all_of(mean_cols), ~ mean(.x, na.rm = FALSE)),
+      min = if (all(is.na(min))) NA else min(min, na.rm = TRUE),
+      max = if (all(is.na(max))) NA else max(max, na.rm = TRUE),
+      .groups = "drop"
+    ) |>
+    as.data.frame()
+}

R/get-descriptive-data.R

@@ -1,14 +1,57 @@
 # get-descriptive-data.R
 # Calculate descriptive statistics for the study population.
-# Ported from DemPoRT-V2-dev (origin/dev).
+# Ported from DemPoRT-V2-dev (origin/dev); extended with survey-weighted
+# statistics (weight_var) per protocol v0.3.0 §3.4.1 / Appendix D.
+
+#' Weighted quantile (midpoint-ECDF with linear interpolation)
+#'
+#' Interpolates linearly between order statistics positioned at the midpoints
+#' of their cumulative weight intervals (the standard weighted analogue of a
+#' continuous sample quantile). With equal weights the weighted median equals
+#' the type-7 sample median; other quantiles are close but not identical to
+#' type 7. Values beyond the first/last midpoint take the boundary value.
+#'
+#' @param x Numeric vector (NAs removed by caller)
+#' @param w Numeric weights aligned with x
+#' @param probs Probabilities in [0, 1]
+#' @return Numeric vector of weighted quantiles
+weighted_quantile <- function(x, w, probs) {
+  if (anyNA(w) || any(w <= 0)) {
+    stop("weighted_quantile: weights must be positive and non-missing ",
+         "(found ", sum(is.na(w)), " NA and ", sum(!is.na(w) & w <= 0),
+         " non-positive of ", length(w), ").", call. = FALSE)
+  }
+  if (length(x) == 0) return(rep(NA_real_, length(probs)))
+  if (length(x) == 1) return(rep(x, length(probs)))
+  ord <- order(x)
+  x <- x[ord]; w <- w[ord]
+  cw <- cumsum(w)
+  midpoints <- (cw - w / 2) / sum(w)
+  stats::approx(midpoints, x, xout = probs, rule = 2, ties = "ordered")$y
+}
 
 get_descriptive_data <- function(
   data,
   variables_sheet,
   variables_details_sheet,
   variables,
-  stratify_config
+  stratify_config,
+  weight_var = NULL
 ) {
+  use_weights <- !is.null(weight_var)
+  if (use_weights) {
+    if (!weight_var %in% colnames(data)) {
+      stop("weight_var '", weight_var, "' not found in data.")
+    }
+    w_all <- data[[weight_var]]
+    if (anyNA(w_all) || any(w_all <= 0, na.rm = TRUE)) {
+      stop("Survey weight '", weight_var, "' has ", sum(is.na(w_all)),
+           " missing and ", sum(!is.na(w_all) & w_all <= 0),
+           " non-positive values — weighted statistics would silently ",
+           "renormalize over a different population than the displayed n. ",
+           "Resolve upstream (harmonization/cleaning) first.", call. = FALSE)
+    }
+  }
   descriptive_data <- data.frame(
     variable   = c(),
     cat        = c(),
@@ -60,14 +103,23 @@ get_descriptive_data <- function(
               current_stratifier_info$stratifier_combination[[strat]][1]
           }
           vals <- current_stratifier_info$data[[variable]]
-          vals <- vals[!is.na(vals)]
+          keep <- !is.na(vals)
+          vals <- vals[keep]
           s    <- summary(vals)
           new_row$median       <- s[[3]]
           new_row$percentile25 <- s[[2]]
           new_row$percentile75 <- s[[5]]
           new_row$min          <- s[[1]]
           new_row$max          <- s[[6]]
           new_row$n            <- length(vals)
+          if (use_weights) {
+            w <- current_stratifier_info$data[[weight_var]][keep]
+            wq <- weighted_quantile(vals, w, c(0.25, 0.5, 0.75))
+            new_row$wtd_percentile25 <- wq[1]
+            new_row$wtd_median       <- wq[2]
+            new_row$wtd_percentile75 <- wq[3]
+            new_row$wtd_percent      <- NA_real_
+          }
           descriptive_data <<- rbind(descriptive_data, new_row)
         }
       )
@@ -99,6 +151,13 @@ get_descriptive_data <- function(
           )
           new_row$n       <- nrow(filtered)
           new_row$percent <- new_row$n / nrow(current_stratifier_info$data)
+          if (use_weights) {
+            new_row$wtd_percentile25 <- NA_real_
+            new_row$wtd_median       <- NA_real_
+            new_row$wtd_percentile75 <- NA_real_
+            new_row$wtd_percent <- sum(filtered[[weight_var]], na.rm = TRUE) /
+              sum(current_stratifier_info$data[[weight_var]], na.rm = TRUE)
+          }
           descriptive_data <<- rbind(descriptive_data, new_row)
         }
       )
@@ -135,6 +194,13 @@ get_descriptive_data <- function(
             )
             new_row$n       <- nrow(filtered)
             new_row$percent <- new_row$n / nrow(current_stratifier_info$data)
+            if (use_weights) {
+              new_row$wtd_percentile25 <- NA_real_
+              new_row$wtd_median       <- NA_real_
+              new_row$wtd_percentile75 <- NA_real_
+              new_row$wtd_percent <- sum(filtered[[weight_var]], na.rm = TRUE) /
+                sum(current_stratifier_info$data[[weight_var]], na.rm = TRUE)
+            }
             descriptive_data <<- rbind(descriptive_data, new_row)
           }
         )