|
| 1 | +# ────────────────────────────────────────────────────────────────── |
| 2 | +# Overlap Weighting (ATO) — exact overlap-weighted effect estimation |
| 3 | +# POST /analysis/overlap-weighting/run |
| 4 | +# |
| 5 | +# Reuses CohortMethod for data extraction + propensity scoring, then applies |
| 6 | +# ATO overlap weights (w = 1-e for treated, w = e for controls; Li-Morgan- |
| 7 | +# Zaslavsky) and a weighted Cox model. Returns the same normalized shape as |
| 8 | +# /analysis/estimation/run so the caller's mapping is unchanged. Negative-control |
| 9 | +# calibration reuses the shared compute_calibration helper. |
| 10 | +# ────────────────────────────────────────────────────────────────── |
| 11 | + |
| 12 | +library(CohortMethod) |
| 13 | +library(FeatureExtraction) |
| 14 | +library(DatabaseConnector) |
| 15 | +library(survival) |
| 16 | +source("/app/R/connection.R") |
| 17 | +source("/app/R/covariates.R") |
| 18 | +source("/app/R/progress.R") |
| 19 | +source("/app/R/results.R") |
| 20 | +source("/app/R/calibration.R") |
| 21 | + |
| 22 | +# ATO-weighted standardized mean differences, computed directly (CohortMethod's |
| 23 | +# computeCovariateBalance does not apply an arbitrary weights column). Covariates |
| 24 | +# from FeatureExtraction are ~all binary indicators, so SD = sqrt(p(1-p)). Returns |
| 25 | +# the top covariates by pre-weighting imbalance in the caller's balance shape. |
| 26 | +compute_ato_balance <- function(cmData, df) { |
| 27 | + covs <- tryCatch(as.data.frame(dplyr::collect(cmData$covariates)), error = function(e) NULL) |
| 28 | + if (is.null(covs) || nrow(covs) == 0) return(list()) |
| 29 | + covs <- covs[covs$rowId %in% df$rowId, , drop = FALSE] |
| 30 | + if (nrow(covs) == 0) return(list()) |
| 31 | + info <- merge(covs, df[, c("rowId", "treatment", "ato_w")], by = "rowId") |
| 32 | + # Restrict to binary indicator covariates (value == 1); continuous covariates |
| 33 | + # (e.g. age) break the p(1-p) SMD formula and are reported as balanced-by-design. |
| 34 | + info <- info[!is.na(info$covariateValue) & info$covariateValue == 1, , drop = FALSE] |
| 35 | + if (nrow(info) == 0) return(list()) |
| 36 | + n1 <- sum(df$treatment == 1); n0 <- sum(df$treatment == 0) |
| 37 | + W1 <- sum(df$ato_w[df$treatment == 1]); W0 <- sum(df$ato_w[df$treatment == 0]) |
| 38 | + if (n1 == 0 || n0 == 0 || W1 == 0 || W0 == 0) return(list()) |
| 39 | + info$wv <- info$ato_w * info$covariateValue |
| 40 | + # Fast group sums via rowsum (C-level) keyed by "covariateId_treatment". |
| 41 | + key <- paste0(info$covariateId, "_", info$treatment) |
| 42 | + sv <- rowsum(info$covariateValue, key) |
| 43 | + swv <- rowsum(info$wv, key) |
| 44 | + get <- function(m, k) if (k %in% rownames(m)) m[k, 1] else 0 |
| 45 | + covRef <- tryCatch(as.data.frame(dplyr::collect(cmData$covariateRef)), error = function(e) NULL) |
| 46 | + name_of <- function(cid) { |
| 47 | + if (!is.null(covRef) && "covariateName" %in% names(covRef)) { |
| 48 | + nm <- covRef$covariateName[covRef$covariateId == cid] |
| 49 | + if (length(nm) > 0) return(substr(as.character(nm[1]), 1, 120)) |
| 50 | + } |
| 51 | + paste0("covariate ", cid) |
| 52 | + } |
| 53 | + out <- list() |
| 54 | + for (cid in unique(info$covariateId)) { |
| 55 | + k1 <- paste0(cid, "_1"); k0 <- paste0(cid, "_0") |
| 56 | + p1u <- get(sv, k1) / n1; p0u <- get(sv, k0) / n0 |
| 57 | + p1w <- get(swv, k1) / W1; p0w <- get(swv, k0) / W0 |
| 58 | + sdu <- sqrt((p1u * (1 - p1u) + p0u * (1 - p0u)) / 2) |
| 59 | + sdw <- sqrt((p1w * (1 - p1w) + p0w * (1 - p0w)) / 2) |
| 60 | + out[[length(out) + 1]] <- list( |
| 61 | + covariate_name = name_of(cid), |
| 62 | + smd_before = round(if (!is.na(sdu) && sdu > 0) (p1u - p0u) / sdu else 0, 4), |
| 63 | + smd_after = round(if (!is.na(sdw) && sdw > 0) (p1w - p0w) / sdw else 0, 4), |
| 64 | + mean_target_before = round(p1u, 4), mean_comp_before = round(p0u, 4), |
| 65 | + mean_target_after = round(p1w, 4), mean_comp_after = round(p0w, 4) |
| 66 | + ) |
| 67 | + } |
| 68 | + ord <- order(sapply(out, function(x) -abs(x$smd_before))) |
| 69 | + out[ord[seq_len(min(60, length(out)))]] |
| 70 | +} |
| 71 | + |
| 72 | +#* Run overlap-weighted (ATO) population-level estimation |
| 73 | +#* @post /analysis/overlap-weighting/run |
| 74 | +#* @serializer unboxedJSON |
| 75 | +function(body, response) { |
| 76 | + spec <- body |
| 77 | + logger <- create_analysis_logger() |
| 78 | + |
| 79 | + if (is.null(spec)) { |
| 80 | + response$status <- 400L |
| 81 | + return(list(status = "error", message = "No specification provided in request body")) |
| 82 | + } |
| 83 | + missing <- setdiff(c("source", "cohorts", "model"), names(spec)) |
| 84 | + if (length(missing) > 0) { |
| 85 | + response$status <- 400L |
| 86 | + return(list(status = "error", message = paste("Missing required fields:", paste(missing, collapse = ", ")))) |
| 87 | + } |
| 88 | + |
| 89 | + safe_execute(response, logger, { |
| 90 | + connectionDetails <- create_hades_connection(spec$source) |
| 91 | + connection <- connect_with_retry(connectionDetails) |
| 92 | + on.exit(safe_disconnect(connection), add = TRUE) |
| 93 | + |
| 94 | + cdmSchema <- spec$source$cdm_schema |
| 95 | + vocabSchema <- spec$source$vocab_schema %||% cdmSchema |
| 96 | + resultsSchema <- spec$source$results_schema |
| 97 | + |
| 98 | + targetId <- as.integer(spec$cohorts$target_cohort_id) |
| 99 | + comparatorId <- as.integer(spec$cohorts$comparator_cohort_id) |
| 100 | + outcomeIds <- as.integer(spec$cohorts$outcome_cohort_ids) |
| 101 | + outcomeNames <- spec$cohorts$outcome_names %||% list() |
| 102 | + ncOutcomeIds <- as.integer(spec$negative_control_outcomes %||% spec$negativeControlOutcomes %||% list()) |
| 103 | + extractOutcomeIds <- unique(c(outcomeIds, ncOutcomeIds)) |
| 104 | + |
| 105 | + covariateSettings <- build_covariate_settings(spec$covariate_settings) |
| 106 | + dataArgs <- CohortMethod::createGetDbCohortMethodDataArgs(covariateSettings = covariateSettings) |
| 107 | + logger$info("Extracting CohortMethod data (ATO)") |
| 108 | + cmData <- CohortMethod::getDbCohortMethodData( |
| 109 | + connectionDetails = connectionDetails, |
| 110 | + cdmDatabaseSchema = cdmSchema, |
| 111 | + targetId = targetId, |
| 112 | + comparatorId = comparatorId, |
| 113 | + outcomeIds = extractOutcomeIds, |
| 114 | + exposureDatabaseSchema = resultsSchema, |
| 115 | + exposureTable = "cohort", |
| 116 | + outcomeDatabaseSchema = resultsSchema, |
| 117 | + outcomeTable = "cohort", |
| 118 | + getDbCohortMethodDataArgs = dataArgs |
| 119 | + ) |
| 120 | + |
| 121 | + tar_start <- as.integer(spec$model$time_at_risk_start %||% spec$model$timeAtRiskStart %||% 1) |
| 122 | + tar_end <- as.integer(spec$model$time_at_risk_end %||% spec$model$timeAtRiskEnd %||% 9999) |
| 123 | + end_anchor <- spec$model$end_anchor %||% spec$model$endAnchor %||% "cohort end" |
| 124 | + |
| 125 | + popArgs <- CohortMethod::createCreateStudyPopulationArgs( |
| 126 | + removeSubjectsWithPriorOutcome = TRUE, |
| 127 | + riskWindowStart = tar_start, startAnchor = "cohort start", |
| 128 | + riskWindowEnd = tar_end, endAnchor = end_anchor, minDaysAtRisk = 1 |
| 129 | + ) |
| 130 | + psArgs <- CohortMethod::createCreatePsArgs( |
| 131 | + maxCohortSizeForFitting = 250000, errorOnHighCorrelation = FALSE, stopOnError = FALSE |
| 132 | + ) |
| 133 | + |
| 134 | + # Fit one ATO-weighted Cox for a given outcome id; returns the model + the |
| 135 | + # PS object (for diagnostics) + the weighted population data frame. |
| 136 | + ato_fit <- function(oid) { |
| 137 | + pop <- CohortMethod::createStudyPopulation( |
| 138 | + cohortMethodData = cmData, population = NULL, outcomeId = oid, |
| 139 | + createStudyPopulationArgs = popArgs |
| 140 | + ) |
| 141 | + pdf <- as.data.frame(pop) |
| 142 | + if (sum(pdf$treatment == 1) < 10 || sum(pdf$treatment == 0) < 10) { |
| 143 | + return(list(ok = FALSE, df = pdf)) |
| 144 | + } |
| 145 | + ps <- CohortMethod::createPs(cohortMethodData = cmData, population = pop, createPsArgs = psArgs) |
| 146 | + df <- as.data.frame(ps) |
| 147 | + e <- pmin(pmax(df$propensityScore, 1e-6), 1 - 1e-6) |
| 148 | + df$ato_w <- ifelse(df$treatment == 1, 1 - e, e) |
| 149 | + df$time <- if ("survivalTime" %in% names(df)) df$survivalTime else df$timeAtRisk |
| 150 | + df$event <- as.integer(df$outcomeCount > 0) |
| 151 | + df <- df[!is.na(df$time) & df$time > 0, , drop = FALSE] |
| 152 | + m <- tryCatch( |
| 153 | + survival::coxph(survival::Surv(time, event) ~ treatment, data = df, weights = ato_w, robust = TRUE), |
| 154 | + error = function(e) NULL |
| 155 | + ) |
| 156 | + list(ok = TRUE, ps = ps, df = df, model = m) |
| 157 | + } |
| 158 | + |
| 159 | + estimates_list <- list() |
| 160 | + ps_auc <- NA_real_; equipoise_val <- NA_real_; ps_dist_data <- NULL |
| 161 | + balance_all <- NULL; n_target <- NA_integer_; n_comparator <- NA_integer_ |
| 162 | + |
| 163 | + for (oid in outcomeIds) { |
| 164 | + logger$info(sprintf("ATO outcome %d", oid)) |
| 165 | + fit <- ato_fit(oid) |
| 166 | + oname <- outcomeNames[[as.character(oid)]] %||% sprintf("Outcome %d", oid) |
| 167 | + if (!isTRUE(fit$ok) || is.null(fit$model)) { |
| 168 | + estimates_list[[length(estimates_list) + 1]] <- list( |
| 169 | + outcome_id = oid, outcome_name = oname, hazard_ratio = NA, ci_95_lower = NA, |
| 170 | + ci_95_upper = NA, p_value = NA, log_hr = NA, se_log_hr = NA, |
| 171 | + target_outcomes = as.integer(sum(fit$df$outcomeCount[fit$df$treatment == 1] > 0)), |
| 172 | + comparator_outcomes = as.integer(sum(fit$df$outcomeCount[fit$df$treatment == 0] > 0)), |
| 173 | + warning = "Insufficient subjects or model did not converge" |
| 174 | + ) |
| 175 | + next |
| 176 | + } |
| 177 | + df <- fit$df |
| 178 | + if (is.na(ps_auc)) { |
| 179 | + ps_auc <- tryCatch(CohortMethod::computePsAuc(fit$ps), error = function(e) NA_real_) |
| 180 | + equipoise_val <- tryCatch(CohortMethod::computeEquipoise(fit$ps), error = function(e) NA_real_) |
| 181 | + ps_dist_data <- tryCatch(extract_ps_distribution(fit$ps), error = function(e) NULL) |
| 182 | + n_target <- as.integer(sum(df$treatment == 1)) |
| 183 | + n_comparator <- as.integer(sum(df$treatment == 0)) |
| 184 | + balance_all <- tryCatch( |
| 185 | + compute_ato_balance(cmData, df), |
| 186 | + error = function(e) { logger$warn(paste("ATO balance failed:", e$message)); list() } |
| 187 | + ) |
| 188 | + } |
| 189 | + cf <- tryCatch(coef(fit$model)[["treatment"]], error = function(e) NA_real_) |
| 190 | + se <- tryCatch(sqrt(diag(vcov(fit$model)))[["treatment"]], error = function(e) NA_real_) |
| 191 | + hr <- exp(cf) |
| 192 | + lo <- exp(cf - 1.96 * se); hi <- exp(cf + 1.96 * se) |
| 193 | + pv <- if (!is.na(cf) && !is.na(se) && se > 0) 2 * pnorm(-abs(cf / se)) else NA_real_ |
| 194 | + estimates_list[[length(estimates_list) + 1]] <- list( |
| 195 | + outcome_id = oid, outcome_name = oname, |
| 196 | + hazard_ratio = round(hr, 4), ci_95_lower = round(lo, 4), ci_95_upper = round(hi, 4), |
| 197 | + p_value = round(pv, 6), log_hr = round(cf, 4), se_log_hr = round(se, 4), |
| 198 | + target_outcomes = as.integer(sum(df$event[df$treatment == 1])), |
| 199 | + comparator_outcomes = as.integer(sum(df$event[df$treatment == 0])) |
| 200 | + ) |
| 201 | + logger$info(sprintf("ATO outcome %d: HR=%.3f [%.3f, %.3f]", oid, hr, lo, hi)) |
| 202 | + } |
| 203 | + |
| 204 | + # Negative controls (ATO-weighted Cox) for empirical calibration. |
| 205 | + nc_estimates <- list() |
| 206 | + for (nc_id in ncOutcomeIds) { |
| 207 | + tryCatch({ |
| 208 | + fit <- ato_fit(nc_id) |
| 209 | + if (isTRUE(fit$ok) && !is.null(fit$model)) { |
| 210 | + cf <- coef(fit$model)[["treatment"]] |
| 211 | + se <- sqrt(diag(vcov(fit$model)))[["treatment"]] |
| 212 | + if (is.finite(cf) && is.finite(se)) { |
| 213 | + nc_estimates[[length(nc_estimates) + 1]] <- list( |
| 214 | + outcome_id = nc_id, log_rr = round(cf, 4), se_log_rr = round(se, 4) |
| 215 | + ) |
| 216 | + } |
| 217 | + } |
| 218 | + }, error = function(e) logger$warn(sprintf("NC %d failed: %s", nc_id, e$message))) |
| 219 | + } |
| 220 | + calibration_data <- if (length(nc_estimates) > 0) { |
| 221 | + tryCatch(compute_calibration(estimates_list, nc_estimates), error = function(e) NULL) |
| 222 | + } else NULL |
| 223 | + |
| 224 | + balance_summary <- if (is.null(balance_all)) list() else balance_all |
| 225 | + max_smd_after <- if (length(balance_summary) > 0) { |
| 226 | + round(max(sapply(balance_summary, function(x) abs(x$smd_after)), na.rm = TRUE), 4) |
| 227 | + } else NA_real_ |
| 228 | + max_smd_before <- if (length(balance_summary) > 0) { |
| 229 | + round(max(sapply(balance_summary, function(x) abs(x$smd_before)), na.rm = TRUE), 4) |
| 230 | + } else NA_real_ |
| 231 | + |
| 232 | + list( |
| 233 | + status = "completed", |
| 234 | + method = "ATO overlap weighting (weighted Cox)", |
| 235 | + summary = list(target_count = n_target, comparator_count = n_comparator, outcome_counts = list()), |
| 236 | + estimates = estimates_list, |
| 237 | + propensity_score = list( |
| 238 | + auc = round(ps_auc, 4), equipoise = round(equipoise_val, 4), |
| 239 | + max_smd_after = max_smd_after, max_smd_before = max_smd_before, |
| 240 | + distribution = ps_dist_data |
| 241 | + ), |
| 242 | + calibration = calibration_data, |
| 243 | + covariate_balance = balance_summary, |
| 244 | + negative_controls = list(estimates = nc_estimates) |
| 245 | + ) |
| 246 | + }) |
| 247 | +} |
0 commit comments