The IPW Family: One Idea, Three Kinds of Missing People
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
Inverse probability weighting fills three gaps in clinical data with one move. Each observed patient is weighted to stand in, given covariates, for similar people who were not observed. Treatment weights stand in for the arm a patient did not receive, censoring weights for lost follow-up, and selection weights for people a study never enrolled. In a simulated hip-fracture registry, censoring weights moved the risk difference for one-year death from -0.181, with lost patients dropped, to -0.174, against -0.175 had nobody been lost. Only treatment and censoring weights together, -0.035, came near the causal -0.038. In the registry's nested trial, selection weights re-targeted the delirium result to the whole registry or non-participants, barely changing the estimate, but standard errors grew more than fourfold. This article concludes that each weight needs its own model, its own assumptions and a named target population, and that a selection weight transports an effect only through the effect modifiers (characteristics that change the effect's size) in its model.
A trial inside a registry, and the people nobody observed
A national hip-fracture registry holds 20,000 older adults, all simulated data. It records surgery within 24 hours of admission (early surgery), postoperative delirium and death within one year. Inside the registry runs a randomised trial of early surgery: 950 patients, allocated 1:1 to early or later surgery.
The trial enrolled mostly younger patients without dementia. Across the registry, 28.5% of patients with dementia were lost to follow-up before one year, against 10.3% of those without. The steering committee asks whether the trial result applies to the whole registry, and whether the losses bias the registry's one-year death comparison.
Behind both questions stand people nobody observed. They are each patient under the other surgery timing, the patients lost before one year, and the registry patients who never entered the trial. One weighting idea can stand in for all three groups, if each use gets its own model, its own assumptions and a named target.
One idea: weight by the inverse probability of being observed
Inverse probability weighting (IPW) analyses the people who were observed, and weights each one to stand in for similar people who were not. In its basic form, the weight is one over the probability of being observed in one's own state, given covariates $X$ (characteristics measured at baseline, such as age and dementia):
$$w_i = \frac{1}{P(i \text{ observed in own state} \mid X_i)}$$Here $w_i$ is the weight of person $i$. Someone unlikely to be observed in their state gets a large weight, because few people like them were observed. The same principle underlies weighting for missing data in general [1].
This form targets everyone, observed or not. When the target is only the people who were not observed, the weight becomes the probability of not being observed divided by the probability of being observed. These are the inverse odds used below for selection.
- Treatment: the missing people are the same patients under the other arm.
- Censoring: the missing people were lost before their outcome was seen.
- Selection: the missing people were never enrolled in the study.
The first kind: the treatment a patient did not receive
Inverse probability of treatment weighting (IPTW) is the best-known member. Let $A$ be the treatment, 1 for early surgery and 0 for later surgery. The propensity score $e(X) = P(A = 1 \mid X)$ is the probability of early surgery given $X$. A patient who had early surgery gets the weight $1/e(X)$, and one who had later surgery gets $1/(1 - e(X))$.
In the weighted sample, called a pseudo-population, the measured covariates no longer predict treatment. Its comparison is causal under consistency (the observed outcome is the outcome under the treatment received) and exchangeability (no unmeasured confounding). It also needs positivity (every covariate pattern has a non-zero chance of each treatment) and a correct propensity model. The series hub, IPTW: Rebuilding the Comparison a Trial Would Have Given You, builds these weights in full.
The second kind: follow-up that was lost
A patient lost to follow-up before one year is censored: the one-year outcome is unknown. Let $C_i(t)$ be the censoring indicator, 1 if patient $i$ has been lost by time $t$ and 0 otherwise.
Dropping the patients who were lost biases a one-year risk even when loss is unrelated to the outcome. A patient who dies early is observed before there is much time to be lost, while a survivor must stay in follow-up for the whole year, so the patients kept over-represent deaths. The Kaplan-Meier curve, the standard survival curve, keeps each lost patient in the analysis until the loss and so repairs that first bias. Censoring is informative when the chance of loss depends on something that also predicts the outcome, such as dementia, and informative loss adds a second bias that Kaplan-Meier cannot repair.
Inverse probability of censoring weighting (IPCW) weights each patient whose outcome was seen by one over the probability of still being followed [2]:
$$w^C_i(t) = \frac{1}{P(C_i(t) = 0 \mid X_i, A_i)}$$Here $t$ is the patient's own end of follow-up (death, or one year for survivors), and $A_i$ is the treatment received. Survivors had the whole year in which to be lost, so they get larger weights than similar patients who died early, which repairs the first bias. With dementia in the censoring model, a patient with dementia who stayed stands in for similar patients who were lost, which repairs the second.
Here, with treatment fixed at baseline, censoring is not a confounder. It is a loss of outcome information after treatment, a form of selection bias that weighting corrects when the predictors of loss are measured [3].
One-year death in the simulated registry
Outside the trial the registry holds 19,050 patients, and 2,859 were lost before one year (simulated data). In this simulation, loss depends on dementia only, and dementia also raises the risk of death. The censoring model is an exponential model (a constant rate of loss over time) on dementia, age and early surgery, and the largest censoring weight is 1.64.
Censoring weights remove censoring bias, not confounding. So IPCW alone is judged against the censoring-free association, the difference between the arms had nobody been lost. Treatment times censoring weights are judged against the causal average treatment effect (ATE): early versus later surgery for every patient outside the trial.
The treatment weights come from a logistic propensity model fitted outside the trial. Its form matches how early surgery was generated in this simulation, so the model is correctly specified here by construction. Its covariates are age, age squared, sex, frailty, dementia, frailty times dementia, anticoagulant use and ASA physical status 3 or higher (severe systemic disease or worse).
One-year death, early minus later surgery: four analyses
| Analysis | Risk difference (95% CI) | Truth to compare with |
|---|---|---|
| Complete case: the 2,859 lost patients dropped | -0.181 (-0.194 to -0.168) | -0.175, censoring-free association |
| Kaplan-Meier risks (early 0.157, later 0.324) | -0.167 | -0.175, censoring-free association |
| IPCW | -0.174 (-0.187 to -0.162) | -0.175, censoring-free association |
| IPTW (ATE weights) times IPCW | -0.035 (-0.060 to -0.010) | -0.038, causal effect (ATE) |
The third kind: people the trial never enrolled
A characteristic that changes the size of the treatment effect, on the scale being reported, is an effect modifier. When effect modifiers are distributed differently in a trial and in the population of interest, the trial result need not hold there.
Let $S$ be the participation indicator: 1 for a trial participant, 0 for a member of the source population (here the registry) who did not take part. The sampling score is the probability of taking part given covariates [4, 5]:
$$s(X) = P(S = 1 \mid X)$$It is usually estimated by a logistic regression of $S$ on $X$, fitted to participants and non-participants together. A selection weight then weights each participant by a function of $s(X)$. That function depends on the target population, so name the target first.
Target: the whole source population, participants included. The weight is
$$w^S_i = \frac{1}{s(X_i)}$$so each participant stands for $1/s(X_i)$ people like them. This is inverse probability of selection weighting (IPSW), also called inverse probability of sampling weighting. It serves generalisability: carrying the result to the population the trial was drawn from [4, 6].
Target: the non-participants, or an external population. The weight is the inverse odds of participation:
$$w^S_i = \frac{1 - s(X_i)}{s(X_i)}$$so each participant stands only for the people like them who stayed out. These inverse odds of sampling weights serve transportability: carrying the result to people outside the trial [7]. For an external target, $s(X)$ comes from the trial data and a sample of that population combined in one dataset.
Hand example: one trial, three target populations
A source population of 1,000 holds 500 younger and 500 older patients. The sampling score is 0.80 for a younger and 0.10 for an older patient, so the trial enrols 400 younger and 50 older patients, 450 in all. In and out of the trial, early surgery changes the risk of delirium by -0.05 in younger and -0.15 in older patients, so age is the effect modifier.
-
The trial's own risk difference
\[ \frac{400 \times (-0.05) + 50 \times (-0.15)}{450} = \frac{-20 - 7.5}{450} = \frac{-27.5}{450} = -0.061 \]
The trial is 88.9% younger, so its result sits near the younger value.
-
Weights for the whole source population
\[ w^S = \frac{1}{0.80} = 1.25, \quad w^S = \frac{1}{0.10} = 10 \]
The first weight is for a younger patient and the second for an older patient. Then 400 × 1.25 = 500 and 50 × 10 = 500, so the weighted trial rebuilds all 1,000 people.
-
Risk difference in the whole source population
\[ \frac{500 \times (-0.05) + 500 \times (-0.15)}{1000} = \frac{-25 - 75}{1000} = \frac{-100}{1000} = -0.100 \]
This is the generalised result.
-
Inverse-odds weights for the non-participants
\[ \frac{1 - 0.80}{0.80} = 0.25, \quad \frac{1 - 0.10}{0.10} = 9 \]
Again the first weight is for a younger patient and the second for an older patient. Then 400 × 0.25 = 100 and 50 × 9 = 450: the 550 people who stayed out, 18.2% younger and 81.8% older.
-
Risk difference in the non-participants
\[ \frac{100 \times (-0.05) + 450 \times (-0.15)}{550} = \frac{-5 - 67.5}{550} = \frac{-72.5}{550} = -0.132 \]
This is the transported result.
Result: One trial gives three results: -0.061 for its participants, -0.100 for the whole source population and -0.132 for the non-participants. Each is right for its own target; the two weighted results are right only because the effect modifier, age, is in the sampling model.
The simulated registry: carrying the trial result to two targets
In the simulated registry, the 950-patient trial measured postoperative delirium. The sampling model regresses trial entry on age, sex, frailty, dementia, anticoagulant use and ASA physical status 3 or higher in all 20,000 patients. By design, the effect of early surgery on delirium differs only by frailty on the odds ratio scale. On the risk difference scale reported here, it also changes with every risk factor for delirium, and all of those that drive trial entry are in the model.
One model gives two weights: $1/s(X)$ for the whole registry and the inverse odds for the 19,050 non-participants. Each risk difference comes from a weighted linear regression with a robust (sandwich) standard error, estimated from the observed residuals instead of the model's variance formula. A sandwich standard error repairs the standard error when the variance assumption is wrong. It does not repair a wrong mean model, and a biased coefficient keeps its bias.
Here the robust standard error also treats the weights as known, ignoring that $s(X)$ was itself estimated. The effective sample size (ESS), $(\sum w)^2 / \sum w^2$, is roughly the number of equally weighted patients that would give the same precision.
Delirium in the nested trial: risk difference by target
| Weight (target) | Risk difference, early minus later (95% CI) | Truth in the target | Robust SE | ESS (of 950) | Largest weight |
|---|---|---|---|---|---|
| None (the 950 participants) | -0.053 (-0.092 to -0.014) | -0.033 | 0.020 | 950 | 1 |
| $1/s(X)$ (whole registry, 20,000) | -0.053 (-0.221 to 0.114) | -0.038 | 0.085 | 199 | 695 |
| Inverse odds (non-participants, 19,050) | -0.052 (-0.226 to 0.122) | -0.038 | 0.089 | 184 | 694 |
What the table shows: a price, not a new answer
The three estimates are almost identical, and the three targets differ by less than one percentage point. With interval half-widths of 0.168 and 0.174 for the weighted estimates, against 0.039 for the trial alone, these data cannot tell the targets apart. What they show is the price of transport: a standard error 4.31 times the trial's for $1/s(X)$ and 4.46 times for the inverse odds.
The price comes from a few heavy weights. Older patients and patients with dementia rarely entered the trial, so the few who did stand for many registry patients, up to 695 each. Weights this large flag patients the trial almost never enrols, a near-violation of positivity for trial entry. The ESS falls to 199 and 184 of 950, about a fifth of the trial.
The Stata and R scripts below are the full simulation code for this registry, shared across this series and other articles. For this article, read only the sampling-score, weight and risk-difference lines. Lines that print CANON only write each result as plain text, and the VERIFY lines only check that the software is installed.
The simulated truth values come from the simulation's own settings file, which is not published. These truth values are fixed by how the data were generated, so treat them as the known answer the estimates are being checked against, not as results you can recompute from the printed output.
Stata: sampling-score model, weights and trial estimate
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
. * ---- 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)
Iteration 0: Log likelihood = -3821.7458
Iteration 1: Log likelihood = -3451.1611
Iteration 2: Log likelihood = -3354.8237
Iteration 3: Log likelihood = -3350.7977
Iteration 4: Log likelihood = -3350.7525
Iteration 5: Log likelihood = -3350.7525
Iteration 6: Log likelihood = -3350.7525
Logistic regression Number of obs = 20,000
LR chi2(6) = 941.99
Prob > chi2 = 0.0000
Log likelihood = -3350.7525 Pseudo R2 = 0.1232
------------------------------------------------------------------------------
trial | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
age_c | -.0474933 .0032768 -14.49 0.000 -.0539157 -.0410709
female | .0114397 .0754341 0.15 0.879 -.1364083 .1592878
frail | -.736586 .0957994 -7.69 0.000 -.9243494 -.5488225
dementia | -1.778538 .1827157 -9.73 0.000 -2.136654 -1.420422
asa3 | -.598664 .0729991 -8.20 0.000 -.7417396 -.4555884
anticoag | -.5789111 .1168301 -4.96 0.000 -.807894 -.3499283
_cons | -2.45805 .0781865 -31.44 0.000 -2.611292 -2.304807
------------------------------------------------------------------------------
. predict double s_hat, pr
. * target: the whole registry
. generate double w_ipsw = 1/s_hat if trial == 1
(19,050 missing values generated)
. * target: the non-participants (inverse odds of participation)
. generate double w_iosw = (1 - s_hat)/s_hat if trial == 1
(19,050 missing values generated)
. regress delirium surg24 if trial == 1, vce(robust)
Linear regression Number of obs = 950
F(1, 948) = 7.11
Prob > F = 0.0078
R-squared = 0.0075
Root MSE = .30333
------------------------------------------------------------------------------
| Robust
delirium | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.0528838 .0198341 -2.67 0.008 -.0918075 -.01396
_cons | .1304348 .0157191 8.30 0.000 .0995866 .161283
------------------------------------------------------------------------------R: transport and one-year death tables
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("Nested trial (n = 950): delirium risk difference moved to a target population, simulated data\n")
Nested trial (n = 950): delirium risk difference moved to a target population, simulated data
> print(round(t_tr, 4))
rd rd_lo rd_hi se se_ratio ess
trial -0.0529 -0.0918 -0.0140 0.0198 1.0000 950.0000
ipsw_whole_registry -0.0531 -0.2207 0.1144 0.0854 4.3054 199.0468
inverse_odds_non_participants -0.0518 -0.2256 0.1220 0.0886 4.4646 183.5473
max_weight target_rd
trial 1.0000 -0.0328
ipsw_whole_registry 695.0629 -0.0381
inverse_odds_non_participants 694.0629 -0.0384
> 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")
One-year death, observational part, simulated data
> print(round(t_cens, 4))
rd rd_lo rd_hi target_rd
naive_complete_case -0.1810 -0.1942 -0.1679 -0.1746
ipcw -0.1742 -0.1870 -0.1615 -0.1746
iptw_x_ipcw_ate -0.0352 -0.0601 -0.0102 -0.0378More than one kind at once: multiply the weights
When an analysis has more than one kind of missing person, the weights multiply:
$$w_i = w^A_i \times w^C_i \times w^S_i$$Here $w^A_i$, $w^C_i$ and $w^S_i$ are the treatment, censoring and selection weights. With $1/s(X)$ as the selection weight, their product is one over the probability of the whole path, from entering the study to staying in follow-up, when the models follow the order of events; the inverse odds replace that factor when the target is the non-participants. So each model conditions on the steps before it: the treatment model is fitted among those selected into the analysis, and the censoring model includes the treatment received as a predictor.
Not every factor is always needed. In a 1:1 randomised trial, $w^A_i$ is the same for everyone, so transporting the trial needs only $w^S_i$. The registry's death analysis has no selection step and uses $w^A_i \times w^C_i$, the last row of the death table.
What each weight assumes
The weights rest on the assumptions below; exchangeability, positivity and a correct model apply to each mechanism separately [6, 8].
- Consistency and no interference. Each patient's observed outcome is the outcome under the treatment actually received, and early surgery means the same intervention in the trial, in the registry and in any target population. One patient's treatment does not change another patient's outcome.
- Exchangeability. Treatment: no unmeasured confounding given $X$. Censoring: given $X$ and treatment, patients lost at any time would have had the same risk of dying from then on as similar patients still followed at that time. Selection: given $X$, participants and target share the same treatment effect, so every effect modifier, on the scale of the reported effect, whose distribution differs between them must be in $X$.
- Positivity. Every covariate pattern needs a non-zero chance of each treatment, of staying in follow-up and of entering the trial.
- Correct models. A wrong model gives wrong weights, so after weighting check that each covariate has the same distribution in the weighted group as in its comparison group, and in the weighted trial as in the target.
The selection assumption is easy to skip. A sampling model built from what predicts enrolment, rather than from what modifies the effect, can balance the wrong variables.
The IPW family side by side
| Weight | Missing people | Probability modelled | Weight of an observed person | Target | Key assumption |
|---|---|---|---|---|---|
| IPTW | Same patients, other treatment | $e(X) = P(A = 1 \mid X)$ | $1/e(X)$ or $1/(1 - e(X))$ | Population analysed (ATE) | No unmeasured confounding |
| IPCW | Patients lost to follow-up | $P(C_i(t) = 0 \mid X_i, A_i)$ | One over that probability | Same population, nobody lost | Loss non-informative given $X$ and $A$ |
| IPSW | Source population not enrolled | $s(X) = P(S = 1 \mid X)$ | $1/s(X)$ | Whole source population | Effect modifiers in $X$ |
| Inverse odds | Non-participants | $s(X)$ | $(1 - s(X))/s(X)$ | Non-participants (or an external population) | Effect modifiers in $X$ |
Common misreadings and their fixes
-
"Weighting by the sampling score makes the trial apply to everyone."
The weights rebalance only what the sampling model contains, toward the one target they were built for.
Fix: Weighting by the inverse sampling score makes the trial resemble the target population only on the covariates in the sampling model, and only when the effect modifiers are among them and every covariate pattern in the target has a chance of being in the trial. Name the target population before choosing the weight.
-
"Censoring is a confounder, so adjust for it."
Censoring happens after treatment and hides outcomes; it does not drive the treatment choice.
Fix: Weight observed patients by one over their modelled probability of staying in follow-up, given the measured predictors of loss and the treatment received.
-
"IPCW makes the registry comparison causal."
In the simulated registry, IPCW alone gives -0.174, near the censoring-free association of -0.175 and far from the causal -0.038.
Fix: Multiply the censoring weights by treatment weights; in the registry this gives -0.035.
-
"A transported estimate is as precise as the trial."
Transport to the whole registry multiplied the standard error by 4.31 and left an ESS of 199 of 950.
Fix: Report the ESS, the largest weight and the interval beside the trial's own estimate.
What to do in your own analysis
- Name the missing people and the target population before fitting any model.
- Fit one model per mechanism: treatment on the confounders, follow-up on the predictors of loss and treatment, trial entry on every suspected effect modifier.
- For each weight, inspect the smallest fitted probabilities, the largest weight and the ESS.
- Consider a robust (sandwich) standard error, or a bootstrap (repeat the whole analysis on resampled patients, refitting every weight model each time).
- Report each weighted estimate beside its unweighted one.
Glossary
- IPCW (การถ่วงน้ำหนักด้วยส่วนกลับของความน่าจะเป็นที่ไม่ถูกเซ็นเซอร์)
- Inverse probability of censoring weighting: followed patients weighted by one over their probability of staying in follow-up.
- IPSW (การถ่วงน้ำหนักด้วยส่วนกลับของความน่าจะเป็นในการถูกคัดเลือก)
- Inverse probability of selection weighting: trial participants weighted by one over the sampling score, which carries the result to the whole source population.
- sampling score (คะแนนการถูกคัดเลือกเข้าการศึกษา)
- s(X) = P(S = 1 given X), the probability of being in the study given covariates.
- inverse odds of sampling weights
- Weights (1 - s(X))/s(X) that carry a trial result to non-participants.
- informative censoring
- Loss to follow-up whose chance depends on the outcome risk, directly or through characteristics such as dementia.
- generalisability
- Whether a result holds in the population the study was drawn from.
- transportability
- Whether a result holds in a population the study did not sample.
References
- Seaman SR, White IR. Review of inverse probability weighting for dealing with missing data. Stat Methods Med Res. 2013;22(3):278-295. doi:10.1177/0962280210395740 https://doi.org/10.1177/0962280210395740
- Robins JM, Finkelstein DM. Correcting for noncompliance and dependent censoring in an AIDS clinical trial with inverse probability of censoring weighted (IPCW) log-rank tests. Biometrics. 2000;56:779-788. doi:10.1111/j.0006-341x.2000.00779.x https://doi.org/10.1111/j.0006-341x.2000.00779.x
- Hernán MA, Hernandez-Diaz S, Robins JM. A structural approach to selection bias. Epidemiology. 2004;15:615-625. doi:10.1097/01.ede.0000135174.63482.43 https://doi.org/10.1097/01.ede.0000135174.63482.43
- Cole SR, Stuart EA. Generalizing evidence from randomized clinical trials to target populations: the ACTG 320 trial. Am J Epidemiol. 2010;172:107-115. doi:10.1093/aje/kwq084 https://doi.org/10.1093/aje/kwq084
- Stuart EA, Cole SR, Bradshaw CP, Leaf PJ. The use of propensity scores to assess the generalizability of results from randomized trials. J R Stat Soc Ser A Stat Soc. 2011;174(2):369-386. doi:10.1111/j.1467-985x.2010.00673.x https://doi.org/10.1111/j.1467-985x.2010.00673.x
- Dahabreh IJ, Robertson SE, Tchetgen EJ, Stuart EA, Hernán MA. Generalizing causal inferences from individuals in randomized trials to all trial-eligible individuals. Biometrics. 2019;75(2):685-694. doi:10.1111/biom.13009 https://doi.org/10.1111/biom.13009
- Westreich D, Edwards JK, Lesko CR, Stuart E, Cole SR. Transportability of trial results using inverse odds of sampling weights. Am J Epidemiol. 2017;186:1010-1014. doi:10.1093/aje/kwx164 https://doi.org/10.1093/aje/kwx164
- Lesko CR, Buchanan AL, Westreich D, Edwards JK, Hudgens MG, Cole SR. Generalizing study results: a potential outcomes perspective. Epidemiology. 2017;28(4):553-561. doi:10.1097/ede.0000000000000664 https://doi.org/10.1097/ede.0000000000000664
Key takeaways
- Inverse probability weighting stands in for unseen people by up-weighting observed people who were unlikely to be observed.
- Treatment, censoring and selection weights share one template but need their own models, assumptions and named targets.
- Censoring weights remove bias from loss to follow-up, including informative loss explained by the measured predictors of loss, but not confounding: in the simulated registry, IPCW alone gave -0.174 against a causal -0.038.
- When the sampling model holds the effect modifiers, one over the sampling score generalises a trial to its whole source population and the inverse odds transport it to non-participants.
- A selection weight transports an effect only through the effect modifiers in its model, at a precision cost the effective sample size makes visible.
Related in the wiki: [[iptw-guide]]
0