When Weighting Leaves Imbalance: Fix the Model Before You Blame the Data
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
When a balance table shows large differences after inverse probability of treatment weighting, the data often get blamed. In 19,050 simulated hip-fracture patients whose surgery timing followed routine care, early (within 24 hours) and later surgery are compared for delirium. A main-terms propensity model (a logistic model of each patient's chance of early surgery) enters each covariate once and linearly. It left age squared and the frailty-by-dementia product imbalanced, with absolute standardised mean differences (SMDs, mean gaps in pooled standard deviations) reaching 0.408. Its weighted risk difference, -0.088, had a confidence interval excluding -0.038, the simulated average treatment effect of early surgery for all versus none. Adding both terms brought every absolute SMD to 0.024 or less and the estimate to -0.028. Because SMDs compare only means, variance ratios, distribution plots and product terms are checked too. This article shows how to read a balance table, tell model misspecification from a positivity problem, revise the model blind to the outcome and report it.
A balance table that will not settle
An analyst is preparing a study from a national hip-fracture registry of older adults. The question is whether surgery within 24 hours of admission, called early surgery here, lowers the risk of postoperative delirium compared with later surgery. Clinicians chose the timing, and frail patients, patients with dementia and patients on anticoagulants tended to wait. A crude comparison therefore flatters early surgery.
The analyst fits a propensity model, weights the patients and opens the balance table. Most rows look reassuring, but two do not settle: the spread of age, captured by a squared age term, and the share of patients with both frailty and dementia. A colleague suggests that the groups are simply too different to compare. This article argues that the propensity model is the first suspect, and shows how to repair it without looking at the outcome.
What a balance table shows
A balance table compares the treatment groups before and after weighting. It has one row per covariate, plus a row for each square or product worth checking, whether or not that term is in the propensity model. It shows each group's mean for a continuous covariate such as age, and each group's proportion for a binary one such as dementia. Two summary columns follow: the standardised mean difference, which compares means, and the variance ratio, which compares spread.
Checking these tables is called balance diagnostics. The weights are judged by the balance they produce in the measured covariates, before any outcome is analysed [1]. The same checks were set out for propensity-score matching, which pairs each treated patient with an untreated patient of similar score [2], and weighting uses weighted versions of them [1].
The standardised mean difference
The standardised mean difference (SMD) of a covariate is
$$\mathrm{SMD} = \frac{\bar x_1 - \bar x_0}{\sqrt{(s_1^2 + s_0^2)/2}}$$Here $\bar x_1$ and $\bar x_0$ are the covariate's means in the treated and untreated groups, and $s_1^2$ and $s_0^2$ are its variances. For a binary covariate the mean is a proportion $p$ and the variance is $p(1 - p)$. The SMD is the difference in means in units of the pooled standard deviation (SD), so age in years and a proportion with dementia share one scale.
Unlike a P value, the SMD does not depend on sample size. As a registry grows, the same small difference gives an ever smaller P value while its SMD stays put [1, 2]. After weighting, $\bar x_1$ and $\bar x_0$ become weighted means.
This article, like the R code below, keeps the unweighted pooled SD in the denominator before and after weighting, so the SMD moves only when the means move. Stata's tebalance summarize, like Austin and Stuart [1], uses weighted SDs after weighting, so its weighted SMDs can differ. It also lists only the terms of the fitted model, so after a model without squares or products, age squared and frailty by dementia must be checked separately.
An absolute SMD below 0.1 is a widely used convention for adequate balance [1, 2]. It is not a statistical test and not a guarantee.
Hand example: what a 12-year age gap means
A balance table reports mean ages of 86 and 74 years in two surgery groups, with an SD of 12 years in each. Let subscript 1 label the older group; swapping the labels flips only the sign of the SMD.
-
Gap in means
\[ \bar x_1 - \bar x_0 = 86 - 74 = 12 \]
The groups differ by 12 years.
-
Pooled SD
\[ \sqrt{(12^2 + 12^2)/2} = \sqrt{(144 + 144)/2} = 12 \]
With equal SDs, the pooled SD is the common SD.
-
Standardised mean difference
\[ \mathrm{SMD} = 12 / 12 = 1.00 \]
The gap is one full SD, ten times the 0.1 convention.
-
A quick check on a printed SMD
\[ s_p = 12 / 0.45 = 26.7 \]
A quick check on any printed SMD: divide the gap in means by it to recover the pooled SD $s_p$ it implies, and compare that with the SDs the table reports. An SMD of 0.45 for this same 12-year gap would imply a pooled SD of 26.7 years, more than twice the 12 years reported.
-
Against the simulated registry
\[ 12 / 11.6 = 1.03 \]
Across the whole simulated registry described below, the SD of age is 11.6 years, so the same gap is about one SD there too.
Result: A 12-year age gap among older hip-fracture patients is an SMD of about 1.0, a very large imbalance. A printed SMD of 0.45 for that gap fails the quick check: the implied SD of 26.7 years is more than twice the 11.6 years in the simulated registry, so the printed value would be wrong.
What the SMD misses
An SMD compares one number per group, the mean. Two groups can share a mean and still differ in three ways.
- Spread. The variance ratio, $\mathrm{VR} = s_1^2 / s_0^2$, uses weighted variances after weighting. A value near 1 means similar spread. Ratios as far from 1 as one half or two have been called far too extreme [3], and one half to two is often quoted as an outer limit [4]. The VR is read mainly for continuous covariates [2], because a binary covariate's variance, $p(1 - p)$, is fixed by its proportion. The SMD of a squared term such as age squared offers a second check on spread [2].
- Shape. The empirical cumulative distribution function (ECDF) gives, at each value, the share of a group at or below it. The Kolmogorov-Smirnov (KS) distance, the largest vertical gap between two groups' ECDFs, catches differences in shape that the mean and variance miss [1].
- Joint distribution. Balance on each covariate does not imply balance on their combinations. The SMD of the frailty $\times$ dementia product checks the patients with both conditions [2].
A love plot shows the SMDs at a glance. It has one row per covariate or term, a dot for the absolute SMD before and after weighting, and a vertical line at 0.1.
Two reasons balance can stay poor
When weighting leaves imbalance, two explanations call for different responses, and both can be present at once. The first is a misspecified propensity model: the way covariates drive treatment includes a curve or an interaction that the model leaves out. Its signature is imbalance concentrated in the omitted terms, which shrinks when those terms are added.
The second is a positivity problem. Positivity means that every covariate pattern has a non-zero chance of each treatment. When the remaining imbalance sits among patients whose scores are close to 0 or 1, no added term can create comparable patients in the other group. Part 3 on extreme weights takes that case up.
The simulated registry: two propensity models
The registry is simulated data. It holds 20,000 adults with hip fracture, with a mean age of about 80 years and an SD of 11.6 years. We set aside 950 patients who were enrolled in a separate small trial. The example uses only the other 19,050, whose surgery timing was chosen in routine care: 8,053 with early and 10,997 with later surgery.
Two propensity models were fitted by logistic regression on these 19,050 patients. The main-terms model enters each covariate once, with age as a straight line through age_c, age minus 80 years. Its covariates are age, sex, frailty, dementia, anticoagulant use and ASA grade 3 or more, which marks severe systemic disease on the American Society of Anesthesiologists (ASA) scale. The revised model adds age_c2, the square of age_c, and fd, frailty times dementia.
In the simulation, treatment depends on every main term and also on age squared and on frailty $\times$ dementia. The target is the ATE among patients like these 19,050, whose surgery timing is chosen in routine care. For early surgery on delirium it is a risk difference of -0.038 (risk ratio 0.89). In real data the true form is unknown, so the revision is guided by clinical reasoning and by the balance table.
Balance before and after weighting, by propensity model
| Term | SMD, no weights | SMD, main-terms weights | SMD, revised weights | VR, no weights | VR, main-terms weights | VR, revised weights |
|---|---|---|---|---|---|---|
| Age, linear (age_c) | -0.486 | -0.046 | 0.007 | 0.647 | 0.589 | 1.029 |
| Age squared (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 |
| Frail $\times$ dementia (fd) | -0.677 | -0.290 | -0.007 | 0.153 | 0.491 | 0.989 |
| ASA grade 3 or more | -0.356 | -0.012 | 0.003 | 1.131 | 1.004 | 0.999 |
| Anticoagulant | -0.402 | 0.000 | -0.016 | 0.448 | 1.000 | 0.973 |
Reading the table
Under the main-terms weights the mean of age is balanced (SMD -0.046), but its VR moved further from 1, from 0.647 before weighting to 0.589. The squared age term, the squared distance of each patient's age from 80 years, is largest for the youngest and the oldest patients and shows the problem plainly. Its SMD is -0.408, larger in absolute value than before weighting (-0.318), and its VR is 0.453, beyond the one-half limit. Weighting with an incomplete model can make an omitted term worse.
Frailty and dementia each look balanced, with SMDs of -0.061 and -0.071. Their product, the patients with both conditions, still has an SMD of -0.290. Under the revised weights the largest absolute SMD is 0.024, and every VR lies between 0.973 and 1.056.
R: SMD and variance ratio under each propensity model
R code w1_sim_r.R
# 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)) Output of the run w1_sim_r.log
> 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.973Delirium estimates under each set of weights
| Weights | Largest absolute SMD | Risk difference (95% CI) | Risk ratio (95% CI) | Simulated target: RD (RR) |
|---|---|---|---|---|
| None (crude comparison) | 0.677 | -0.262 (-0.275 to -0.250) | 0.39 (0.37 to 0.41) | Association: -0.264 (0.38) |
| Main-terms model | 0.408 | -0.088 (-0.101 to -0.075) | 0.74 (0.70 to 0.77) | ATE: -0.038 (0.89) |
| Revised model | 0.024 | -0.028 (-0.046 to -0.010) | 0.92 (0.86 to 0.97) | ATE: -0.038 (0.89) |
What the imbalance does to the estimate
A real study would estimate the effect only under the final model; both are shown here because the simulation knows the truth. The main-terms estimate, -0.088, sits between the crude difference (-0.262) and the simulated truth (-0.038), and its 95% confidence interval, -0.101 to -0.075, excludes the truth. Patients with both frailty and dementia carry a much higher simulated delirium risk. They stay over-represented in the later group, which pulls the estimate toward the crude comparison.
The revised estimate is -0.028, with a 95% confidence interval of -0.046 to -0.010 that covers the truth. Its risk ratio is 0.92 against the simulated 0.89, while the main-terms risk ratio is 0.74. The revision worked because the added terms were exactly those in the simulated truth.
Stata: the main-terms and the revised propensity models
Stata code w1_sim.do
* 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 (ข้อมูลจำลอง)."Output of the run w1_sim.log
. * ---- 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
------------------------------------------------------------------------------Revise the model without looking at the outcome
Revising the propensity model is part of the study design, so it happens before any outcome is analysed. If balance is inadequate, the model is changed, for example with nonlinear terms or interactions, and balance is checked again [1, 2]. One pass has four steps.
- List the covariates and terms that remain imbalanced.
- Add squares, splines or interactions of the covariates already in the model that plausibly shape the treatment decision, such as a squared term or a spline for age, or the frailty $\times$ dementia interaction. New variables that predict treatment but not the outcome mainly add variance [4]. A spline lets a covariate's effect bend smoothly at a few chosen values.
- Refit, recompute the weights and recheck balance on every covariate, its square and the plausible interactions.
- Stop when every check meets thresholds set before the loop began, for example an absolute SMD below 0.1 and a VR close to 1 and never beyond one half or two. If no plausible term closes the gap, check where in the score distribution the imbalance sits before going further.
The outcome stays out of sight throughout [3, 4]. If the estimate is visible while terms are tried, the model can drift toward the expected result, knowingly or not. Record each pass: the terms added, the reason and the resulting balance.
Recheck everything after the revision
A new model gives new scores and new weights, so every check is repeated. That means the weighted SMDs, the VRs of continuous terms and the overlap plot, which draws the scores of each group on one axis.
The weights get one more check. The effective sample size, $\mathrm{ESS} = (\sum w)^2 / \sum w^2$, is approximately the number of equally weighted patients that would carry the same information as a weighted group. With the revised weights it is 2,997 of 8,053 early-surgery patients and 9,539 of 10,997 later-surgery patients. A drop as large as the early group's points to highly variable weights, often driven by a few heavy ones, the subject of Part 3.
| Weights | ESS, early surgery | ESS, later surgery |
|---|---|---|
| Revised model | 2,997 of 8,053 | 9,539 of 10,997 |
Methods that aim at balance directly
Two methods build balance into the estimation instead of checking it afterwards. The covariate balancing propensity score (CBPS) fits the logistic propensity model so that weighted covariate means balance, rather than only maximising the likelihood [5]. Entropy balancing chooses the weights directly, so that chosen moments such as means and variances match exactly while the weights stay as close to equal as possible [6].
Both balance only what they are told to, so a product such as frailty $\times$ dementia must still be listed. In its original form, entropy balancing reweights the untreated group toward the treated group. That targets the average treatment effect in the treated (ATT) rather than the ATE.
What to report
A reader can judge the weighting only from what is reported, which usually includes the following.
- The final propensity model and why each added term was plausible.
- The balance table before and after weighting, with SMDs for every covariate, square and interaction checked, and VRs for continuous covariates.
- The threshold used, stated as a convention, for example an absolute SMD below 0.1.
- The number of revisions, and that the outcome was not examined while they were made.
- The ESS in each group, beside the effect estimate.
Common misreadings and their fixes
-
"The groups are too different, so the data cannot be analysed."
Leftover imbalance often reflects the propensity model. In the simulated registry, adding two omitted terms removed it.
Fix: Find the terms that stay imbalanced, add them and recheck. Suspect positivity mainly when the imbalance sits where scores are close to 0 or 1, and check the overlap plot either way.
-
"Mean age is balanced, so age is balanced."
Under the main-terms weights the SMD for age was -0.046, but the VR for age had moved further from 1, and the squared term had an SMD of -0.408 and a VR of 0.453.
Fix: Check the VR and the squared term of every continuous covariate, and the ECDF when shape matters.
-
"Frailty and dementia are each balanced, so patients with both are too."
Each had an absolute SMD below 0.1 under the main-terms weights, yet their product had an SMD of -0.290.
Fix: Check plausible products separately from their components.
-
"Revising the model should make the propensity scores more alike."
A model that captures more of how clinicians chose the timing usually separates the two groups' unweighted scores further. That is the model doing its job.
Fix: Judge the revision by covariate balance after weighting, and read the overlap plot for positivity, not as a score of the model.
-
"The better-fitting model is the better propensity model."
Fit statistics such as the C statistic, which measures how well the model separates treated from untreated patients, describe prediction. The weights exist to balance covariates.
Fix: Judge each version of the model by the balance its weights produce [4].
-
"Try terms until the effect estimate looks plausible."
Choosing the model by its estimate turns the design into a search for a result, and the confidence interval loses its meaning.
Fix: Keep the outcome out of sight until the model is frozen, and record every revision [3, 4].
-
"All SMDs are below 0.1, so the estimate is unbiased."
In the simulated registry the revised weights brought every absolute SMD to 0.024 or less. The estimate landed near the simulated truth only because the simulation has no unmeasured confounder, every patient has a chance of either timing, and the revised model has the true form.
Fix: Small SMDs show balance on the measured covariates in the form checked. They say nothing about unmeasured confounders or about covariate forms that were not checked.
What to do in your own analysis
- Consider writing the balance thresholds and the terms to check, including squares and plausible interactions, into the analysis plan before fitting any propensity model.
- Keep the outcome out of the working dataset until the propensity model is frozen.
- When a term stays imbalanced, try adding it or a spline before concluding that the groups cannot be compared. Check, too, where in the score distribution the imbalance sits.
- Put the balance table, the VRs and the ESS next to the effect estimate in the report.
Glossary
- propensity score (คะแนนแนวโน้มการได้รับการรักษา)
- The probability of receiving treatment given the measured baseline covariates.
- standardised mean difference (ผลต่างค่าเฉลี่ยมาตรฐาน)
- A difference in covariate means between groups, divided by a pooled standard deviation.
- variance ratio
- A covariate's variance in the treated group over its variance in the untreated group.
- balance diagnostics
- Checks that the weighted groups have similar covariate distributions, made before any outcome is analysed.
- love plot
- A dot plot of the absolute SMD of each covariate before and after weighting, with a line at 0.1.
- positivity
- Every covariate pattern has a non-zero chance of each treatment.
- effective sample size (ขนาดตัวอย่างประสิทธิผล)
- Approximately the number of equally weighted patients that would carry the same information as a weighted group.
- covariate balancing propensity score
- A propensity model fitted so that weighted covariate means balance.
- entropy balancing
- Weights chosen so that chosen covariate moments match exactly, with the weights as close to equal as possible.
References
- Austin PC, Stuart EA. Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies. Stat Med. 2015;34(28):3661-3679. doi:10.1002/sim.6607 https://doi.org/10.1002/sim.6607
- Austin PC. Balance diagnostics for comparing the distribution of baseline covariates between treatment groups in propensity-score matched samples. Stat Med. 2009;28(25):3083-3107. doi:10.1002/sim.3697 https://doi.org/10.1002/sim.3697
- Rubin DB. Using propensity scores to help design observational studies: application to the tobacco litigation. Health Serv Outcomes Res Methodol. 2001;2:169-188. doi:10.1023/a:1020363010465 https://doi.org/10.1023/a:1020363010465
- Stuart EA. Matching methods for causal inference: a review and a look forward. Stat Sci. 2010;25(1):1-21. doi:10.1214/09-STS313 https://doi.org/10.1214/09-STS313
- Imai K, Ratkovic M. Covariate balancing propensity score. J R Stat Soc Series B Stat Methodol. 2014;76(1):243-263. doi:10.1111/rssb.12027 https://doi.org/10.1111/rssb.12027
- Hainmueller J. Entropy balancing for causal effects: a multivariate reweighting method to produce balanced samples in observational studies. Polit Anal. 2012;20(1):25-46. doi:10.1093/pan/mpr025 https://doi.org/10.1093/pan/mpr025
Key takeaways
- Imbalance left after weighting points first to the propensity model, not to the data.
- The SMD compares means in pooled-SD units and does not grow with sample size; an absolute SMD below 0.1 is a convention, not a test.
- Variance ratios, distribution plots and the balance of squares and products catch what each covariate's SMD misses.
- Revise the model with the outcome out of sight, stop at thresholds set in advance and record each pass.
- Small SMDs show balance on the measured covariates in the form checked; they say nothing about unmeasured confounders.
Related in the wiki: [[iptw-guide]]
0