hinna0818/ukbanalytica_v2/ukbsci-propensity

ukbsci-propensity — Propensity scores, matching, and IPTW on UKB cohorts

0. RAP guardrails

Strict local-agent boundary: this skill is for script generation, workflow planning, package guidance, and interpretation of aggregate outputs. The agent must not read, inspect, summarize, or process real UK Biobank RAP participant-level data, including de-identified row-level tables, raw RAP fields, exact dates, per-row predictions, row-level SHAP matrices, screenshots, tracebacks, or logs containing row-level values. Generate scripts for the user to run inside RAP; only aggregate results or rendered figures may be shared back with the agent. See ../references/agent-privacy-boundary.md.

The matched / weighted dataset is participant-level and stays in RAP project storage. Shareable artefacts are aggregate: balance tables, Love plots, PS distribution figures, and final weighted-effect estimates.

1. When to load

• Estimate a propensity score for a binary treatment / exposure.
• Build a 1:N PSM-matched sub-cohort.
• Construct IPTW weights (ATE / ATT / ATC, stabilized, optionally truncated).
• Assess covariate balance before vs after matching / weighting.
• Run a weighted Cox / logistic / linear analysis with robust SEs.

2. When NOT to load

• Vanilla unadjusted / multivariable Cox → ukbsci-regression.
• Multiple-imputation pooling of weighted estimates → use this skill plus

ukbsci-imputation.

• Mediation analysis → ukbsci-mediation.

3. Prerequisites

library(UKBAnalytica)
library(data.table)
 
# Suggested packages
# - MatchIt for matching
# - gbm for gradient-boosted PS
# - sandwich, lmtest for HC0 robust SEs
# - survival for weighted Cox
 
stopifnot("treatment" %in% colnames(cohort))
stopifnot(all(cohort$treatment %in% c(0, 1, NA)))

4. Pipeline (5 phases)

Phase 1 — Estimate PS

covars_ps <- c("age", "sex", "bmi", "smoking_status",
               "Hypertension_history", "Diabetes_history",
               "ldl_cholesterol", "sbp")
 
ps_data <- estimate_propensity_score(
  data       = cohort,
  treatment  = "treatment",
  covariates = covars_ps,
  method     = "logistic"             # "logistic" or "gbm"
)
# Adds column `ps`, bounded to [0.001, 0.999]

For non-linear / non-additive relationships, prefer method = "gbm".

Phase 2 — Choose matching OR weighting

Matching path:

matched <- match_propensity(
  data        = ps_data,
  ps_col      = "ps",
  treatment   = "treatment",
  ratio       = 1,
  caliper     = 0.2,
  method      = "nearest",            # or "optimal"
  replace     = FALSE,
  exact_match = c("sex")              # optional exact-match strata
)
# Adds `match_id`, `match_distance`

Weighting path:

weighted <- calculate_weights(
  data         = ps_data,
  ps_col       = "ps",
  treatment    = "treatment",
  weight_type  = "ATE",               # "ATE", "ATT", or "ATC"
  stabilized   = TRUE,
  truncate     = c(0.01, 0.99)        # weight quantile clipping
)
# Adds `weight`

Phase 3 — Assess balance (before / after)

b_pre  <- assess_balance(
  data       = cohort,
  treatment  = "treatment",
  covariates = covars_ps,
  method     = "unmatched"
)
 
b_post_m <- assess_balance(
  data       = matched,
  treatment  = "treatment",
  covariates = covars_ps,
  method     = "matched"
)
 
b_post_w <- assess_balance(
  data       = weighted,
  treatment  = "treatment",
  covariates = covars_ps,
  method     = "weighted",
  weight_col = "weight",
  threshold  = 0.1
)
 
p_love <- plot_balance(
  balance_before = b_pre,
  balance_after  = b_post_m,             # or b_post_w for IPTW
  threshold      = 0.1,
  title          = "Love plot — before vs. after matching"
)
 
p_ps <- plot_ps_distribution(
  data      = ps_data,
  ps_col    = "ps",
  treatment = "treatment",
  type      = "mirror"               # "histogram", "density", "mirror"
)

The standard target: SMD ≤ 0.1 on every covariate post-adjustment.

Phase 4 — Weighted / matched analysis

# Cox under matched cohort (no weighting; matching already balanced)
cox_matched <- runmulti_cox(
  data       = matched,
  main_var   = "treatment",
  covariates = NULL,
  endpoint   = c("outcome_surv_time", "outcome_status")
)
 
# Or IPTW-weighted Cox with robust SE
wt_cox <- run_weighted_analysis(
  data       = weighted,
  exposure   = "treatment",
  covariates = NULL,
  weight_col = "weight",
  model_type = "cox",
  endpoint   = c("outcome_surv_time", "outcome_status"),
  robust_se  = TRUE
)

run_weighted_analysis() also supports model_type = "logistic" / "linear". Robust standard errors via sandwich + lmtest.

Phase 5 — Export

All output goes to /mnt/project/<area>/04-results/:

fwrite(b_pre,    "04-results/05a-balance_unmatched.csv")
fwrite(b_post_m, "04-results/05b-balance_matched.csv")
fwrite(b_post_w, "04-results/05c-balance_weighted.csv")
fwrite(wt_cox,   "04-results/02d-cox_iptw.csv")
 
ggplot2::ggsave("05-figs/Fig07-love_plot.pdf", p_love, width = 7, height = 6)
ggplot2::ggsave("05-figs/Fig08-ps_distribution.pdf", p_ps, width = 7, height = 5)

5. Common pitfalls

1. Bounded PS at `[0.001, 0.999]`. Prevents extreme weights but

shrinks the tails. Verify with plot_ps_distribution() that no group has piled-up boundary mass.

1. Caliper too tight. Default 0.2 × SD(PS) is generous; tightening

below 0.05 × SD can drop most of the treated arm.

1. Replace = TRUE. Increases effective sample but inflates variance;

warn the user and confirm.

1. `weight_type = "ATE"` vs `"ATT"` changes the estimand. PSM by default

targets ATT (the matched cohort). Mixing in IPTW with ATE and reporting them side-by-side without flagging is misleading.

1. Truncation defaults at 1st/99th percentile. Aggressive — if the user

wants to keep extremes, set truncate = c(0, 1).

1. SMD on factors. assess_balance() converts factors to numeric for

SMD; multi-level factors are summarized as a single SMD. For per-level SMD, expand the factor into dummies first.

1. Robust SE for Cox. `run_weighted_analysis(model_type = "cox",

robust_se = TRUE) uses survival::coxph(..., robust = TRUE) semantics — appropriate for IPTW. For matched Cox, prefer **cluster-robust SE on match_id`** (build the formula manually).

1. Missing `weight` column for `method = "weighted"` triggers a clear

error from assess_balance(). Pass weight_col = "weight" explicitly.

6. Key functions

See `references/functions.md` for full argument tables.

7. Related skills

8. References

• `references/functions.md`
• `references/examples.md` — three examples

(PSM with caliper + exact match, IPTW with stabilized ATE weights, weighted Cox with robust SE)