Extreme Weights: What Stabilizing, Truncating and Trimming Really Change
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
Inverse probability of treatment weighting (IPTW) gives each patient a weight of one over the probability of the treatment they received, so a patient treated against the odds can carry a weight of 50 and dominate the estimate. Four common responses change different things: stabilising, capping weights (truncation), removing patients (trimming) and overlap weights. In a point-treatment analysis, stabilising multiplies every weight within an arm by the same constant, so an extreme patient stays extreme relative to everyone else in that arm. Stabilisation matters most in marginal structural models, where weights are multiplied over many time points. Truncation trades bias for variance; trimming changes the population the estimate describes. Overlap weights are bounded and target the patients for whom both choices were plausible. In a simulated hip-fracture registry, truncating the weights at the 1st and 99th percentiles moved the delirium risk difference from -0.028 to -0.075, about twice the true -0.038. This article concludes that every fix should be reported with the estimand it targets.
One patient with a weight of 50
A hip-fracture team is reviewing a registry analysis of early surgery, meaning an operation within 24 hours of admission, and postoperative delirium. The analyst has sorted the weights and stopped at one row. A frail woman with dementia, on anticoagulants, was operated on early, although patients like her usually wait. Her propensity score is 0.02, so her weight is $1/0.02 = 50$ and she counts as fifty patients in her arm's weighted risk.
The team asks which of the usual fixes to apply: stabilise the weights, cap them (truncation), remove patients like her (trimming), or switch to another kind of weight such as overlap weights. Each one changes something different: the scale of the weights, the balance of bias and variance, the population described, or the estimand, the precise quantity the analysis sets out to estimate. This article takes the fixes one at a time, first on a five-patient hand example and then on a simulated registry.
The woman and her weight of 50 form the hand example, where she carries almost all of her arm's weight. In the simulated registry the largest weight is bigger still, and it is the many large weights together that cost the early-surgery arm most of its information.
Why some weights explode
A weight becomes large only when a patient received the treatment that was unlikely for them. A treated patient with $e(X)$ near 0 gets $1/e(X)$, and an untreated patient with $e(X)$ near 1 gets $1/(1 - e(X))$; both are huge. A treated patient with $e(X)$ near 1, or an untreated one near 0, has a weight close to 1. In the weighted analysis, a patient with a huge weight stands in for the many similar patients who received the other treatment.
Huge weights signal a near-violation of positivity, the requirement that every combination of covariates has a real chance of either treatment, $0 < e(X) < 1$ [1]. A near-violation is a property of the data and the question, not a software fault. If frail patients with dementia almost never have early surgery, the data hold little information about what early surgery does for them. Weighting can only reuse the patients who were observed; it cannot create comparable ones who were not.
The simulated registry
A simulated hip-fracture registry shows the problem at scale (simulated data). It holds 19,050 patients whose timing was chosen by their clinicians: 8,053 had early surgery and 10,997 waited. A further 950 patients in a small randomised substudy of the same registry are left out, because their timing was not chosen by clinicians. The outcome $Y$ is postoperative delirium, coded 1 if it occurred.
The propensity model is a logistic regression on age, age squared, sex, frailty, dementia, a frailty-by-dementia term, an American Society of Anesthesiologists (ASA) physical status of 3 or more, and anticoagulant use. That model has the same form as the one that generated the data, so what follows is not caused by a wrong model.
The estimated scores run from 0.0045 to 0.721, and 1,198 patients (6.3%) score below 0.05. In the true scores the lowest value is 0.0042, and 6.2% fall below 0.05. The lowest score among early-surgery patients, 0.0072, produces the largest weight in the registry, 138.5. With a correct model, extreme weights are the data's own positivity problem.
Effective sample size: what the weights are worth
One number summarises how much a set of weights costs in information, the effective sample size (ESS):
$$\mathrm{ESS} = \frac{\left(\sum_i w_i\right)^2}{\sum_i w_i^2}$$The sums run over the patients of one arm, and $w_i$ is the weight of patient $i$. With equal weights the ESS equals the number of patients, and the more unequal the weights, the smaller it becomes. It is computed per arm, because each arm's weighted risk rests on that arm's weights alone.
In the registry, the 8,053 early-surgery patients have an ESS of 2,997 under the raw weights. The 10,997 who waited keep an ESS of 9,539. The early-surgery arm holds the extreme weights, so it carries most of the loss.
Stabilisation: same patients, a new scale
A stabilised weight replaces the 1 in the numerator with the overall share of patients who received the same treatment [2]:
$$sw_i = \frac{P(A = a_i)}{P(A = a_i \mid X_i)}$$Here $a_i$ is the treatment patient $i$ actually received, so a treated patient gets $P(A = 1)/e(X_i)$ and an untreated patient gets $P(A = 0)/(1 - e(X_i))$. In the registry $P(A = 1)$ is 42.3% and $P(A = 0)$ is 57.7%, so the largest weight falls from 138.5 to 58.5. But every early-surgery weight was multiplied by the same number, so the patient with the largest weight still outweighs everyone else in that arm by the same factor.
Unstabilised ATE weights add up to about twice the number of patients, because each arm is weighted up to the size of the whole sample. Stabilised weights add up to about the original number [3]. Stabilisation changes the scale of the weights, and so the output of software that treats weights as patient counts, but not who dominates.
Treatment in this registry is decided once, at admission, so this is a point-treatment analysis. A marginal structural model (MSM) instead describes outcomes under treatment strategies that can change over time, with a weight built at every time point and multiplied along the way.
Why the estimate does not move
A weighted regression of the outcome on treatment alone estimates each arm's risk as a normalised weighted mean, the weighted share of patients with delirium:
$$\hat{\mu}_a = \frac{\sum_{i:\,A_i = a} w_i Y_i}{\sum_{i:\,A_i = a} w_i}$$Here $\hat{\mu}_a$ is the estimated risk under treatment $a$, and the sums run over the patients who received $a$. Multiplying every weight in that arm by one constant multiplies the top and the bottom by it, so $\hat{\mu}_a$ does not move. The same constant cancels from the arm's ESS, and the risk difference (RD), $\hat{\mu}_1 - \hat{\mu}_0$, is unchanged as well. This holds for a weighted regression on treatment alone; add covariates to the weighted model and rescaling the arms by different constants can move the estimate.
Hand example: stabilising five weights
Hand example. Take five early-surgery patients with weights 1, 1, 1, 1 and 50. The four weights of 1 stand for patients who were almost certain to have early surgery, rounded for easy arithmetic; the fifth is the woman from the opening scene, with $e(X) = 0.02$. In this small example half of all patients had early surgery, so $P(A = 1) = 0.50$.
-
Her raw weight
\[ w = \frac{1}{e(X)} = \frac{1}{0.02} = 50 \]
She stands in for many similar patients who waited.
-
ESS of the raw weights
\[ \mathrm{ESS} = \frac{(1 + 1 + 1 + 1 + 50)^2}{1^2 + 1^2 + 1^2 + 1^2 + 50^2} = \frac{54^2}{2504} = \frac{2916}{2504} = 1.16 \]
Five patients carry the information of 1.16 equally weighted patients.
-
Her stabilised weight
\[ sw = P(A = 1) \times \frac{1}{e(X)} = 0.50 \times 50 = 25 \]
Each of the other four becomes 0.5.
-
ESS of the stabilised weights
\[ \mathrm{ESS} = \frac{(4 \times 0.5 + 25)^2}{4 \times 0.5^2 + 25^2} = \frac{27^2}{626} = \frac{729}{626} = 1.16 \]
Unchanged: the factor 0.50 cancels from the top and the bottom.
-
Her share of the arm's total weight
\[ \frac{50}{54} = \frac{25}{27} \]
Her share is the same before and after stabilising.
Result: Stabilising halved every weight in this arm and left the ESS at 1.16. The extreme patient is exactly as dominant as before.
What stabilisation did in the registry
A robust (sandwich) standard error is computed from the spread of the data around the fitted model rather than from the model's own variance formula. A weighted regression of delirium on early surgery alone uses such a standard error, treating the weights as known numbers. With the raw weights it gives an RD of -0.028 (95% CI -0.052 to -0.005) and a standard error of 0.0120. With the stabilised weights the RD, the interval and the standard error are identical, and so is the ESS of each arm.
Stata's teffects ipw command, run by the same script shown in the code section below, reports the same RD with a narrower interval, -0.046 to -0.010, and a standard error of 0.0091. That interval also accounts for the uncertainty in estimating $e(X)$, which here makes it narrower; its width is not an effect of stabilisation. The ESS of both arms pooled does shift after stabilising, but only because the two arms are multiplied by different constants. Report the ESS per arm, and compare weighting choices on one variance estimator.
Where stabilisation matters: treatment over time
When treatment is decided more than once, for example a drug given or withheld at several time points in intensive care, each patient's weight is a product of one factor per time point. Unstabilised products can grow very large after only a few time points. The stabilised numerator, the probability of each treatment decision given past treatment and baseline covariates only, keeps the product closer to 1 and narrows the spread of the weights [2]; the marginal structural model must then include those baseline covariates. The series part on marginal structural models works through that case.
Truncation: capping the weights
Weight truncation caps the weights at chosen values, often percentiles of the weight distribution such as the 1st and 99th [2]. Weights above the upper cap are set to it, and weights below the lower cap are raised to it. Every patient stays in the analysis, but the extreme ones lose influence. Truncation trades bias for variance: the variance falls, but a capped patient no longer stands in for everyone they represented.
Some papers, including a simulation study of weighting methods, call this step weight trimming [4]. In this article trimming means removing patients, as in the next section, so check which one a paper means.
In the five-patient hand example, a cap of 10 turns the weights 1, 1, 1, 1 and 50 into 1, 1, 1, 1 and 10:
$$\mathrm{ESS} = \frac{(4 + 10)^2}{4 \times 1^2 + 10^2} = \frac{14^2}{104} = \frac{196}{104} = 1.88$$The ESS rises from 1.16 to 1.88. The cost is that patients with the woman's profile are now under-represented in the weighted early-surgery arm.
What truncation did in the registry
In the registry, the 1st and 99th percentiles of all the weights, both arms together, are 1.02 and 8.10. The highest estimated score is 0.721, so no later-surgery weight comes near the upper cap, and that arm barely changes. In the early-surgery arm the ESS rises from 2,997 to 5,937. The RD, however, moves from -0.028 to -0.075 (95% CI -0.091 to -0.059), about twice the true RD of -0.038, which is known because the data are simulated.
The raw estimate sat a little on the near-zero side of the truth. The truncated one lands beyond it on the other side, and its interval no longer contains the true value. The direction follows from who was capped: early-surgery patients with low scores, who are more often frail with dementia and have a high risk of delirium. Capping them shrinks that high-risk group inside the weighted early-surgery arm and re-creates part of the confounding the weights had removed.
Here truncation bought precision at the price of a much larger error, a poor trade.
The weights themselves: largest weight and ESS per arm
| ATE weights | Largest weight | ESS, early surgery (patients, of 8,053) | ESS, later surgery (patients, of 10,997) |
|---|---|---|---|
| Raw | 138.5 | 2,997 | 9,539 |
| Stabilised | 58.5 | 2,997 | 9,539 |
| Truncated at the 1st and 99th percentiles | 8.1 | 5,937 | 9,540 |
Trimming: changing who is analysed
Trimming removes patients whose propensity score lies outside a chosen range, for example 0.1 to 0.9 [5], or beyond chosen percentiles of the score in each arm [6]. The patients who remain all had a real chance of either treatment, so their weights stay moderate. The estimate, though, now describes only them, and that population must be described in the paper.
In the registry, keeping scores from 0.1 to 0.9 removes 2,190 patients (11.5%) and keeps 16,860. Because the highest score is 0.721, every removed patient sits at the low end, among the profiles for which early surgery was rare. The RD is -0.030 (95% CI -0.046 to -0.015). Its target is the true RD in the trimmed population, -0.042, not the ATE of -0.038.
One caveat applies to that true value. The simulation trims on the true propensity score while the analysis trims on the estimated one, so the two trimmed populations are close but not identical.
Overlap weights: a bounded weight and a new question
Overlap weights give a treated patient the weight $1 - e(X)$ and an untreated patient the weight $e(X)$: the probability of the treatment they did not receive [7]. Every overlap weight lies between 0 and 1, so no weight can grow without limit the way $1/e(X)$ or $1/(1 - e(X))$ can; the woman from the opening scene gets $1 - 0.02 = 0.98$ instead of 50. When $e(X)$ comes from a logistic regression, overlap weights balance the mean of every covariate in that model exactly between the arms [7].
The price is a new estimand. Overlap weights target the ATO, the average treatment effect in the overlap population, in which each patient counts in proportion to $e(X)\,(1 - e(X))$ [8]. That population leans toward patients for whom both choices were plausible. The woman's own weight is near the top of the range, but the many patients with her profile who waited get weights of only 0.02, so her profile counts for little.
In the registry the overlap-weighted RD is -0.033 (95% CI -0.046 to -0.019), against a true ATO of -0.043. Its interval contains its own target, and that target is not the ATE.
Five estimates, three questions
The table below collects the registry estimates, from all 19,050 patients or, after trimming, the 16,860 kept, each beside the true value for its own estimand. All five come from the same weighted regression with the same robust standard error, so their intervals are comparable.
| Weighting choice | Estimand (who the answer is about) | RD | 95% CI | True RD for that estimand |
|---|---|---|---|---|
| Raw ATE weights | ATE: the population the 19,050 patients represent | -0.028 | -0.052 to -0.005 | -0.038 |
| Stabilised ATE weights | ATE: the population the 19,050 patients represent | -0.028 | -0.052 to -0.005 | -0.038 |
| Truncated at the 1st and 99th percentiles | ATE, as intended | -0.075 | -0.091 to -0.059 | -0.038 |
| Trimmed to scores 0.1 to 0.9 | The trimmed population (16,860 patients kept) | -0.030 | -0.046 to -0.015 | -0.042 (trimmed on the true score) |
| Overlap weights | ATO: the overlap population | -0.033 | -0.046 to -0.019 | -0.043 |
Reading the five estimates
In this registry the three true values lie closer together than the sampling error of any one estimate. Apart from truncation, the gaps between the estimates are therefore within sampling error; only the truncated estimate lies well outside that range, for the reason given above.
A change of estimand is still a change of question, even when, as here, the answers happen to be close. In other data the targets can lie far apart.
The code that produced the table follows, in Stata and in R. The registry file held 20,000 rows; 950 were in the randomised substudy, left out here, leaving 19,050.
Stata: raw ATE weights through a robust standard error
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
. * 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
Linear regression Number of obs = 19,050
F(1, 19048) = 5.47
Prob > F = 0.0193
R-squared = 0.0009
Root MSE = .46651
------------------------------------------------------------------------------
| Robust
delirium | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.0281481 .0120344 -2.34 0.019 -.0517366 -.0045596
_cons | .3346588 .0045191 74.06 0.000 .325801 .3435165
------------------------------------------------------------------------------R: largest weight, ESS per arm and the RD for every weighting choice
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("ATE weights: largest weight and effective sample size per arm, simulated data\n")
ATE weights: largest weight and effective sample size per arm, simulated data
> print(round(t_w, 4))
max ess_early ess_later
raw 138.4954 2997.349 9539.221
stabilised 58.5461 2997.349 9539.221
truncated 8.0983 5937.391 9539.829
> 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")
Delirium risk difference by weighting choice, observational part, simulated data
> print(round(t_ext, 4))
rd rd_lo rd_hi target_rd
raw -0.0281 -0.0517 -0.0046 -0.0384
stabilised -0.0281 -0.0517 -0.0046 -0.0384
truncated_p1_p99 -0.0749 -0.0908 -0.0590 -0.0384
trimmed_0.1_0.9 -0.0302 -0.0457 -0.0147 -0.0424
overlap_ato -0.0327 -0.0463 -0.0191 -0.0426Four fixes, four different changes
| Method | What it changes | Estimand | Bias | Variance | When to consider it |
|---|---|---|---|---|---|
| Stabilisation | The scale of each arm's weights | Unchanged | Unchanged | Unchanged for normalised means with a robust standard error | Usually reasonable as a default (see the box on what stabilisation does and does not do) |
| Truncation | The largest and smallest weights | Still the ATE in name | Can grow; here the estimate moved further from the truth | Falls | As a sensitivity analysis at several prespecified caps |
| Trimming | Who is analysed | The trimmed population | Judged against the new population | Usually falls | When the question can be limited to patients with a real choice |
| Overlap weights | Every weight, now between 0 and 1 | ATO | Judged against the ATO | Usually low, since no weight is extreme | When the question is about patients for whom both choices were plausible |
Reporting: one estimand per number
The five estimates in the table answer three different questions. Raw, stabilised and truncated weights target the ATE, trimming targets the trimmed population, and overlap weights target the ATO. Each number belongs beside the name of its estimand.
For inverse probability weights, whether raw, stabilised or capped, also report the largest weight and the ESS of each arm, as the weights table does. Showing several choices side by side lets a reader see whether a conclusion rests on a few patients.
After any change to the weights, check covariate balance again; the balance-checking part of this series shows how. The series hub builds the weights themselves from the start.
Common misreadings and their fixes
-
"Stabilised weights fix extreme weights."
The largest weight shrinks in size, from 138.5 to 58.5 in the registry, but its share of its arm does not change.
Fix: In a point-treatment analysis, stabilising multiplies every weight within an arm by the same constant, so an extreme patient stays extreme relative to everyone else in that arm. Stabilisation matters most in marginal structural models, where weights are multiplied over many time points. Truncation trades bias for variance; trimming changes the population the estimate describes.
-
"Stabilising narrowed the confidence interval."
In the registry the narrower interval comes from Stata's teffects ipw, which also accounts for estimating e(X). On the same variance estimator, raw and stabilised weights give the same standard error, 0.0120.
Fix: Compare weighting choices on one variance estimator, and label an interval that accounts for estimating e(X) as such.
-
"Truncation only removes noise."
Capping changes how much the capped patients count, and in the registry it moved the estimate away from the truth.
Fix: Report the caps and the estimate at several caps, and read a large shift as a warning about positivity, not as a better estimate.
-
"Trimming is a sensitivity analysis of the same effect."
Trimming removes patients, so the estimate describes a different population. In the registry it removed 2,190 patients, all from the low end of the score.
Fix: Describe who was removed and who remains, and name the estimand of the trimmed analysis.
-
"Overlap weights are always better."
They remove extreme weights by changing the question to the ATO.
Fix: Overlap weights answer a question about the patients for whom both choices were plausible. If the question is about the whole population, they answer a different question.
-
"Extreme weights prove the propensity model is wrong."
The registry's propensity model has the correct form, and its weights are still extreme.
Fix: Check the model, but read extreme weights first as a positivity problem in the data and the question.
What to do in your own analysis
- Before opening the outcome, plot the propensity score by arm and list the largest weights with the patients behind them.
- For each set of inverse probability weights, report the largest weight and the ESS of each arm, not only the pooled ESS.
- Stabilised weights are a sensible default, but expect them to change little at a single time point.
- If you truncate, consider prespecifying the caps and reporting the estimate at several of them; a large shift is a warning about positivity rather than a better estimate.
- If you trim or use overlap weights, describe the population that remains and name its estimand wherever the result appears.
- Compare weighting choices on one variance estimator, and label an interval that accounts for estimating $e(X)$ as such.
Glossary
- propensity score
- The probability of receiving treatment given baseline covariates, written e(X).
- positivity
- The requirement that every covariate pattern has a real chance of each treatment.
- pseudo-population
- The weighted population in which treatment no longer depends on the measured covariates.
- effective sample size
- The number of equally weighted patients that would carry the same information as a weighted arm: the squared sum of the weights over the sum of the squared weights.
- stabilised weight
- An inverse probability weight multiplied by the overall probability of the treatment received.
- normalised weighted mean
- A weighted sum of outcomes divided by the sum of the weights in that arm, so rescaling the arm's weights leaves it unchanged.
- marginal structural model
- A model for outcomes under treatment strategies, fitted with inverse probability weights, usually when treatment changes over time.
- weight truncation
- Capping weights at chosen values, often percentiles, which trades bias for variance.
- trimming
- Removing patients whose propensity score lies outside a chosen range, which changes the population the estimate describes.
- overlap weight
- A weight equal to the probability of the treatment not received, always between 0 and 1.
- ATO
- The average treatment effect in the overlap population, where each patient counts in proportion to e(X)(1 - e(X)).
- estimand
- The precise quantity an analysis sets out to estimate, including the population it describes.
- robust (sandwich) standard error
- A standard error computed from the spread of the data around the fitted model rather than from the model's own variance formula.
References
- Petersen ML, Porter KE, Gruber S, Wang Y, van der Laan MJ. Diagnosing and responding to violations in the positivity assumption. Stat Methods Med Res. 2012;21(1):31-54. doi:10.1177/0962280210386207 https://doi.org/10.1177/0962280210386207
- Cole SR, Hernán MA. Constructing inverse probability weights for marginal structural models. Am J Epidemiol. 2008;168(6):656-664. doi:10.1093/aje/kwn164 https://doi.org/10.1093/aje/kwn164
- Xu S, Ross C, Raebel MA, Shetterly S, Blanchette C, Smith D. Use of stabilized inverse propensity scores as weights to directly estimate relative risk and its confidence intervals. Value Health. 2010;13(2):273-277. doi:10.1111/j.1524-4733.2009.00671.x https://doi.org/10.1111/j.1524-4733.2009.00671.x
- Lee BK, Lessler J, Stuart EA. Weight trimming and propensity score weighting. PLoS One. 2011;6(3):e18174. doi:10.1371/journal.pone.0018174 https://doi.org/10.1371/journal.pone.0018174
- Crump RK, Hotz VJ, Imbens GW, Mitnik OA. Dealing with limited overlap in estimation of average treatment effects. Biometrika. 2009;96(1):187-199. doi:10.1093/biomet/asn055 https://doi.org/10.1093/biomet/asn055
- Stürmer T, Rothman KJ, Avorn J, Glynn RJ. Treatment effects in the presence of unmeasured confounding: dealing with observations in the tails of the propensity score distribution - a simulation study. Am J Epidemiol. 2010;172(7):843-854. doi:10.1093/aje/kwq198 https://doi.org/10.1093/aje/kwq198
- Li F, Morgan KL, Zaslavsky AM. Balancing covariates via propensity score weighting. J Am Stat Assoc. 2018;113(521):390-400. doi:10.1080/01621459.2016.1260466 https://doi.org/10.1080/01621459.2016.1260466
- Li F, Thomas LE, Li F. Addressing extreme propensity scores via the overlap weights. Am J Epidemiol. 2019;188(1):250-257. doi:10.1093/aje/kwy201 https://doi.org/10.1093/aje/kwy201
Key takeaways
- Extreme weights come from patients who received a treatment that was unlikely for them, a near-violation of positivity that a correct propensity model does not remove.
- Stabilising shrank the largest weight in the registry but left the ESS of each arm, and the estimate and robust standard error of a weighted regression on treatment alone, unchanged.
- Truncation trades bias for variance, and in the simulated registry it moved the risk difference to about twice its true value.
- Trimming and overlap weights tame extreme weights by changing the population, so each answers its own named question.
- Report every weighting choice side by side with its estimand, and every set of inverse probability weights with its largest weight and the ESS of each arm.
Related in the wiki: [[iptw-guide]] [[propensity-weighting-balance-check-and-model-revision]]
0