ถ่วงน้ำหนักแล้วยังไม่สมดุล: แก้แบบจำลองก่อนจะโทษข้อมูล

# Simulated hip-fracture registry of older adults: surgery within 24 hours (surg24), delirium, one-year
# death, a small nested randomised trial (trial = 1) and a delirium risk model with albumin partly missing.
# Simulated data: not evidence about any real drug or patient.
# The simulated registry files (not published) are read from a folder two levels above this script.

set.seed(202610)
# number of imputations m: at least the percentage of incomplete rows (albumin is missing in about 30 percent)
M_IMP <- 40

# ---- packages: check, install into the user library when missing, report ----
need <- c("WeightIt", "cobalt", "survival", "sandwich", "lmtest", "marginaleffects", "mice")
for (p in need) {
  if (!requireNamespace(p, quietly = TRUE)) {
    install.packages(p, repos = "https://cloud.r-project.org", quiet = TRUE)
  }
  cat(sprintf("VERIFY %s %s %s\n", p,
              if (requireNamespace(p, quietly = TRUE)) "available" else "missing",
              if (requireNamespace(p, quietly = TRUE)) as.character(packageVersion(p)) else ""))
}
suppressPackageStartupMessages({
  library(WeightIt); library(cobalt); library(survival); library(sandwich)
  library(lmtest); library(marginaleffects); library(mice)
})

# ---- helpers: print each result on its own CANON line, rounded to 4 decimals ----
canon <- function(key, x) cat(sprintf("CANON w1.%s %.4f\n", key, x))
canon_n <- function(key, x) cat(sprintf("CANON w1.%s %d\n", key, as.integer(x)))
canon_ci <- function(key, est, lo, hi) {
  canon(key, est); canon(paste0(key, ".lo"), lo); canon(paste0(key, ".hi"), hi)
}
z <- qnorm(0.975)
# weighted linear model with robust (HC1) standard errors and t-based CI, as Stata's regress [pw], vce(robust)
# returns estimate, lower, upper, robust SE
wls_rd <- function(f, data, w) {
  data$wt_ <- w
  fit <- lm(f, data = data, weights = wt_)
  V <- vcovHC(fit, type = "HC1")
  ci <- coefci(fit, vcov. = V)
  c(coef(fit)[2], ci[2, ], sqrt(V[2, 2]))
}
# weighted log-link Poisson model for a risk ratio, robust SE treating the weights as known,
# as Stata's glm [pw], family(poisson) link(log) vce(robust) (sandwich times n/(n-1)); returns RR, lower, upper
wpois_rr <- function(f, data, w) {
  data$wt_ <- w
  fit <- glm(f, family = quasipoisson(link = "log"), data = data, weights = wt_)
  n <- nobs(fit)
  b <- coef(fit)[[2]]; se <- sqrt(sandwich(fit)[2, 2] * n / (n - 1))
  c(exp(b), exp(b - z * se), exp(b + z * se))
}
ess <- function(w) sum(w)^2 / sum(w^2)
# simulated truth: the values the simulation was built to produce, read from its settings file (not
# published) and printed as w1.truth.<name>
tj <- jsonlite::read_json(file.path("..", "..", "datasets", "W1", "truth.json"))
truth_echo <- function(nm) canon(paste0("truth.", nm), tj$canon_echo[[nm]])

# ---- data: the simulated registry file and the simulated validation file (not published) ----
d <- read.csv(file.path("..", "..", "datasets", "W1", "W1.csv"))
v <- read.csv(file.path("..", "..", "datasets", "W1", "W1_validation.csv"))
d$age_c <- d$age - 80; d$age_c2 <- d$age_c^2; d$fd <- d$frail * d$dementia
v$age_c <- v$age - 80
d0 <- d[d$trial == 0, ]                     # observational part: treatment chosen by clinicians
canon_n("n", nrow(d)); canon_n("n.trial", sum(d$trial)); canon_n("n.obs", nrow(d0))
canon_n("n.albumin_missing", sum(is.na(d$albumin))); canon_n("n.validation", nrow(v))
# share of registry rows with albumin missing, and the number of imputations m used below
canon("n.albumin_missing_frac", mean(is.na(d$albumin))); canon_n("e1.mi_m", M_IMP)
canon_n("n.obs.treated", sum(d0$surg24)); canon_n("n.obs.died", sum(d0$died)); canon_n("n.obs.lost", sum(d0$lost))

# ---- crude comparison of delirium (observational part) ----
r <- wls_rd(delirium ~ surg24, d0, rep(1, nrow(d0)))
canon_ci("crude.del.rd", r[1], r[2], r[3])

# ---- propensity score models: main effects only versus the true form ----
f_mis <- surg24 ~ age_c + female + frail + dementia + asa3 + anticoag
f_cor <- surg24 ~ age_c + age_c2 + female + frail + dementia + fd + asa3 + anticoag
ps_mis <- fitted(glm(f_mis, family = binomial, data = d0))
ps_cor <- fitted(glm(f_cor, family = binomial, data = d0))
a <- d0$surg24
w_mis <- ifelse(a == 1, 1 / ps_mis, 1 / (1 - ps_mis))
w_ate <- ifelse(a == 1, 1 / ps_cor, 1 / (1 - ps_cor))
w_att <- ifelse(a == 1, 1, ps_cor / (1 - ps_cor))
w_ato <- ifelse(a == 1, 1 - ps_cor, ps_cor)
canon("ps.min", min(ps_cor)); canon("ps.max", max(ps_cor)); canon_n("ps.n_below_005", sum(ps_cor < 0.05))

# ---- standardised mean differences: weighted means, unweighted pooled SD in the denominator ----
vars <- c("age_c", "age_c2", "female", "frail", "dementia", "fd", "asa3", "anticoag")
smd <- function(x, w) {
  m1 <- weighted.mean(x[a == 1], w[a == 1]); m0 <- weighted.mean(x[a == 0], w[a == 0])
  (m1 - m0) / sqrt((var(x[a == 1]) + var(x[a == 0])) / 2)
}
for (lab in c("raw", "mis", "cor")) {
  w <- switch(lab, raw = rep(1, nrow(d0)), mis = w_mis, cor = w_ate)
  s <- sapply(vars, function(x) smd(d0[[x]], w))
  for (x in vars) canon(sprintf("smd.%s.%s", lab, x), s[[x]])
  canon(sprintf("smd.%s.maxabs", lab), max(abs(s)))
}
# the same balance tables as cobalt reports them (display only)
W_mis <- weightit(f_mis, data = d0, method = "glm", estimand = "ATE")
W_cor <- weightit(f_cor, data = d0, method = "glm", estimand = "ATE")
print(bal.tab(W_mis, data = d0, addl = ~ age_c2 + fd, un = TRUE, s.d.denom = "pooled"))
print(bal.tab(W_cor, un = TRUE, s.d.denom = "pooled"))

# ---- weights: raw, stabilised, truncated at the 1st and 99th percentiles ----
pa <- mean(a)
w_sw <- w_ate * ifelse(a == 1, pa, 1 - pa)                 # stabilised: one constant per arm
q <- quantile(w_ate, c(0.01, 0.99), type = 2)              # type 2 matches Stata's _pctile
w_tr <- pmin(pmax(w_ate, q[1]), q[2])
canon("wt.trunc.p01", q[1]); canon("wt.trunc.p99", q[2])
for (lab in c("raw", "sw", "trunc")) {
  w <- switch(lab, raw = w_ate, sw = w_sw, trunc = w_tr)
  canon(sprintf("wt.%s.max", lab), max(w)); canon(sprintf("wt.%s.ess", lab), ess(w))
  canon(sprintf("wt.%s.ess_treated", lab), ess(w[a == 1])); canon(sprintf("wt.%s.ess_control", lab), ess(w[a == 0]))
}

# ---- IPTW effects on delirium ----
# ATE and ATT: weighted outcome model with M-estimation SE that accounts for the estimated score
fit_ate <- lm_weightit(delirium ~ surg24, data = d0, weightit = W_cor)
fit_mis <- lm_weightit(delirium ~ surg24, data = d0, weightit = W_mis)
W_att <- weightit(f_cor, data = d0, method = "glm", estimand = "ATT")
fit_att <- lm_weightit(delirium ~ surg24, data = d0, weightit = W_att)
for (nm in c("ate", "att", "ate_mis")) {
  fit <- switch(nm, ate = fit_ate, att = fit_att, ate_mis = fit_mis)
  b <- coef(fit)[["surg24"]]; se <- sqrt(vcov(fit)["surg24", "surg24"])
  canon_ci(paste0("ipw.del.", nm), b, b - z * se, b + z * se)
  # SE of the ATE that accounts for estimating e(X) (M-estimation)
  if (nm == "ate") canon("ipw.del.ate.se", se)
}
# ATE risk ratio from the two weighted risks (log link on the weighted means, same M-estimation SE)
fit_rr <- glm_weightit(delirium ~ surg24, data = d0, weightit = W_cor, family = quasipoisson(link = "log"))
b <- coef(fit_rr)[["surg24"]]; se <- sqrt(vcov(fit_rr)["surg24", "surg24"])
canon("ipw.del.ate.risk1", weighted.mean(d0$delirium[a == 1], w_ate[a == 1]))
canon("ipw.del.ate.risk0", weighted.mean(d0$delirium[a == 0], w_ate[a == 0]))
canon_ci("ipw.del.ate.rr", exp(b), exp(b - z * se), exp(b + z * se))
# misspecified score: the same risk ratio and M-estimation SE under the main-effects-only propensity model
fit_rr_mis <- glm_weightit(delirium ~ surg24, data = d0, weightit = W_mis, family = quasipoisson(link = "log"))
b <- coef(fit_rr_mis)[["surg24"]]; se <- sqrt(vcov(fit_rr_mis)["surg24", "surg24"])
canon_ci("ipw.del.ate_mis.rr", exp(b), exp(b - z * se), exp(b + z * se))
# ATO, stabilised, truncated and trimmed: weighted regression, robust SE treating weights as known
# raw ATE weights through the same weights-known estimator, so raw and stabilised compare like for like
r <- wls_rd(delirium ~ surg24, d0, w_ate); canon_ci("ipw.del.ate_rawreg", r[1], r[2], r[3])
canon("ipw.del.ate_rawreg.se", r[4])                         # weights-known robust SE, raw weights
r <- wls_rd(delirium ~ surg24, d0, w_ato); canon_ci("ipw.del.ato", r[1], r[2], r[3])
r <- wls_rd(delirium ~ surg24, d0, w_sw); canon_ci("ipw.del.ate_sw", r[1], r[2], r[3])
canon("ipw.del.ate_sw.se", r[4])                             # weights-known robust SE, stabilised weights
r <- wls_rd(delirium ~ surg24, d0, w_tr); canon_ci("ipw.del.ate_trunc", r[1], r[2], r[3])
keep <- ps_cor >= 0.1 & ps_cor <= 0.9                       # trimming changes the population
canon_n("n.trim", sum(keep))
r <- wls_rd(delirium ~ surg24, d0[keep, ], w_ate[keep]); canon_ci("ipw.del.ate_trim", r[1], r[2], r[3])
# risk ratios beside each weights-known risk difference (log-link Poisson, robust SE, weights known)
r <- wpois_rr(delirium ~ surg24, d0, w_ate); canon_ci("ipw.del.ate_rawreg.rr", r[1], r[2], r[3])
r <- wpois_rr(delirium ~ surg24, d0, w_sw); canon_ci("ipw.del.ate_sw.rr", r[1], r[2], r[3])
r <- wpois_rr(delirium ~ surg24, d0, w_tr); canon_ci("ipw.del.ate_trunc.rr", r[1], r[2], r[3])
r <- wpois_rr(delirium ~ surg24, d0[keep, ], w_ate[keep]); canon_ci("ipw.del.ate_trim.rr", r[1], r[2], r[3])
r <- wpois_rr(delirium ~ surg24, d0, w_ato); canon_ci("ipw.del.ato.rr", r[1], r[2], r[3])
# the true values these delirium estimates are compared with (simulated truth)
for (nm in c("del_rd_ate_trial0", "del_rr_ate_trial0", "del_rd_att_trial0", "del_rr_att_trial0",
             "del_rd_ato_trial0", "del_rr_ato_trial0", "del_rd_trimmed_on_true_ps_trial0",
             "del_rr_trimmed_on_true_ps_trial0", "del_assoc_trial0_rd", "del_assoc_trial0_rr")) truth_echo(nm)

# ---- 1-year death: crude Kaplan-Meier risk and crude Cox ----
km <- summary(survfit(Surv(fu_months, died) ~ surg24, data = d0), times = 12)
risk <- 1 - km$surv; se_km <- km$std.err                    # strata order: surg24 = 0, then 1
canon("km.risk0", risk[1]); canon("km.risk1", risk[2])
rd <- risk[2] - risk[1]; se <- sqrt(sum(se_km^2))
canon_ci("crude.death.rd", rd, rd - z * se, rd + z * se)
cx <- coxph(Surv(fu_months, died) ~ surg24, data = d0, ties = "breslow")
b <- coef(cx)[[1]]; se <- sqrt(vcov(cx)[1, 1])
canon_ci("crude.death.hr", exp(b), exp(b - z * se), exp(b + z * se))

# ---- inverse probability of censoring weights: exponential model for loss to follow-up ----
cm <- survreg(Surv(fu_months, lost) ~ dementia + age_c + surg24, data = d0, dist = "exponential")
rate_c <- exp(-predict(cm, type = "lp"))                    # survreg is on the log-time scale
G <- exp(-rate_c * d0$fu_months)                            # P(still followed at own end time | X)
obs <- d0$lost == 0                                         # vital status at 12 months is known
ipcw <- 1 / G
canon("ipcw.max", max(ipcw[obs]))
r <- wls_rd(died ~ surg24, d0[obs, ], rep(1, sum(obs))); canon_ci("naive.death.rd", r[1], r[2], r[3])
r <- wls_rd(died ~ surg24, d0[obs, ], ipcw[obs]); canon_ci("ipcw.death.rd", r[1], r[2], r[3])
# truth for the IPCW-only contrast: censoring-free ASSOCIATIONAL risks by arm in trial = 0 (no confounding control)
for (nm in c("death_assoc_trial0_risk1", "death_assoc_trial0_risk0", "death_assoc_trial0_rd",
             "death_assoc_trial0_rr")) truth_echo(nm)
# truth for the IPTW x IPCW contrasts below: the causal 1-year risk difference and risk ratio (ATE, trial = 0)
for (nm in c("death_rd_ate_trial0", "death_rr_ate_trial0")) truth_echo(nm)
# ---- IPTW times IPCW: the 1-year risk difference for ATE, ATT and ATO ----
for (nm in c("ate", "att", "ato")) {
  w <- switch(nm, ate = w_ate, att = w_att, ato = w_ato) * ipcw
  r <- wls_rd(died ~ surg24, d0[obs, ], w[obs]); canon_ci(paste0("ipw.death.", nm), r[1], r[2], r[3])
}
# ---- IPTW Cox model (ATE weights, robust SE) ----
cw <- coxph(Surv(fu_months, died) ~ surg24, data = d0, weights = w_ate, robust = TRUE, ties = "breslow")
b <- coef(cw)[[1]]; se <- sqrt(vcov(cw)[1, 1])
canon_ci("ipw.death.hr", exp(b), exp(b - z * se), exp(b + z * se))

# ---- link and family: saturated versus adjusted fits (all registry rows) ----
rr_glm <- function(f, fam, data, robust = FALSE, start = NULL) {
  fit <- glm(f, family = fam, data = data, start = start)
  V <- if (robust) vcovHC(fit, type = "HC0") else vcov(fit)
  b <- coef(fit)[["surg24"]]; se <- sqrt(V["surg24", "surg24"])
  list(fit = fit, est = exp(b), lo = exp(b - z * se), hi = exp(b + z * se), se = se,
       se_naive = sqrt(vcov(fit)["surg24", "surg24"]))
}
# saturated fits (surg24 only): log-binomial, modified Poisson, and Gaussian family with a log link
cb <- rr_glm(delirium ~ surg24, binomial(link = "log"), d)
canon_ci("c2.rr_crude_logbin", cb$est, cb$lo, cb$hi)          # model-based SE (identical bread in both languages here)
canon("c2.se_crude_logbin", cb$se)                            # log-RR SE, log-binomial, model-based
cp <- rr_glm(delirium ~ surg24, poisson(link = "log"), d, TRUE)
canon_ci("c2.rr_crude_poisson", cp$est, cp$lo, cp$hi)         # modified Poisson, robust CI
canon("c2.se_crude_poisson_naive", cp$se_naive)               # Poisson model-based log-RR SE (too large for a binary outcome)
canon("c2.se_crude_poisson_robust", cp$se)                    # sandwich log-RR SE
canon("c2.rr_crude_poisson_naive.lo", exp(log(cp$est) - z * cp$se_naive))   # naive Poisson 95% CI, lower
canon("c2.rr_crude_poisson_naive.hi", exp(log(cp$est) + z * cp$se_naive))   # naive Poisson 95% CI, upper
cg <- rr_glm(delirium ~ surg24, gaussian(link = "log"), d, TRUE, start = coef(cp$fit))
canon_ci("c2.rr_crude_gaussian", cg$est, cg$lo, cg$hi)        # Gaussian log link, robust CI (saturated: same bread both languages)
# adjusted fits: the log-binomial needs starting values (taken from the Poisson fit) to converge
p2 <- rr_glm(delirium ~ surg24 + age_c + female, poisson(link = "log"), d, TRUE)
l2 <- rr_glm(delirium ~ surg24 + age_c + female, binomial(link = "log"), d, start = coef(p2$fit))
canon("c2.rr_adj2_logbin", l2$est)
# R's glm uses the expected-information SE for the non-canonical log link (Stata's ML glm: observed)
canon("c2.rr_adj2_logbin.lo.r", l2$lo); canon("c2.rr_adj2_logbin.hi.r", l2$hi)
canon("c2.se_adj2_logbin.r", l2$se)                           # log-RR SE, adjusted log-binomial (expected information)
canon("c2.rr_adj2_poisson", p2$est)
canon("c2.rr_adj2_poisson.lo", p2$lo); canon("c2.rr_adj2_poisson.hi", p2$hi)   # modified Poisson, robust CI
canon("c2.se_adj2_poisson_naive", p2$se_naive)                # Poisson model-based log-RR SE, adjusted
canon("c2.se_adj2_poisson_robust", p2$se)                     # sandwich log-RR SE, adjusted
g2 <- rr_glm(delirium ~ surg24 + age_c + female, gaussian(link = "log"), d, TRUE, start = coef(p2$fit))
canon("c2.rr_adj2_gaussian", g2$est)                          # Gaussian log link, adjusted for age and sex
canon("c2.rr_adj2_gaussian.lo.r", g2$lo); canon("c2.rr_adj2_gaussian.hi.r", g2$hi)   # expected-information bread
# adjusted for age and frailty: the Poisson fit gives fitted risks above 1, so the log-binomial cannot start
paf <- rr_glm(delirium ~ surg24 + age_c + frail, poisson(link = "log"), d, TRUE)
canon("c2.rr_adjaf_poisson", paf$est)
canon("c2.rr_adjaf_poisson.lo", paf$lo); canon("c2.rr_adjaf_poisson.hi", paf$hi)
canon("c2.se_adjaf_poisson_naive", paf$se_naive)              # Poisson model-based log-RR SE, age and frailty
canon("c2.se_adjaf_poisson_robust", paf$se)                   # sandwich log-RR SE, age and frailty
canon("c2.adjaf_poisson_max_fitted", max(fitted(paf$fit)))    # largest fitted risk: above 1 means the log-binomial wall
lbaf <- tryCatch(glm(delirium ~ surg24 + age_c + frail, family = binomial(link = "log"), data = d, start = coef(paf$fit)),
                 error = function(e) { cat("NOTE log-binomial (age, frailty) stopped:", conditionMessage(e), "\n"); NULL })
canon_n("c2.adjaf_logbin_converged.r", !is.null(lbaf) && isTRUE(lbaf$converged))   # 1 = converged

# ---- the risk-ratio ladder for a common outcome (all registry rows, same covariates) ----
f_out <- delirium ~ surg24 + age_c + female + frail + dementia + asa3 + anticoag
# rung 1: log-binomial; it stops when a covariate pattern would need a risk above 1
lb <- tryCatch(glm(f_out, family = binomial(link = "log"), data = d),
               error = function(e) { cat("NOTE log-binomial stopped:", conditionMessage(e), "\n"); NULL })
canon_n("c3.logbin_converged", !is.null(lb) && isTRUE(lb$converged))
# rung 2: modified Poisson with robust SE
mp <- rr_glm(f_out, poisson(link = "log"), d, TRUE)
canon_ci("c3.rr_poisson", mp$est, mp$lo, mp$hi)
canon("c3.poisson_max_fitted", max(fitted(mp$fit)))
canon_n("c3.poisson_n_fitted_above1", sum(fitted(mp$fit) > 1))
# rung 3: Gaussian family with a LOG link and robust SE (starts from the Poisson fit)
gl <- rr_glm(f_out, gaussian(link = "log"), d, TRUE, start = coef(mp$fit))
# R's sandwich uses the expected-information bread, so this robust CI differs from Stata's ML glm
canon("c3.rr_gaussian", gl$est); canon("c3.rr_gaussian.lo.r", gl$lo); canon("c3.rr_gaussian.hi.r", gl$hi)
# rung 4: logistic model, then marginal standardisation for the risk ratio and risk difference
lg <- glm(f_out, family = binomial, data = d)
b <- coef(lg)[["surg24"]]; se <- sqrt(vcov(lg)["surg24", "surg24"])
canon_ci("c5.or_cond", exp(b), exp(b - z * se), exp(b + z * se))
std <- function(cmp) avg_comparisons(lg, variables = list(surg24 = c(0, 1)), comparison = cmp)
rr <- std("lnratioavg"); canon_ci("c3.rr_std", exp(rr$estimate), exp(rr$conf.low), exp(rr$conf.high))
rd <- std("differenceavg"); canon_ci("c3.rd_std", rd$estimate, rd$conf.low, rd$conf.high)
# the same model averaged over the cohort gives the marginal odds ratio (non-collapsibility)
om <- std("lnoravg"); canon_ci("c5.or_marg", exp(om$estimate), exp(om$conf.low), exp(om$conf.high))
# the realised nested trial (n = 950): one noisy draw. Chance covariate imbalance in a trial this small can
# outweigh non-collapsibility, so these keys are named "realised"; the large-sample pair is printed further down
tr <- d[d$trial == 1, ]
for (nm in c("crude", "adj")) {
  f <- if (nm == "crude") delirium ~ surg24 else f_out
  ft <- glm(f, family = binomial, data = tr)
  b <- coef(ft)[["surg24"]]; se <- sqrt(vcov(ft)["surg24", "surg24"])
  canon_ci(paste0("c5.trial_realised.or_", nm), exp(b), exp(b - z * se), exp(b + z * se))
}
# trial-standardised marginal OR: the adjusted trial model averaged over the trial participants
ltr <- glm(f_out, family = binomial, data = tr)
om_t <- avg_comparisons(ltr, variables = list(surg24 = c(0, 1)), comparison = "lnoravg")
canon_ci("c5.trial_realised.or_std", exp(om_t$estimate), exp(om_t$conf.low), exp(om_t$conf.high))
# large-sample simulated truth: conditional OR within frailty strata versus marginal ORs
# in trial participants, and the large-sample value of the six-covariate adjusted model in the trial
for (nm in c("c5_or_cond_nonfrail", "c5_or_cond_frail", "c5_or_marg_trial_nonfrail", "c5_or_marg_trial_frail",
             "c5_or_marg_trial", "c5_or_adj_pseudo_trial", "del_rr_whole_population", "del_rd_whole_population"))
  truth_echo(nm)

# ---- prediction model for delirium with albumin missing: complete case versus MI without and with Y ----
f_pred <- delirium ~ age_c + female + frail + dementia + asa3 + anticoag + albumin
auc <- function(y, s) {
  rk <- rank(s); n1 <- as.numeric(sum(y == 1)); n0 <- as.numeric(sum(y == 0))
  (sum(rk[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
val_auc <- function(b) auc(v$delirium, as.vector(model.matrix(f_pred, v) %*% b))
# calibration of a linear predictor lp on the outcome y: the slope is the coefficient of lp in a logistic
# regression of y on lp; calibration-in-the-large (CITL) is the intercept when lp enters as an offset
# (slope fixed at 1). Ideal values: slope 1, CITL 0.
cal_fit <- function(y, lp) {
  fs <- glm(y ~ lp, family = binomial)
  fc <- glm(y ~ 1, offset = lp, family = binomial)
  c(slope = coef(fs)[[2]], slope_se = sqrt(vcov(fs)[2, 2]), citl = coef(fc)[[1]], citl_se = sqrt(vcov(fc)[1, 1]))
}
# slope and CITL with normal 95% CIs: (slope, lo, hi, CITL, lo, hi)
cal_normal <- function(y, lp) {
  r <- cal_fit(y, lp)
  c(r[["slope"]] + c(0, -z, z) * r[["slope_se"]], r[["citl"]] + c(0, -z, z) * r[["citl_se"]])
}
# Rubin's rules over the m rows of (slope, SE, CITL, SE): pooled estimate, total variance W + (1 + 1/m) B,
# and a t-based 95% CI with the large-sample Rubin df; returns (slope, lo, hi, CITL, lo, hi)
rubin_cal <- function(X) {
  M <- nrow(X); out <- numeric(0)
  for (j in c(1, 3)) {
    q <- mean(X[, j]); w <- mean(X[, j + 1]^2); b <- var(X[, j])
    tv <- w + (1 + 1 / M) * b
    df <- (M - 1) * (1 + w / ((1 + 1 / M) * b))^2
    cq <- qt(0.975, df)
    out <- c(out, q, q - cq * sqrt(tv), q + cq * sqrt(tv))
  }
  out
}
# print slope and CITL as <key>.slope, <key>.citl (and .lo, .hi) and keep a row for the table
cal_tab <- NULL
cal_print <- function(key, row, sfx = "") {
  nm <- c("slope", "slope.lo", "slope.hi", "citl", "citl.lo", "citl.hi")
  for (k in seq_along(nm)) canon(paste0(key, ".", nm[k], sfx), row[k])
  cal_tab <<- rbind(cal_tab, setNames(row, c("slope", "slope_lo", "slope_hi", "citl", "citl_lo", "citl_hi")))
}
cc <- glm(f_pred, family = binomial, data = d)              # glm drops rows with albumin missing
canon_n("e1.cc.n", nobs(cc))
b <- coef(cc)[["albumin"]]; se <- sqrt(vcov(cc)["albumin", "albumin"])
canon_ci("e1.cc.b_albumin", b, b - z * se, b + z * se)
print(round(coef(summary(cc)), 4))                          # complete-case coefficient table
auc_val <- c(cc = val_auc(coef(cc)))
auc_dev <- c(cc = auc(cc$y, cc$linear.predictors))          # apparent AUROC on the complete-case development rows
canon("e1.cc.auroc", auc_val[["cc"]])
canon("e1.cc.dev_auroc", auc_dev[["cc"]])
xvars <- c("age_c", "female", "frail", "dementia", "asa3", "anticoag")
mi_fit <- function(with_y) {
  cols <- c(xvars, "albumin", "delirium")
  dd <- d[, cols]
  pm <- make.predictorMatrix(dd); pm[, ] <- 0
  pm["albumin", xvars] <- 1
  if (with_y) pm["albumin", "delirium"] <- 1
  imp <- mice(dd, m = M_IMP, method = "pmm", donors = 10, predictorMatrix = pm, seed = 202610, printFlag = FALSE)
  s <- summary(pool(with(imp, glm(delirium ~ age_c + female + frail + dementia + asa3 + anticoag + albumin,
                                  family = binomial))), conf.int = TRUE)
  list(s = s, imp = imp)
}
lp_val <- list(cc = as.vector(model.matrix(f_pred, v) %*% coef(cc)))
for (nm in c("mi_noy", "mi_y")) {
  res <- mi_fit(nm == "mi_y")
  s <- res$s
  cat(sprintf("Pooled logistic model, imputation %s the outcome (m = %d)\n", if (nm == "mi_y") "with" else "without", M_IMP))
  print(data.frame(term = s$term, estimate = round(s$estimate, 4), se = round(s$std.error, 4),
                   lo = round(s$`2.5 %`, 4), hi = round(s$`97.5 %`, 4)))
  bb <- setNames(s$estimate, as.character(s$term))
  i <- which(s$term == "albumin")
  canon(sprintf("e1.%s.b_albumin.r", nm), s$estimate[i])
  canon(sprintf("e1.%s.b_albumin.lo.r", nm), s$`2.5 %`[i]); canon(sprintf("e1.%s.b_albumin.hi.r", nm), s$`97.5 %`[i])
  bv <- bb[colnames(model.matrix(f_pred, v))]
  auc_val[[nm]] <- val_auc(bv)
  canon(sprintf("e1.%s.auroc.r", nm), auc_val[[nm]])
  lp_val[[nm]] <- as.vector(model.matrix(f_pred, v) %*% bv)
  # apparent development AUROC of the pooled model, averaged over the m completed development sets
  dev_auc <- mean(sapply(seq_len(M_IMP), function(m) {
    dm <- complete(res$imp, m)
    auc(dm$delirium, as.vector(model.matrix(f_pred, dm) %*% bv))
  }))
  auc_dev[[nm]] <- dev_auc
  canon(sprintf("e1.%s.dev_auroc.r", nm), dev_auc)
}
cat("AUROC of each fitted model: development rows (apparent) and validation set, simulated data\n")
print(round(cbind(development = auc_dev, validation = auc_val), 4))
# calibration slope and calibration-in-the-large of each fitted model on the validation file, normal 95% CIs
cal_print("e1.cc.cal", cal_normal(v$delirium, lp_val$cc))
cal_print("e1.mi_noy.cal", cal_normal(v$delirium, lp_val$mi_noy), ".r")
cal_print("e1.mi_y.cal", cal_normal(v$delirium, lp_val$mi_y), ".r")

# validation with missing albumin: the validation file has albumin fully observed, so a reproducible MAR mask
# is applied in memory (the file is not changed). The mask uses the missingness model the registry was
# simulated with and a deterministic uniform u = frac(id x 0.6180339887498949), identical in Stata and R.
pmiss <- tj$parameters$albumin_missing
u <- (v$id * 0.6180339887498949) %% 1
p_m <- plogis(pmiss$intercept + pmiss$dementia * v$dementia + pmiss$asa3 * v$asa3 + pmiss$age_c * v$age_c +
                pmiss$frail * v$frail)
vm <- v
vm$albumin[u < p_m] <- NA
canon_n("e1.val.n_masked", sum(is.na(vm$albumin)))           # validation rows whose albumin is masked
b_cc <- coef(cc)                                              # one fixed model is scored: the complete-case fit
lp_of <- function(dd) as.vector(model.matrix(f_pred, dd) %*% b_cc)
auc_mask <- c(full = auc(v$delirium, lp_of(v)))              # albumin fully observed (same as e1.cc.auroc)
canon("e1.val.auroc_full", auc_mask[["full"]])
cal_print("e1.val.cal_full", cal_normal(v$delirium, lp_of(v)))
# deployment-style single regression imputation fitted on the development data, no outcome anywhere
ri <- lm(albumin ~ age_c + female + frail + dementia + asa3 + anticoag, data = d)
vr <- vm
vr$albumin[is.na(vr$albumin)] <- predict(ri, newdata = vr[is.na(vr$albumin), ])
auc_mask[["regimp"]] <- auc(vr$delirium, lp_of(vr))          # Y-free regression imputation from development data
canon("e1.val.auroc_regimp", auc_mask[["regimp"]])
# single imputation: these calibration CIs ignore the uncertainty of the filled values
cal_print("e1.val.cal_regimp", cal_normal(vr$delirium, lp_of(vr)))
# multiple imputation inside the validation sample, with versus without the validation outcomes:
# AUROC averaged over the m imputations, slope and CITL pooled with Rubin's rules
val_mi <- function(with_y) {
  dd <- vm[, c(xvars, "albumin", "delirium")]
  pm <- make.predictorMatrix(dd); pm[, ] <- 0
  pm["albumin", xvars] <- 1
  if (with_y) pm["albumin", "delirium"] <- 1
  imp <- mice(dd, m = M_IMP, method = "pmm", donors = 10, predictorMatrix = pm, seed = 202611, printFlag = FALSE)
  per <- t(sapply(seq_len(M_IMP), function(m) {
    dm <- complete(imp, m); lp <- lp_of(dm)
    c(auc = auc(dm$delirium, lp), cal_fit(dm$delirium, lp))
  }))
  list(auc = mean(per[, "auc"]), cal = rubin_cal(per[, c("slope", "slope_se", "citl", "citl_se")]))
}
for (wy in c("with", "without")) {
  r <- val_mi(wy == "with")
  auc_mask[[if (wy == "with") "mi_y" else "mi_noy"]] <- r$auc
  canon(sprintf("e1.val.auroc_mi_%s_y.r", wy), r$auc)        # validation albumin imputed with / without the outcomes
  cal_print(sprintf("e1.val.cal_mi_%s_y", wy), r$cal, ".r")
}
# calibration table: the three fitted models on the validation file, then the fixed complete-case model
# after the masked validation albumin was filled four ways (slope 1 and CITL 0 mean perfect calibration);
# val_mi_y and val_mi_noy: multiple imputation in the validation file with and without the validation outcomes
rownames(cal_tab) <- c("cc", "mi_noy", "mi_y", "val_full", "val_regimp", "val_mi_y", "val_mi_noy")
cat("AUROC of the complete-case model after the masked validation albumin was filled four ways, simulated data\n")
print(round(auc_mask, 4))
cat("Calibration on the validation set, simulated data\n")
print(round(cal_tab, 4))
# the reference value for the albumin coefficient: this working model fitted to a very large simulated
# population with every albumin value recorded, and that reference model's AUROC on the validation file
canon("truth.pseudo_true_albumin", tj$prediction_model_delirium$pseudo_true_coefficients$albumin)
canon("truth.auroc_pseudo_true_on_validation", tj$prediction_model_delirium$auroc_pseudo_true_on_validation)

# ---- nested trial: sampling-score weights to move the trial result to a target population ----
sm <- glm(trial ~ age_c + female + frail + dementia + asa3 + anticoag, family = binomial, data = d)
s_hat <- fitted(sm)[d$trial == 1]
w_ipsw <- 1 / s_hat                                         # target: the whole registry
w_iosw <- (1 - s_hat) / s_hat                               # target: the non-participants (inverse odds)
r <- wls_rd(delirium ~ surg24, tr, rep(1, nrow(tr))); canon_ci("a2.trial.rd", r[1], r[2], r[3])
se_trial <- r[4]; canon("a2.trial.rd.se", se_trial)          # robust SE of the unweighted trial difference
r <- wls_rd(delirium ~ surg24, tr, w_ipsw); canon_ci("a2.ipsw.rd", r[1], r[2], r[3])
canon("a2.ipsw.rd.se", r[4])                                 # robust SE after weighting to the whole registry
canon("a2.ipsw.se_ratio", r[4] / se_trial)                   # variance cost: IPSW SE over the trial SE
r <- wls_rd(delirium ~ surg24, tr, w_iosw); canon_ci("a2.iosw.rd", r[1], r[2], r[3])
canon("a2.iosw.rd.se", r[4])                                 # robust SE after weighting to the non-participants
canon("a2.iosw.se_ratio", r[4] / se_trial)                   # variance cost: inverse-odds SE over the trial SE
canon("a2.ipsw.max", max(w_ipsw)); canon("a2.ipsw.ess", ess(w_ipsw))
canon("a2.iosw.max", max(w_iosw)); canon("a2.iosw.ess", ess(w_iosw))
canon("a2.ipsw.ess_frac", ess(w_ipsw) / nrow(tr))            # ESS as a share of the 950 trial participants
canon("a2.iosw.ess_frac", ess(w_iosw) / nrow(tr))
# risk ratios beside the transported risk differences (log-link Poisson, robust SE, weights known)
r <- wpois_rr(delirium ~ surg24, tr, rep(1, nrow(tr))); canon_ci("a2.trial.rr", r[1], r[2], r[3])
r <- wpois_rr(delirium ~ surg24, tr, w_ipsw); canon_ci("a2.ipsw.rr", r[1], r[2], r[3])
r <- wpois_rr(delirium ~ surg24, tr, w_iosw); canon_ci("a2.iosw.rr", r[1], r[2], r[3])
# their targets: trial participants, the whole registry population, the non-participants (trial = 0)
for (nm in c("del_rd_trial_participants", "del_rr_trial_participants", "del_rd_obs_part_trial0",
             "del_rr_obs_part_trial0")) truth_echo(nm)
cat("All numbers above are from simulated data (ข้อมูลจำลอง).\n")

# ---- weighting comparisons: crude risks, simulated truth and registry facts (observational part) ----
# crude delirium risk in each arm and the crude risk ratio (log-link Poisson, robust SE)
canon("crude.del.risk1", mean(d0$delirium[a == 1])); canon("crude.del.risk0", mean(d0$delirium[a == 0]))
r <- wpois_rr(delirium ~ surg24, d0, rep(1, nrow(d0))); canon_ci("crude.del.rr", r[1], r[2], r[3])
# simulated truth: delirium risk by the arm actually received, with no confounding control
for (nm in c("del_assoc_trial0_risk1", "del_assoc_trial0_risk0")) truth_echo(nm)
# simulated truth: delirium risk if every patient had early surgery (risk1) or later surgery (risk0)
canon("truth.del_risk1_ate_trial0", tj$true_delirium$ate_trial0$risk1)
canon("truth.del_risk0_ate_trial0", tj$true_delirium$ate_trial0$risk0)
# the TRUE propensity score of each patient (the form the simulation used, with age squared and
# frailty x dementia): its range and the share of patients below 0.05
tp <- tj$parameters$ps
e_true <- plogis(tp$intercept + tp$age_c * d0$age_c + tp$age_c_sq * d0$age_c2 + tp$frail * d0$frail +
                   tp$dementia * d0$dementia + tp$frail_x_dementia * d0$fd + tp$asa3 * d0$asa3 +
                   tp$anticoag * d0$anticoag + tp$female * d0$female)
canon("truth.ps_min_trial0", min(e_true)); canon("truth.ps_max_trial0", max(e_true))
canon("truth.ps_below_005_trial0", mean(e_true < 0.05))
# share operated within 24 hours: the constant that stabilises the early-surgery weights (1 minus it for the rest)
canon("wt.sw.p_treated", pa); canon("wt.sw.p_control", 1 - pa)
# registry facts (all 20,000 rows): age SD, and loss to follow-up with and without dementia
canon("reg.age_sd", sd(d$age))
canon("reg.lost_dementia1", mean(d$lost[d$dementia == 1])); canon("reg.lost_dementia0", mean(d$lost[d$dementia == 0]))
# half-width of each 95% CI in the nested trial: the precision a transported estimate gives up
r <- wls_rd(delirium ~ surg24, tr, rep(1, nrow(tr))); canon("a2.trial.rd.halfwidth", (r[3] - r[2]) / 2)
r <- wls_rd(delirium ~ surg24, tr, w_ipsw); canon("a2.ipsw.rd.halfwidth", (r[3] - r[2]) / 2)
r <- wls_rd(delirium ~ surg24, tr, w_iosw); canon("a2.iosw.rd.halfwidth", (r[3] - r[2]) / 2)
# arm size of the later-surgery group, and the ATT weighted risks (early-surgery arm as observed, later-surgery
# arm reweighted to look like it) beside their simulated truth
canon_n("n.obs.control", sum(a == 0))
canon("ipw.del.att.risk1", mean(d0$delirium[a == 1]))
canon("ipw.del.att.risk0", weighted.mean(d0$delirium[a == 0], w_att[a == 0]))
canon("truth.del_risk1_att_trial0", tj$true_delirium$att_trial0$risk1)
canon("truth.del_risk0_att_trial0", tj$true_delirium$att_trial0$risk0)
# variance ratio, treated over control, from weighted variances (weights scaled to sum to the arm size,
# divisor n - 1): before weighting, after the main-effects model, after the revised model
wvar <- function(x, w) {
  w <- w * length(w) / sum(w)
  m <- sum(w * x) / sum(w)
  sum(w * (x - m)^2) / (length(x) - 1)
}
vr <- function(x, w) wvar(x[a == 1], w[a == 1]) / wvar(x[a == 0], w[a == 0])
for (lab in c("raw", "mis", "cor")) {
  w <- switch(lab, raw = rep(1, nrow(d0)), mis = w_mis, cor = w_ate)
  for (x in vars) canon(sprintf("vr.%s.%s", lab, x), vr(d0[[x]], w))
}
# overlap: quantiles of the estimated propensity score (revised model) in each arm
qs <- c(min = 0, p1 = 0.01, p5 = 0.05, p50 = 0.5, p95 = 0.95, p99 = 0.99, max = 1)
ps_q <- t(sapply(c(early = 1, later = 0), function(g) quantile(ps_cor[a == g], qs, type = 2)))
colnames(ps_q) <- names(qs)
for (arm in rownames(ps_q)) for (k in names(qs)) canon(sprintf("ps.%s.%s", arm, k), ps_q[arm, k])
# balance table: standardised mean differences (weighted means over the unweighted pooled SD) and variance
# ratios, before weighting, after the main-effects model and after the revised model (age squared, frail x dementia)
one <- rep(1, nrow(d0))
bal <- t(sapply(vars, function(x) c(
  smd_before = smd(d0[[x]], one), smd_main = smd(d0[[x]], w_mis), smd_revised = smd(d0[[x]], w_ate),
  vr_before = vr(d0[[x]], one), vr_main = vr(d0[[x]], w_mis), vr_revised = vr(d0[[x]], w_ate))))
cat("Covariate balance in the observational part, simulated data\n")
print(round(bal, 3))

# ---- summary tables of the weighting results (simulated data) ----
# target_rd is the value each estimate aims at: the associational difference for the crude comparison and
# the naive or IPCW-only death contrasts, the causal effect for the weighted ones (simulated truth)
te <- tj$canon_echo
mest <- function(fit) {                      # estimate and 95% CI, M-estimation SE (accounts for e(X))
  b <- coef(fit)[["surg24"]]; se <- sqrt(vcov(fit)["surg24", "surg24"])
  c(b, b - z * se, b + z * se)
}
# overlap: the estimated propensity score (revised model) by arm
cat("Estimated propensity score by arm, observational part, simulated data\n")
print(round(ps_q, 4))
# delirium: crude, IPTW ATE and IPTW ATT
t_iptw <- rbind(
  crude = c(mean(d0$delirium[a == 1]), mean(d0$delirium[a == 0]),
            wls_rd(delirium ~ surg24, d0, one)[1:3], te$del_assoc_trial0_rd),
  iptw_ate = c(weighted.mean(d0$delirium[a == 1], w_ate[a == 1]), weighted.mean(d0$delirium[a == 0], w_ate[a == 0]),
               mest(fit_ate), te$del_rd_ate_trial0),
  iptw_att = c(mean(d0$delirium[a == 1]), weighted.mean(d0$delirium[a == 0], w_att[a == 0]),
               mest(fit_att), te$del_rd_att_trial0))
colnames(t_iptw) <- c("risk_early", "risk_later", "rd", "rd_lo", "rd_hi", "target_rd")
cat("Delirium, early versus later surgery, observational part, simulated data\n")
print(round(t_iptw, 4))
# extreme weights: the weights themselves, then the risk difference under each choice (robust SE, weights known)
t_w <- t(sapply(list(raw = w_ate, stabilised = w_sw, truncated = w_tr), function(w)
  c(max = max(w), ess_early = ess(w[a == 1]), ess_later = ess(w[a == 0]))))
cat("ATE weights: largest weight and effective sample size per arm, simulated data\n")
print(round(t_w, 4))
t_ext <- rbind(
  raw = c(wls_rd(delirium ~ surg24, d0, w_ate)[1:3], te$del_rd_ate_trial0),
  stabilised = c(wls_rd(delirium ~ surg24, d0, w_sw)[1:3], te$del_rd_ate_trial0),
  truncated_p1_p99 = c(wls_rd(delirium ~ surg24, d0, w_tr)[1:3], te$del_rd_ate_trial0),
  trimmed_0.1_0.9 = c(wls_rd(delirium ~ surg24, d0[keep, ], w_ate[keep])[1:3], te$del_rd_trimmed_on_true_ps_trial0),
  overlap_ato = c(wls_rd(delirium ~ surg24, d0, w_ato)[1:3], te$del_rd_ato_trial0))
colnames(t_ext) <- c("rd", "rd_lo", "rd_hi", "target_rd")
cat("Delirium risk difference by weighting choice, observational part, simulated data\n")
print(round(t_ext, 4))
# moving the nested-trial result to a target population: precision cost of the sampling weights
tr_row <- function(w, target) {
  r <- wls_rd(delirium ~ surg24, tr, w)
  c(r[1:4], se_ratio = r[[4]] / se_trial, ess = ess(w), max_weight = max(w), target_rd = target)
}
t_tr <- rbind(trial = tr_row(rep(1, nrow(tr)), te$del_rd_trial_participants),
              ipsw_whole_registry = tr_row(w_ipsw, te$del_rd_whole_population),
              inverse_odds_non_participants = tr_row(w_iosw, te$del_rd_obs_part_trial0))
colnames(t_tr)[1:4] <- c("rd", "rd_lo", "rd_hi", "se")
cat("Nested trial (n = 950): delirium risk difference moved to a target population, simulated data\n")
print(round(t_tr, 4))
# one-year death with loss to follow-up (rows with known vital status at 12 months)
t_cens <- rbind(
  naive_complete_case = c(wls_rd(died ~ surg24, d0[obs, ], rep(1, sum(obs)))[1:3], te$death_assoc_trial0_rd),
  ipcw = c(wls_rd(died ~ surg24, d0[obs, ], ipcw[obs])[1:3], te$death_assoc_trial0_rd),
  iptw_x_ipcw_ate = c(wls_rd(died ~ surg24, d0[obs, ], (w_ate * ipcw)[obs])[1:3], te$death_rd_ate_trial0))
colnames(t_cens) <- c("rd", "rd_lo", "rd_hi", "target_rd")
cat("One-year death, observational part, simulated data\n")
print(round(t_cens, 4))
cat("Simulated data (ข้อมูลจำลอง): not evidence about any real patient.\n")

# ---- definition of early surgery, registry facts, positivity tail and trimming (simulated data) ----
# early surgery (surg24 = 1) means surgery within this many hours of admission
canon_n("design.surgery_window_hours", 24)
# registry facts (all 20,000 rows): mean age and the share with delirium
canon("reg.age_mean", mean(d$age)); canon("reg.delirium_risk", mean(d$delirium))
# loss to follow-up with and without dementia, unrounded from the simulation's settings file (not published)
canon6 <- function(key, x) cat(sprintf("CANON w1.%s %.6f\n", key, x))
canon6("truth.registry.lost_dementia1", tj$registry_summary$lost_dementia1)
canon6("truth.registry.lost_dementia0", tj$registry_summary$lost_dementia0)
# share of the observational part with an estimated propensity score below 0.05 (revised model)
canon("ps.frac_below_005", mean(ps_cor < 0.05))
# trimming to 0.1 <= e(X) <= 0.9: how many patients leave the analysis, and their share
canon_n("n.trim_removed", sum(!keep)); canon("n.trim_removed_frac", mean(!keep))

> cat("Covariate balance in the observational part, simulated data\n")
Covariate balance in the observational part, simulated data

> print(round(bal, 3))
         smd_before smd_main smd_revised vr_before vr_main vr_revised
age_c        -0.486   -0.046       0.007     0.647   0.589      1.029
age_c2       -0.318   -0.408       0.024     0.589   0.453      1.056
female        0.045    0.001       0.001     0.961   0.999      0.999
frail        -0.516   -0.061      -0.006     0.820   0.977      0.998
dementia     -0.535   -0.071      -0.002     0.509   0.920      0.998
fd           -0.677   -0.290      -0.007     0.153   0.491      0.989
asa3         -0.356   -0.012       0.003     1.131   1.004      0.999
anticoag     -0.402    0.000      -0.016     0.448   1.000      0.973

* Simulated hip-fracture registry of older adults: surgery within 24 hours (surg24), delirium, one-year
* death, a small nested randomised trial (trial = 1) and a delirium risk model with albumin partly missing.
* Simulated data: not evidence about any real drug or patient.
* The simulated registry files (not published) are read from a folder two levels above this script.
version 18
clear all
set more off
set linesize 120
set seed 202610
* number of imputations m: at least the percentage of incomplete rows (albumin is missing in about 30 percent)
global MI_M 40

* ---- commands this file relies on (all ship with Stata) ----
foreach c in teffects tebalance glm margins nlcom stcox streg mi roctab {
    capture which `c'
    display "VERIFY `c' " cond(_rc == 0, "available", "missing")
}

* ---- helpers: print each result on its own CANON line, rounded to 4 decimals ----
program define canon
    gettoken key 0 : 0
    display "CANON w1.`key' " strtrim(string(`0', "%20.4f"))
end
program define canonn
    gettoken key 0 : 0
    display "CANON w1.`key' " strtrim(string(`0', "%20.0f"))
end
* estimate and 95% CI from the last regression: normal (z) or t with the model's residual df
program define ciz
    args key coef tr
    scalar sc_b = _b[`coef']
    scalar sc_se = _se[`coef']
    scalar sc_q = invnormal(0.975)
    if "`tr'" == "t" {
        scalar sc_q = invttail(e(df_r), 0.025)
    }
    if "`tr'" == "exp" {
        canon `key' exp(scalar(sc_b))
        canon `key'.lo exp(scalar(sc_b) - scalar(sc_q)*scalar(sc_se))
        canon `key'.hi exp(scalar(sc_b) + scalar(sc_q)*scalar(sc_se))
    }
    else {
        canon `key' scalar(sc_b)
        canon `key'.lo scalar(sc_b) - scalar(sc_q)*scalar(sc_se)
        canon `key'.hi scalar(sc_b) + scalar(sc_q)*scalar(sc_se)
    }
end
* Kish effective sample size, (sum w)^2 / sum w^2 (left in scalar sc_ess as well)
program define essprint
    syntax varname [if], key(string)
    tempvar sq
    quietly generate double `sq' = `varlist'^2 `if'
    quietly summarize `varlist' `if', meanonly
    scalar sc_s1 = r(sum)
    quietly summarize `sq' `if', meanonly
    scalar sc_ess = scalar(sc_s1)^2 / r(sum)
    canon `key' scalar(sc_ess)
end
* risk ratio from a weighted log-link Poisson model, robust SE treating the weights as known
program define wpoisrr
    syntax varname [if], key(string)
    quietly glm delirium surg24 [pw = `varlist'] `if', family(poisson) link(log) vce(robust)
    glm
    ciz `key' surg24 exp
end
* risk ratio of the two weighted risks after teffects ipw ..., pomeans (delta method on the log ratio,
* joint M-estimation variance that accounts for estimating e(X))
program define pomrr
    args key
    matrix B = e(b)
    matrix V = e(V)
    matrix G = (-1/B[1,1], 1/B[1,2])
    matrix VG = G*V[1..2,1..2]*G'
    scalar sc_b = ln(B[1,2]/B[1,1])
    scalar sc_se = sqrt(VG[1,1])
    canon `key' exp(scalar(sc_b))
    canon `key'.lo exp(scalar(sc_b) - invnormal(0.975)*scalar(sc_se))
    canon `key'.hi exp(scalar(sc_b) + invnormal(0.975)*scalar(sc_se))
end
* simulated truth: the values the simulation was built to produce, read from its settings file (not
* published) and printed as w1.truth.<name>; an optional anchor names the block the value sits in
mata:
real scalar truthval(string scalar f, string scalar k, | string scalar anchor)
{
    string colvector L
    string scalar pat, s
    real scalar i, p, i0
    L = cat(f)
    i0 = 1
    if (args() < 3) anchor = ""
    if (anchor != "") {
        for (i = 1; i <= rows(L); i++) {
            if (strpos(L[i], char(34) + anchor + char(34) + ":") > 0) {
                i0 = i
                break
            }
        }
    }
    pat = char(34) + k + char(34) + ":"
    for (i = i0; i <= rows(L); i++) {
        p = strpos(L[i], pat)
        if (p > 0) {
            s = substr(L[i], p + strlen(pat), .)
            s = subinstr(s, ",", "", .)
            return(strtoreal(strtrim(s)))
        }
    }
    return(.)
}
end
program define truthecho
    local tj ../../datasets/W1/truth.json
    foreach nm of local 0 {
        mata: st_numscalar("sc_tv", truthval("`tj'", "`nm'"))
        canon truth.`nm' scalar(sc_tv)
    }
end
* one value inside a named block of the settings file: truthat <block> <name in the block> <printed name>
program define truthat
    args block nm key
    local tj ../../datasets/W1/truth.json
    mata: st_numscalar("sc_tv", truthval("`tj'", "`nm'", "`block'"))
    canon truth.`key' scalar(sc_tv)
end
* calibration of a linear predictor lp on the outcome: the slope is the coefficient of lp in a logistic
* regression of delirium on lp; calibration-in-the-large (CITL) is the intercept when lp enters as an
* offset (slope fixed at 1). Ideal values: slope 1, CITL 0.
program define calfit, rclass
    args lp
    quietly logit delirium `lp'
    return scalar slope = _b[`lp']
    return scalar slope_se = _se[`lp']
    quietly logit delirium, offset(`lp')
    return scalar citl = _b[_cons]
    return scalar citl_se = _se[_cons]
end
* print slope and CITL with normal 95% CIs as <key>.slope, <key>.citl (and .lo, .hi), keep a row for the table
program define calprint
    args key lp sfx
    calfit `lp'
    scalar sc_q = invnormal(0.975)
    scalar sc_s = r(slope)
    scalar sc_sse = r(slope_se)
    scalar sc_c = r(citl)
    scalar sc_cse = r(citl_se)
    matrix CALROW = (scalar(sc_s), scalar(sc_s) - scalar(sc_q)*scalar(sc_sse), scalar(sc_s) + scalar(sc_q)*scalar(sc_sse), scalar(sc_c), scalar(sc_c) - scalar(sc_q)*scalar(sc_cse), scalar(sc_c) + scalar(sc_q)*scalar(sc_cse))
    calcanon `key' "`sfx'"
end
* the CANON lines and the table row from the 1 x 6 matrix CALROW (slope, lo, hi, CITL, lo, hi)
program define calcanon
    args key sfx
    canon `key'.slope`sfx' CALROW[1,1]
    canon `key'.slope.lo`sfx' CALROW[1,2]
    canon `key'.slope.hi`sfx' CALROW[1,3]
    canon `key'.citl`sfx' CALROW[1,4]
    canon `key'.citl.lo`sfx' CALROW[1,5]
    canon `key'.citl.hi`sfx' CALROW[1,6]
    matrix CAL = nullmat(CAL) \ CALROW
end
* Rubin's rules over the m rows of a matrix holding (slope, SE, CITL, SE) per imputed dataset: pooled
* estimate, total variance W + (1 + 1/m) B, and a t-based 95% CI with the large-sample Rubin df; the
* result goes to CALROW (slope, lo, hi, CITL, lo, hi)
mata:
void rubin_cal(string scalar nm)
{
    real matrix X
    real rowvector out
    real scalar j, M, q, w, b, t, df, c
    X = st_matrix(nm)
    M = rows(X)
    out = J(1, 6, .)
    for (j = 0; j <= 1; j++) {
        q = sum(X[., 2*j + 1]) / M
        w = sum(X[., 2*j + 2]:^2) / M
        b = sum((X[., 2*j + 1] :- q):^2) / (M - 1)
        t = w + (1 + 1/M)*b
        df = (M - 1)*(1 + w/((1 + 1/M)*b))^2
        c = invttail(df, 0.025)
        out[3*j + 1] = q
        out[3*j + 2] = q - c*sqrt(t)
        out[3*j + 3] = q + c*sqrt(t)
    }
    st_matrix("CALROW", out)
}
end
* linear predictor of the complete-case delirium model with a chosen albumin variable
program define lpcc
    args newvar albvar
    generate double `newvar' = scalar(cc__cons) + scalar(cc_age_c)*age_c + scalar(cc_female)*female + scalar(cc_frail)*frail + scalar(cc_dementia)*dementia + scalar(cc_asa3)*asa3 + scalar(cc_anticoag)*anticoag + scalar(cc_albumin)*`albvar'
end
* apparent development AUROC of a pooled MI model (posted by mi estimate, post), averaged over the m completed sets
program define mi_devauc
    args key
    tempvar lp0 lpm
    quietly generate double `lp0' = _b[_cons] + _b[age_c]*age_c + _b[female]*female + _b[frail]*frail + _b[dementia]*dementia + _b[asa3]*asa3 + _b[anticoag]*anticoag
    scalar sc_auc = 0
    forvalues m = 1/$MI_M {
        capture drop `lpm'
        quietly generate double `lpm' = `lp0' + _b[albumin]*_`m'_albumin
        quietly roctab delirium `lpm'
        scalar sc_auc = scalar(sc_auc) + r(area)/$MI_M
    }
    canon `key' scalar(sc_auc)
end

* ---- data: the simulated registry file (not published) ----
import delimited using ../../datasets/W1/W1.csv, clear varnames(1) asdouble
generate double age_c = age - 80
generate double age_c2 = age_c^2
generate double fd = frail*dementia
tempfile full
quietly save `full'
canonn n _N
quietly count if trial == 1
canonn n.trial r(N)
quietly count if trial == 0
canonn n.obs r(N)
quietly count if missing(albumin)
scalar sc_nmiss = r(N)
canonn n.albumin_missing scalar(sc_nmiss)
* share of registry rows with albumin missing, and the number of imputations m used below
canon n.albumin_missing_frac scalar(sc_nmiss)/_N
canonn e1.mi_m $MI_M
* the simulated validation file (not published)
preserve
import delimited using ../../datasets/W1/W1_validation.csv, clear varnames(1) asdouble
canonn n.validation _N
restore

* observational part: treatment chosen by clinicians
keep if trial == 0
quietly count if surg24 == 1
canonn n.obs.treated r(N)
quietly count if died == 1
canonn n.obs.died r(N)
quietly count if lost == 1
canonn n.obs.lost r(N)

* ---- crude comparison of delirium ----
regress delirium surg24, vce(robust)
ciz crude.del.rd surg24 t

* ---- propensity score models: main effects only versus the true form (tight tolerances) ----
logit surg24 age_c female frail dementia asa3 anticoag, tolerance(1e-10) ltolerance(1e-12) nrtolerance(1e-10)
predict double ps_mis, pr
logit surg24 age_c age_c2 female frail dementia fd asa3 anticoag, tolerance(1e-10) ltolerance(1e-12) nrtolerance(1e-10)
predict double ps_cor, pr
generate double w_raw = 1
generate double w_mis = cond(surg24 == 1, 1/ps_mis, 1/(1 - ps_mis))
generate double w_cor = cond(surg24 == 1, 1/ps_cor, 1/(1 - ps_cor))
generate double w_att = cond(surg24 == 1, 1, ps_cor/(1 - ps_cor))
generate double w_ato = cond(surg24 == 1, 1 - ps_cor, ps_cor)
summarize ps_cor
canon ps.min r(min)
canon ps.max r(max)
quietly count if ps_cor < 0.05
canonn ps.n_below_005 r(N)

* ---- standardised mean differences: weighted means, unweighted pooled SD in the denominator ----
foreach lab in raw mis cor {
    scalar sc_max = 0
    foreach x of varlist age_c age_c2 female frail dementia fd asa3 anticoag {
        quietly summarize `x' if surg24 == 1
        scalar sc_v1 = r(Var)
        quietly summarize `x' if surg24 == 0
        scalar sc_v0 = r(Var)
        quietly summarize `x' [aw = w_`lab'] if surg24 == 1
        scalar sc_m1 = r(mean)
        quietly summarize `x' [aw = w_`lab'] if surg24 == 0
        scalar sc_m0 = r(mean)
        scalar sc_smd = (sc_m1 - sc_m0)/sqrt((sc_v1 + sc_v0)/2)
        canon smd.`lab'.`x' scalar(sc_smd)
        scalar sc_max = max(sc_max, abs(sc_smd))
    }
    canon smd.`lab'.maxabs scalar(sc_max)
}

* ---- weights: raw, stabilised, truncated at the 1st and 99th percentiles ----
summarize surg24, meanonly
scalar sc_pa = r(mean)
generate double w_sw = w_cor*cond(surg24 == 1, sc_pa, 1 - sc_pa)
_pctile w_cor, p(1 99)
scalar sc_p01 = r(r1)
scalar sc_p99 = r(r2)
canon wt.trunc.p01 scalar(sc_p01)
canon wt.trunc.p99 scalar(sc_p99)
generate double w_trunc = min(max(w_cor, sc_p01), sc_p99)
foreach lab in cor sw trunc {
    local nm = cond("`lab'" == "cor", "raw", "`lab'")
    summarize w_`lab', meanonly
    canon wt.`nm'.max r(max)
    essprint w_`lab', key(wt.`nm'.ess)
    essprint w_`lab' if surg24 == 1, key(wt.`nm'.ess_treated)
    essprint w_`lab' if surg24 == 0, key(wt.`nm'.ess_control)
}

* ---- IPTW effects on delirium: teffects (M-estimation SE that accounts for estimating e(X)) ----
teffects ipw (delirium) (surg24 age_c age_c2 female frail dementia fd asa3 anticoag, logit), ate
matrix list e(b)
ciz ipw.del.ate ATE:r1vs0.surg24
* this interval accounts for estimating e(X); its SE is printed for comparison with the weights-known SE below
canon ipw.del.ate.se _se[ATE:r1vs0.surg24]
tebalance summarize
teffects ipw (delirium) (surg24 age_c age_c2 female frail dementia fd asa3 anticoag, logit), atet
ciz ipw.del.att ATET:r1vs0.surg24
teffects ipw (delirium) (surg24 age_c female frail dementia asa3 anticoag, logit), ate
ciz ipw.del.ate_mis ATE:r1vs0.surg24
tebalance summarize
* the two weighted risks and their ratio (delta method on the log ratio, joint M-estimation variance)
teffects ipw (delirium) (surg24 age_c age_c2 female frail dementia fd asa3 anticoag, logit), pomeans
matrix list e(b)
matrix B = e(b)
canon ipw.del.ate.risk1 B[1,2]
canon ipw.del.ate.risk0 B[1,1]
pomrr ipw.del.ate.rr
* misspecified score: the same risk ratio and M-estimation SE under the main-effects-only propensity model
teffects ipw (delirium) (surg24 age_c female frail dementia asa3 anticoag, logit), pomeans
pomrr ipw.del.ate_mis.rr
* ATO, stabilised, truncated and trimmed: weighted regression, robust SE treating weights as known
* (each weighted fit runs quietly and is then replayed, which shows the same table without the weight-sum note)
* raw ATE weights through the same weights-known estimator, so raw and stabilised compare like for like
quietly regress delirium surg24 [pw = w_cor], vce(robust)
regress
ciz ipw.del.ate_rawreg surg24 t
* weights-known robust SE, raw weights
canon ipw.del.ate_rawreg.se _se[surg24]
quietly regress delirium surg24 [pw = w_ato], vce(robust)
regress
ciz ipw.del.ato surg24 t
quietly regress delirium surg24 [pw = w_sw], vce(robust)
regress
ciz ipw.del.ate_sw surg24 t
* weights-known robust SE, stabilised weights
canon ipw.del.ate_sw.se _se[surg24]
quietly regress delirium surg24 [pw = w_trunc], vce(robust)
regress
ciz ipw.del.ate_trunc surg24 t
* trimming changes the population: keep 0.1 <= e(X) <= 0.9
quietly count if ps_cor >= 0.1 & ps_cor <= 0.9
canonn n.trim r(N)
quietly regress delirium surg24 [pw = w_cor] if ps_cor >= 0.1 & ps_cor <= 0.9, vce(robust)
regress
ciz ipw.del.ate_trim surg24 t
* risk ratios beside each weights-known risk difference (log-link Poisson, robust SE, weights known)
wpoisrr w_cor, key(ipw.del.ate_rawreg.rr)
wpoisrr w_sw, key(ipw.del.ate_sw.rr)
wpoisrr w_trunc, key(ipw.del.ate_trunc.rr)
wpoisrr w_cor if ps_cor >= 0.1 & ps_cor <= 0.9, key(ipw.del.ate_trim.rr)
wpoisrr w_ato, key(ipw.del.ato.rr)
* the true values these delirium estimates are compared with (simulated truth)
truthecho del_rd_ate_trial0 del_rr_ate_trial0 del_rd_att_trial0 del_rr_att_trial0 del_rd_ato_trial0 del_rr_ato_trial0
truthecho del_rd_trimmed_on_true_ps_trial0 del_rr_trimmed_on_true_ps_trial0 del_assoc_trial0_rd del_assoc_trial0_rr

* ---- 1-year death: crude Kaplan-Meier risk and crude Cox ----
stset fu_months, failure(died)
sts generate km = s, by(surg24)
sts generate kmse = se(s), by(surg24)
summarize km if surg24 == 0 & _t >= 12, meanonly
scalar sc_r0 = 1 - r(mean)
summarize kmse if surg24 == 0 & _t >= 12, meanonly
scalar sc_se0 = r(mean)
summarize km if surg24 == 1 & _t >= 12, meanonly
scalar sc_r1 = 1 - r(mean)
summarize kmse if surg24 == 1 & _t >= 12, meanonly
scalar sc_se1 = r(mean)
canon km.risk0 scalar(sc_r0)
canon km.risk1 scalar(sc_r1)
scalar sc_rd = sc_r1 - sc_r0
scalar sc_se = sqrt(sc_se0^2 + sc_se1^2)
canon crude.death.rd scalar(sc_rd)
canon crude.death.rd.lo scalar(sc_rd) - invnormal(0.975)*scalar(sc_se)
canon crude.death.rd.hi scalar(sc_rd) + invnormal(0.975)*scalar(sc_se)
stcox surg24, nohr
ciz crude.death.hr surg24 exp

* ---- inverse probability of censoring weights: exponential model for loss to follow-up ----
stset fu_months, failure(lost)
streg dementia age_c surg24, distribution(exponential) nohr
predict double xb_c, xb
* probability of still being followed at one's own end time, given X
generate double ipcw = 1/exp(-exp(xb_c)*fu_months)
summarize ipcw if lost == 0
canon ipcw.max r(max)
regress died surg24 if lost == 0, vce(robust)
ciz naive.death.rd surg24 t
quietly regress died surg24 [pw = ipcw] if lost == 0, vce(robust)
regress
ciz ipcw.death.rd surg24 t
* truth for the IPCW-only contrast: censoring-free ASSOCIATIONAL risks by arm in trial = 0 (no confounding control)
truthecho death_assoc_trial0_risk1 death_assoc_trial0_risk0 death_assoc_trial0_rd death_assoc_trial0_rr
* truth for the IPTW x IPCW contrasts below: the causal 1-year risk difference and risk ratio (ATE, trial = 0)
truthecho death_rd_ate_trial0 death_rr_ate_trial0
* ---- IPTW times IPCW: the 1-year risk difference for ATE, ATT and ATO ----
foreach lab in cor att ato {
    local nm = cond("`lab'" == "cor", "ate", "`lab'")
    generate double wd_`lab' = w_`lab'*ipcw
    quietly regress died surg24 [pw = wd_`lab'] if lost == 0, vce(robust)
    regress
    ciz ipw.death.`nm' surg24 t
}
* ---- IPTW Cox model (ATE weights; pweights give a robust SE) ----
stset fu_months [pw = w_cor], failure(died)
stcox surg24, nohr
ciz ipw.death.hr surg24 exp

* ---- link and family: saturated versus adjusted fits (all registry rows) ----
use `full', clear
* saturated fits (surg24 only): log-binomial, modified Poisson, and Gaussian family with a log link
glm delirium surg24, family(binomial) link(log)
* model-based CI (identical bread in both languages for a saturated model)
ciz c2.rr_crude_logbin surg24 exp
* log-RR SE, log-binomial, model-based
canon c2.se_crude_logbin _se[surg24]
quietly glm delirium surg24, family(poisson) link(log)
scalar sc_pn = _se[surg24]
glm delirium surg24, family(poisson) link(log) vce(robust)
matrix bpc = e(b)
* modified Poisson, robust CI
ciz c2.rr_crude_poisson surg24 exp
* Poisson model-based log-RR SE (too large for a binary outcome) versus the sandwich log-RR SE
canon c2.se_crude_poisson_naive scalar(sc_pn)
canon c2.se_crude_poisson_robust _se[surg24]
* naive Poisson 95% CI
canon c2.rr_crude_poisson_naive.lo exp(_b[surg24] - invnormal(0.975)*scalar(sc_pn))
canon c2.rr_crude_poisson_naive.hi exp(_b[surg24] + invnormal(0.975)*scalar(sc_pn))
* Gaussian log link, robust CI (saturated: same bread in both languages)
glm delirium surg24, family(gaussian) link(log) vce(robust) from(bpc)
ciz c2.rr_crude_gaussian surg24 exp
* adjusted for age and sex
quietly glm delirium surg24 age_c female, family(poisson) link(log)
scalar sc_p2n = _se[surg24]
glm delirium surg24 age_c female, family(poisson) link(log) vce(robust)
matrix bp2 = e(b)
scalar sc_p2 = exp(_b[surg24])
scalar sc_p2lo = exp(_b[surg24] - invnormal(0.975)*_se[surg24])
scalar sc_p2hi = exp(_b[surg24] + invnormal(0.975)*_se[surg24])
scalar sc_p2r = _se[surg24]
* adjusted fits: the log-binomial starts from the Poisson fit
glm delirium surg24 age_c female, family(binomial) link(log) from(bp2)
canon c2.rr_adj2_logbin exp(_b[surg24])
* Stata's ML glm uses the observed-information SE for the non-canonical log link (R's glm: expected)
canon c2.rr_adj2_logbin.lo.stata exp(_b[surg24] - invnormal(0.975)*_se[surg24])
canon c2.rr_adj2_logbin.hi.stata exp(_b[surg24] + invnormal(0.975)*_se[surg24])
* log-RR SE, adjusted log-binomial (observed information)
canon c2.se_adj2_logbin.stata _se[surg24]
canon c2.rr_adj2_poisson scalar(sc_p2)
* modified Poisson, robust CI; model-based versus sandwich log-RR SE, adjusted
canon c2.rr_adj2_poisson.lo scalar(sc_p2lo)
canon c2.rr_adj2_poisson.hi scalar(sc_p2hi)
canon c2.se_adj2_poisson_naive scalar(sc_p2n)
canon c2.se_adj2_poisson_robust scalar(sc_p2r)
* Gaussian log link, adjusted for age and sex (observed-information bread here)
glm delirium surg24 age_c female, family(gaussian) link(log) vce(robust) from(bp2)
canon c2.rr_adj2_gaussian exp(_b[surg24])
canon c2.rr_adj2_gaussian.lo.stata exp(_b[surg24] - invnormal(0.975)*_se[surg24])
canon c2.rr_adj2_gaussian.hi.stata exp(_b[surg24] + invnormal(0.975)*_se[surg24])
* adjusted for age and frailty: the Poisson fit gives fitted risks above 1, so the log-binomial cannot start
quietly glm delirium surg24 age_c frail, family(poisson) link(log)
scalar sc_pafn = _se[surg24]
glm delirium surg24 age_c frail, family(poisson) link(log) vce(robust)
matrix bpaf = e(b)
ciz c2.rr_adjaf_poisson surg24 exp
* model-based versus sandwich log-RR SE, age and frailty
canon c2.se_adjaf_poisson_naive scalar(sc_pafn)
canon c2.se_adjaf_poisson_robust _se[surg24]
predict double mu_af, mu
summarize mu_af
* largest fitted risk: above 1 means the log-binomial wall
canon c2.adjaf_poisson_max_fitted r(max)
capture noisily glm delirium surg24 age_c frail, family(binomial) link(log) from(bpaf) iterate(100)
scalar sc_afconv = 0
if _rc == 0 {
    scalar sc_afconv = e(converged)
}
* 1 = converged
canonn c2.adjaf_logbin_converged.stata scalar(sc_afconv)

* ---- the risk-ratio ladder for a common outcome (all registry rows, same covariates) ----
* rung 1: log-binomial; it stops when a covariate pattern would need a risk above 1
capture noisily glm delirium surg24 age_c female frail dementia asa3 anticoag, family(binomial) link(log) iterate(100)
scalar sc_rc = _rc
scalar sc_conv = 0
if sc_rc == 0 {
    scalar sc_conv = e(converged)
}
display "NOTE log-binomial return code " sc_rc
canonn c3.logbin_converged scalar(sc_conv)
* rung 2: modified Poisson with robust SE
glm delirium surg24 age_c female frail dementia asa3 anticoag, family(poisson) link(log) vce(robust) eform
ciz c3.rr_poisson surg24 exp
matrix bp = e(b)
predict double mu_p, mu
summarize mu_p
canon c3.poisson_max_fitted r(max)
quietly count if mu_p > 1
canonn c3.poisson_n_fitted_above1 r(N)
* rung 3: Gaussian family with a LOG link and robust SE (starts from the Poisson fit)
glm delirium surg24 age_c female frail dementia asa3 anticoag, family(gaussian) link(log) vce(robust) from(bp) eform
* Stata's ML glm uses the observed-information bread here (non-canonical link), so its robust CI differs from R's
canon c3.rr_gaussian exp(_b[surg24])
canon c3.rr_gaussian.lo.stata exp(_b[surg24] - invnormal(0.975)*_se[surg24])
canon c3.rr_gaussian.hi.stata exp(_b[surg24] + invnormal(0.975)*_se[surg24])
* rung 4: logistic model, then marginal standardisation for the risk ratio and risk difference
logit delirium i.surg24 age_c female frail dementia asa3 anticoag
ciz c5.or_cond 1.surg24 exp
margins surg24, post
nlcom (rd: _b[1.surg24] - _b[0.surg24]) (lnrr: ln(_b[1.surg24]/_b[0.surg24])) (lnor: ln((_b[1.surg24]/(1 - _b[1.surg24]))/(_b[0.surg24]/(1 - _b[0.surg24])))), post
ciz c3.rr_std lnrr exp
ciz c3.rd_std rd
* the same model averaged over the cohort gives the marginal odds ratio (non-collapsibility)
ciz c5.or_marg lnor exp
* the realised nested trial (n = 950): one noisy draw. Chance covariate imbalance in a trial this small can
* outweigh non-collapsibility, so these keys are named "realised"; the large-sample pair is printed further down
logit delirium surg24 if trial == 1
ciz c5.trial_realised.or_crude surg24 exp
logit delirium i.surg24 age_c female frail dementia asa3 anticoag if trial == 1
ciz c5.trial_realised.or_adj 1.surg24 exp
* trial-standardised marginal OR: the adjusted trial model averaged over the trial participants
margins surg24, post
nlcom (lnor: ln((_b[1.surg24]/(1 - _b[1.surg24]))/(_b[0.surg24]/(1 - _b[0.surg24])))), post
ciz c5.trial_realised.or_std lnor exp
* large-sample simulated truth: conditional OR within frailty strata versus marginal ORs
* in trial participants, and the large-sample value of the six-covariate adjusted model in the trial
truthecho c5_or_cond_nonfrail c5_or_cond_frail c5_or_marg_trial_nonfrail c5_or_marg_trial_frail c5_or_marg_trial
truthecho c5_or_adj_pseudo_trial del_rr_whole_population del_rd_whole_population

* ---- prediction model for delirium with albumin missing: complete case versus MI without and with Y ----
logit delirium age_c female frail dementia asa3 anticoag albumin
canonn e1.cc.n e(N)
ciz e1.cc.b_albumin albumin
matrix b_cc = e(b)
foreach v in age_c female frail dementia asa3 anticoag albumin _cons {
    scalar cc_`v' = _b[`v']
}
predict double lp_ccdev if e(sample), xb
roctab delirium lp_ccdev
scalar auc_dev_cc = r(area)
* apparent AUROC on the complete-case development rows
canon e1.cc.dev_auroc scalar(auc_dev_cc)
* deployment-style imputation model for albumin fitted on the development data, no outcome in it
regress albumin age_c female frail dementia asa3 anticoag
matrix b_ri = e(b)
keep id delirium age_c female frail dementia asa3 anticoag albumin
tempfile dev
quietly save `dev'
mi set wide
mi register imputed albumin
mi register regular delirium age_c female frail dementia asa3 anticoag
* imputation model WITHOUT the outcome
mi impute chained (pmm, knn(10)) albumin = age_c female frail dementia asa3 anticoag, add($MI_M) rseed(202610)
mi estimate, post: logit delirium age_c female frail dementia asa3 anticoag albumin
matrix T = r(table)
canon e1.mi_noy.b_albumin.stata T[1,7]
canon e1.mi_noy.b_albumin.lo.stata T[5,7]
canon e1.mi_noy.b_albumin.hi.stata T[6,7]
matrix b_noy = e(b)
* apparent development AUROC of the pooled model, averaged over the m completed development sets
mi_devauc e1.mi_noy.dev_auroc.stata
scalar auc_dev_noy = sc_auc
use `dev', clear
mi set wide
mi register imputed albumin
mi register regular delirium age_c female frail dementia asa3 anticoag
* imputation model WITH the outcome
mi impute chained (pmm, knn(10)) albumin = age_c female frail dementia asa3 anticoag delirium, add($MI_M) rseed(202610)
mi estimate, post: logit delirium age_c female frail dementia asa3 anticoag albumin
matrix T = r(table)
canon e1.mi_y.b_albumin.stata T[1,7]
canon e1.mi_y.b_albumin.lo.stata T[5,7]
canon e1.mi_y.b_albumin.hi.stata T[6,7]
matrix b_y = e(b)
* apparent development AUROC of the pooled model, averaged over the m completed development sets
mi_devauc e1.mi_y.dev_auroc.stata
scalar auc_dev_y = sc_auc
* discrimination and calibration of each fitted model in the simulated validation file (albumin fully observed)
import delimited using ../../datasets/W1/W1_validation.csv, clear varnames(1) asdouble
generate double age_c = age - 80
matrix score double lp_cc = b_cc
matrix score double lp_noy = b_noy
matrix score double lp_y = b_y
roctab delirium lp_cc
scalar auc_val_cc = r(area)
canon e1.cc.auroc scalar(auc_val_cc)
roctab delirium lp_noy
scalar auc_val_noy = r(area)
canon e1.mi_noy.auroc.stata scalar(auc_val_noy)
roctab delirium lp_y
scalar auc_val_y = r(area)
canon e1.mi_y.auroc.stata scalar(auc_val_y)
* AUROC table: development rows (apparent) and the validation file, one row per fitted model
matrix AUC = (scalar(auc_dev_cc), scalar(auc_val_cc) \ scalar(auc_dev_noy), scalar(auc_val_noy) \ scalar(auc_dev_y), scalar(auc_val_y))
matrix rownames AUC = cc mi_noy mi_y
matrix colnames AUC = development validation
matrix list AUC, format(%9.4f) title(AUROC of each fitted model, simulated data)
* calibration slope and calibration-in-the-large of each fitted model, normal 95% CIs
calprint e1.cc.cal lp_cc
calprint e1.mi_noy.cal lp_noy .stata
calprint e1.mi_y.cal lp_y .stata
* validation with missing albumin: the validation file has albumin fully observed, so a reproducible MAR mask
* is applied in memory (the file is not changed). The mask uses the missingness model the registry was
* simulated with and a deterministic uniform u = frac(id x 0.6180339887498949), identical in R.
generate double u_mask = mod(id*0.6180339887498949, 1)
generate double p_mask = invlogit(-1.57 + 0.9*dementia + 0.5*asa3 + 0.02*age_c + 0.3*frail)
generate double alb_m = albumin
replace alb_m = . if u_mask < p_mask
quietly count if missing(alb_m)
* validation rows whose albumin is masked
canonn e1.val.n_masked r(N)
* one fixed model is scored: the complete-case fit; albumin fully observed first (same as e1.cc.auroc)
lpcc lpv_full albumin
roctab delirium lpv_full
scalar auc_m_full = r(area)
canon e1.val.auroc_full scalar(auc_m_full)
calprint e1.val.cal_full lpv_full
* Y-free regression imputation fitted on the development data (deployment-style)
matrix score double alb_hat = b_ri
generate double alb_ri = cond(missing(alb_m), alb_hat, alb_m)
lpcc lpv_ri alb_ri
roctab delirium lpv_ri
scalar auc_m_regimp = r(area)
canon e1.val.auroc_regimp scalar(auc_m_regimp)
* single imputation: these calibration CIs ignore the uncertainty of the filled values
calprint e1.val.cal_regimp lpv_ri
* multiple imputation inside the validation sample, with versus without the validation outcomes
keep id delirium age_c female frail dementia asa3 anticoag alb_m
tempfile valm
quietly save `valm'
foreach wy in with without {
    use `valm', clear
    local yv = cond("`wy'" == "with", "delirium", "")
    mi set wide
    mi register imputed alb_m
    mi register regular delirium age_c female frail dementia asa3 anticoag
    mi impute chained (pmm, knn(10)) alb_m = age_c female frail dementia asa3 anticoag `yv', add($MI_M) rseed(202611)
    scalar sc_auc = 0
    matrix CALM = J($MI_M, 4, .)
    forvalues m = 1/$MI_M {
        quietly lpcc lpm _`m'_alb_m
        quietly roctab delirium lpm
        scalar sc_auc = scalar(sc_auc) + r(area)/$MI_M
        calfit lpm
        matrix CALM[`m', 1] = r(slope)
        matrix CALM[`m', 2] = r(slope_se)
        matrix CALM[`m', 3] = r(citl)
        matrix CALM[`m', 4] = r(citl_se)
        drop lpm
    }
    * validation AUROC averaged over the m imputations, validation outcomes `wy' in the imputation model
    canon e1.val.auroc_mi_`wy'_y.stata scalar(sc_auc)
    scalar auc_m_`wy' = sc_auc
    * slope and CITL pooled over the m imputations with Rubin's rules (t-based 95% CI)
    mata: rubin_cal("CALM")
    calcanon e1.val.cal_mi_`wy'_y .stata
}
* AUROC of the fixed complete-case model after the masked validation albumin was filled four ways
matrix AUCM = (scalar(auc_m_full) \ scalar(auc_m_regimp) \ scalar(auc_m_with) \ scalar(auc_m_without))
matrix rownames AUCM = full regimp mi_y mi_noy
matrix colnames AUCM = auroc
matrix list AUCM, format(%9.4f) title(AUROC after filling masked validation albumin, simulated data)
* calibration table: the three fitted models on the validation file, then the fixed complete-case model
* after the masked validation albumin was filled four ways (slope 1 and CITL 0 mean perfect calibration);
* val_mi_y and val_mi_noy: multiple imputation in the validation file with and without the validation outcomes
matrix rownames CAL = cc mi_noy mi_y val_full val_regimp val_mi_y val_mi_noy
matrix colnames CAL = slope slope_lo slope_hi citl citl_lo citl_hi
matrix list CAL, format(%9.4f) title(Calibration on the validation set, simulated data)
* the reference value for the albumin coefficient: this working model fitted to a very large simulated
* population with every albumin value recorded, and that reference model's AUROC on the validation file
truthat pseudo_true_coefficients albumin pseudo_true_albumin
truthecho auroc_pseudo_true_on_validation

* ---- nested trial: sampling-score weights to move the trial result to a target population ----
use `full', clear
logit trial age_c female frail dementia asa3 anticoag, tolerance(1e-10) ltolerance(1e-12) nrtolerance(1e-10)
predict double s_hat, pr
* target: the whole registry
generate double w_ipsw = 1/s_hat if trial == 1
* target: the non-participants (inverse odds of participation)
generate double w_iosw = (1 - s_hat)/s_hat if trial == 1
regress delirium surg24 if trial == 1, vce(robust)
ciz a2.trial.rd surg24 t
* robust SE of the unweighted trial difference
scalar sc_setr = _se[surg24]
canon a2.trial.rd.se scalar(sc_setr)
quietly regress delirium surg24 [pw = w_ipsw] if trial == 1, vce(robust)
regress
ciz a2.ipsw.rd surg24 t
* robust SE after weighting to the whole registry, and the variance cost: IPSW SE over the trial SE
canon a2.ipsw.rd.se _se[surg24]
canon a2.ipsw.se_ratio _se[surg24]/scalar(sc_setr)
quietly regress delirium surg24 [pw = w_iosw] if trial == 1, vce(robust)
regress
ciz a2.iosw.rd surg24 t
* robust SE after weighting to the non-participants, and its ratio to the trial SE
canon a2.iosw.rd.se _se[surg24]
canon a2.iosw.se_ratio _se[surg24]/scalar(sc_setr)
quietly count if trial == 1
scalar sc_ntr = r(N)
summarize w_ipsw, meanonly
canon a2.ipsw.max r(max)
essprint w_ipsw if trial == 1, key(a2.ipsw.ess)
* ESS as a share of the 950 trial participants
canon a2.ipsw.ess_frac scalar(sc_ess)/scalar(sc_ntr)
summarize w_iosw, meanonly
canon a2.iosw.max r(max)
essprint w_iosw if trial == 1, key(a2.iosw.ess)
canon a2.iosw.ess_frac scalar(sc_ess)/scalar(sc_ntr)
* risk ratios beside the transported risk differences (log-link Poisson, robust SE, weights known)
generate double w_one = 1
wpoisrr w_one if trial == 1, key(a2.trial.rr)
wpoisrr w_ipsw if trial == 1, key(a2.ipsw.rr)
wpoisrr w_iosw if trial == 1, key(a2.iosw.rr)
* their targets: trial participants, the whole registry population, the non-participants (trial = 0)
truthecho del_rd_trial_participants del_rr_trial_participants del_rd_obs_part_trial0 del_rr_obs_part_trial0
display "All numbers above are from simulated data (ข้อมูลจำลอง)."

. * ---- propensity score models: main effects only versus the true form (tight tolerances) ----
. logit surg24 age_c female frail dementia asa3 anticoag, tolerance(1e-10) ltolerance(1e-12) nrtolerance(1e-10)

Iteration 0:  Log likelihood = -12976.056
Iteration 1:  Log likelihood = -11381.133
Iteration 2:  Log likelihood = -11362.234
Iteration 3:  Log likelihood = -11362.195
Iteration 4:  Log likelihood = -11362.195
Iteration 5:  Log likelihood = -11362.195

Logistic regression                                    Number of obs =  19,050
                                                       LR chi2(6)    = 3227.72
                                                       Prob > chi2   =  0.0000
Log likelihood = -11362.195                            Pseudo R2     =  0.1244

------------------------------------------------------------------------------
      surg24 | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       age_c |  -.0201935   .0015235   -13.25   0.000    -.0231796   -.0172074
      female |   .0855296   .0350511     2.44   0.015     .0168307    .1542285
       frail |  -.6993138   .0345164   -20.26   0.000    -.7669647    -.631663
    dementia |  -1.025651   .0407161   -25.19   0.000    -1.105453   -.9458486
        asa3 |  -.4517041   .0332777   -13.57   0.000    -.5169273   -.3864809
    anticoag |  -1.189063    .047972   -24.79   0.000    -1.283087    -1.09504
       _cons |    .592086    .038296    15.46   0.000     .5170272    .6671449
------------------------------------------------------------------------------

. predict double ps_mis, pr

. logit surg24 age_c age_c2 female frail dementia fd asa3 anticoag, tolerance(1e-10) ltolerance(1e-12) nrtolerance(1e-10
> )

Iteration 0:  Log likelihood = -12976.056
Iteration 1:  Log likelihood = -10984.436
Iteration 2:  Log likelihood = -10905.294
Iteration 3:  Log likelihood = -10904.587
Iteration 4:  Log likelihood = -10904.586
Iteration 5:  Log likelihood = -10904.586

Logistic regression                                    Number of obs =  19,050
                                                       LR chi2(8)    = 4142.94
                                                       Prob > chi2   =  0.0000
Log likelihood = -10904.586                            Pseudo R2     =  0.1596

------------------------------------------------------------------------------
      surg24 | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       age_c |    -.02911   .0016808   -17.32   0.000    -.0324042   -.0258157
      age_c2 |  -.0024792   .0001106   -22.41   0.000     -.002696   -.0022624
      female |   .0885217   .0357782     2.47   0.013     .0183977    .1586456
       frail |  -.3874434   .0389039    -9.96   0.000    -.4636936   -.3111933
    dementia |  -.3835303   .0534843    -7.17   0.000    -.4883576    -.278703
          fd |  -1.583228   .0917753   -17.25   0.000    -1.763105   -1.403352
        asa3 |  -.4718794   .0339289   -13.91   0.000    -.5383789   -.4053799
    anticoag |  -1.220435   .0483475   -25.24   0.000    -1.315194   -1.125675
       _cons |   .7777036   .0415885    18.70   0.000     .6961916    .8592156
------------------------------------------------------------------------------