Marginal Structural Models: เมื่อตัวกวนวันนี้คือผลของการรักษาเมื่อวาน

* Simulated ICU cohort, three treatment days: treatment-confounder feedback.
* Simulated data: not evidence about any real drug or patient.
* Question: what does each treated day of drug X do to the risk of death by day 28?
* Methods: naive regressions (with and without L_t) versus a marginal structural model (MSM)
* fitted with cumulative stabilised inverse probability of treatment weights.
* Variables: L_t = 1 when severe organ dysfunction is present on day t, A_t = 1 when the drug is given on day t,
* V = 1 for age 65 or older, Y = 1 for death by day 28.
version 18
clear all
set more off

* Check the built-in commands this file uses.
capture which logit
display "VERIFY logit " cond(_rc == 0, "available", "missing")
capture which glm
display "VERIFY glm " cond(_rc == 0, "available", "missing")

* Read the simulated cohort file (not published): one row per patient-day (days 0, 2, 4).
import delimited using ../../datasets/W2/W2.csv, clear case(preserve) varnames(1)
sort id t
count if t == 0
display "CANON w2.n " strtrim(string(r(N), "%12.0f"))
count
display "CANON w2.n_rows " strtrim(string(r(N), "%12.0f"))
count if t == 2 & Y == 1
display "CANON w2.deaths " strtrim(string(r(N), "%12.0f"))
count if A == 1
display "CANON w2.treated_days " strtrim(string(r(N), "%12.0f"))

* Naive 1: pool all patient-days and regress death on that day's treatment (linear risk model), CI clustered by patient.
glm Y A, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.naive.unadj.b " strtrim(string(_b[A], "%12.4f"))
display "CANON w2.naive.unadj.se " strtrim(string(_se[A], "%12.4f"))
display "CANON w2.naive.unadj.lo " strtrim(string(_b[A] - invnormal(0.975) * _se[A], "%12.4f"))
display "CANON w2.naive.unadj.hi " strtrim(string(_b[A] + invnormal(0.975) * _se[A], "%12.4f"))

* Naive 2: the same pooled regression adjusted for that day's organ dysfunction L_t and age V.
glm Y A L V, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.naive.adj.b " strtrim(string(_b[A], "%12.4f"))
display "CANON w2.naive.adj.se " strtrim(string(_se[A], "%12.4f"))
display "CANON w2.naive.adj.lo " strtrim(string(_b[A] - invnormal(0.975) * _se[A], "%12.4f"))
display "CANON w2.naive.adj.hi " strtrim(string(_b[A] + invnormal(0.975) * _se[A], "%12.4f"))

* Cumulative treatment through each day, and organ dysfunction on days 0, 2 and 4 copied to every row.
by id (t): generate cumA = sum(A)
by id (t): generate L0 = L[1]
by id (t): generate L1 = L[2]
by id (t): generate L2 = L[3]
count if t == 2 & cumA == 3
display "CANON w2.always_treated " strtrim(string(r(N), "%12.0f"))
count if t == 2 & cumA == 0
display "CANON w2.never_treated " strtrim(string(r(N), "%12.0f"))

* Naive 3: one row per patient (day 4), death on the number of treated days (no adjustment).
glm Y cumA if t == 2, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.cum.crude.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.cum.crude.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.cum.crude.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.cum.crude.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))

* Naive 4: the same, adjusted for age and organ dysfunction on days 0, 2 and 4.
glm Y cumA V L0 L1 L2 if t == 2, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.cum.adj.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.cum.adj.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.cum.adj.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.cum.adj.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))

* Treatment models, pooled over days: the denominator uses history (L_t, A_{t-1}, V).
logit A L A_prev V i.t
predict double pden, pr
* Numerator for the stabilised weight: same model without L_t.
logit A A_prev V i.t
predict double pnum, pr
* Numerator for the marginal MSM: without L_t and without V.
logit A A_prev i.t
predict double pnumm, pr

* Probability of the treatment each patient actually received on each day.
generate double fden = cond(A == 1, pden, 1 - pden)
generate double fnum = cond(A == 1, pnum, 1 - pnum)
generate double fnumm = cond(A == 1, pnumm, 1 - pnumm)

* Cumulative product over days within each patient (as a running sum of logs) gives the weight at day 4.
by id (t): generate double w = exp(sum(ln(1 / fden)))
by id (t): generate double sw = exp(sum(ln(fnum / fden)))
by id (t): generate double swm = exp(sum(ln(fnumm / fden)))

* Weight summaries at day 4; the effective sample size is (sum of weights)^2 / sum of squared weights.
summarize sw if t == 2, detail
display "CANON w2.sw.mean " strtrim(string(r(mean), "%12.4f"))
display "CANON w2.sw.sd " strtrim(string(r(sd), "%12.4f"))
display "CANON w2.sw.min " strtrim(string(r(min), "%12.4f"))
display "CANON w2.sw.p1 " strtrim(string(r(p1), "%12.4f"))
display "CANON w2.sw.p99 " strtrim(string(r(p99), "%12.4f"))
display "CANON w2.sw.max " strtrim(string(r(max), "%12.4f"))
generate double sw_sq = sw^2
summarize sw if t == 2
scalar sum_sw = r(sum)
summarize sw_sq if t == 2
display "CANON w2.sw.ess " strtrim(string(sum_sw^2 / r(sum), "%12.4f"))
generate double w_sq = w^2
summarize w if t == 2
display "CANON w2.w.mean " strtrim(string(r(mean), "%12.4f"))
display "CANON w2.w.max " strtrim(string(r(max), "%12.4f"))
scalar sum_w = r(sum)
summarize w_sq if t == 2
display "CANON w2.w.ess " strtrim(string(sum_w^2 / r(sum), "%12.4f"))
summarize swm if t == 2
display "CANON w2.swm.mean " strtrim(string(r(mean), "%12.4f"))
display "CANON w2.swm.max " strtrim(string(r(max), "%12.4f"))

* MSM: weighted regression of death on cumulative treatment and age, CI clustered by patient.
glm Y cumA V if t == 2 [pweight = sw], family(gaussian) link(identity) vce(cluster id)
display "CANON w2.msm.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.msm.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.msm.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.msm.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))
* Always treated versus never treated is three treated days: three times the per-day effect.
display "CANON w2.msm.rd_always_never " strtrim(string(3 * _b[cumA], "%12.4f"))
display "CANON w2.msm.rd_always_never.lo " strtrim(string(3 * (_b[cumA] - invnormal(0.975) * _se[cumA]), "%12.4f"))
display "CANON w2.msm.rd_always_never.hi " strtrim(string(3 * (_b[cumA] + invnormal(0.975) * _se[cumA]), "%12.4f"))

* Marginal MSM (numerator without V, no V in the model): gives the risks under never and always treated.
glm Y cumA if t == 2 [pweight = swm], family(gaussian) link(identity) vce(cluster id)
display "CANON w2.msm_marg.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.msm_marg.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.msm_marg.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.msm_marg.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.msm_marg.risk_never " strtrim(string(_b[_cons], "%12.4f"))
display "CANON w2.msm_marg.risk_always " strtrim(string(_b[_cons] + 3 * _b[cumA], "%12.4f"))

display "Simulated data (ข้อมูลจำลอง): no number above is evidence about any real drug."

* ---- The true values built into the simulated cohort (exact: computed from the design, no sampling) ----
* The cohort was simulated from these probabilities. U is an unmeasured patient trait; L_6 is the day-6
* value of L (after the last treatment day, never recorded).
*   P(U = 1) = pU and P(V = 1) = pV
*   P(L_t = 1) = lBase + lU*U + lA*A_{t-1}, with A_{t-1} = 0 on day 0
*   logit P(A_t = 1) = aBase + aL*L_t + aA*A_{t-1} + aV*V
*   P(death by day 28) = yBase + yU*U + yV*V + yL*(L_0 + L_1 + L_2 + L_6) + yA*(A_0 + A_1 + A_2)
* The true values need no patient data, so start from an empty dataset.
clear
scalar pU = 0.40
scalar pV = 0.45
scalar lBase = 0.32
scalar lU = 0.45
scalar lA = -0.30
scalar aBase = -1.2
scalar aL = 2.0
scalar aA = 1.5
scalar aV = -0.5
scalar yBase = 0.05
scalar yU = 0.25
scalar yV = 0.05
scalar yL = 0.15
scalar yA = -0.005
display "CANON w2.truth.design.p_u " strtrim(string(pU, "%12.4f"))
display "CANON w2.truth.design.p_v " strtrim(string(pV, "%12.4f"))
display "CANON w2.truth.design.l_base " strtrim(string(lBase, "%12.4f"))
display "CANON w2.truth.design.l_u " strtrim(string(lU, "%12.4f"))
display "CANON w2.truth.design.l_a " strtrim(string(lA, "%12.4f"))
display "CANON w2.truth.design.a_base " strtrim(string(aBase, "%12.4f"))
display "CANON w2.truth.design.a_l " strtrim(string(aL, "%12.4f"))
display "CANON w2.truth.design.a_a " strtrim(string(aA, "%12.4f"))
display "CANON w2.truth.design.a_v " strtrim(string(aV, "%12.4f"))
display "CANON w2.truth.design.y_base " strtrim(string(yBase, "%12.4f"))
display "CANON w2.truth.design.y_u " strtrim(string(yU, "%12.4f"))
display "CANON w2.truth.design.y_v " strtrim(string(yV, "%12.4f"))
display "CANON w2.truth.design.y_l " strtrim(string(yL, "%12.4f"))
display "CANON w2.truth.design.y_a " strtrim(string(yA, "%12.4f"))

* Each treated day lowers the risk of death directly (yA) and by lowering the next day's L (yL x lA).
scalar g1 = yA + yL * lA
display "CANON w2.truth.msm.g1_per_treated_day " strtrim(string(g1, "%12.4f"))
display "CANON w2.truth.msm.g1_direct_part " strtrim(string(yA, "%12.4f"))
display "CANON w2.truth.msm.g1_through_l_part " strtrim(string(yL * lA, "%12.4f"))
* Risk of death by day 28 if every patient were treated on k of the three days (k = 0 never, k = 3 always).
* EL is P(L_t = 1) after an untreated day; L_0, L_1, L_2 and L_6 each add yL x EL, and each treated day adds g1.
scalar EL = lBase + lU * pU
scalar rNever = yBase + yU * pU + yV * pV + yL * 4 * EL
display "CANON w2.truth.msm.risk_never_treated " strtrim(string(rNever, "%12.4f"))
display "CANON w2.truth.msm.risk_cum1 " strtrim(string(rNever + g1, "%12.4f"))
display "CANON w2.truth.msm.risk_cum2 " strtrim(string(rNever + 2 * g1, "%12.4f"))
display "CANON w2.truth.msm.risk_always_treated " strtrim(string(rNever + 3 * g1, "%12.4f"))
display "CANON w2.truth.msm.rd_always_vs_never " strtrim(string(3 * g1, "%12.4f"))
* The model with age V (the one fitted above): intercept at V = 0 and the coefficient of V.
display "CANON w2.truth.msm.g0_conditional " strtrim(string(yBase + yU * pU + yL * 4 * EL, "%12.4f"))
display "CANON w2.truth.msm.g2_v " strtrim(string(yV, "%12.4f"))

* Where each naive regression lands in an infinitely large cohort: list all 512 combinations of
* U, V, L_0, A_0, L_1, A_1, L_2, A_2 and L_6 with their probabilities under the observed treatment
* process, then fit the same regressions to the expected risk, weighting each combination by its probability.
set obs 512
generate int cell = _n - 1
generate byte U = mod(cell, 2)
generate byte V = mod(floor(cell / 2), 2)
generate byte L0 = mod(floor(cell / 4), 2)
generate byte A0 = mod(floor(cell / 8), 2)
generate byte L1 = mod(floor(cell / 16), 2)
generate byte A1 = mod(floor(cell / 32), 2)
generate byte L2 = mod(floor(cell / 64), 2)
generate byte A2 = mod(floor(cell / 128), 2)
generate byte L6 = mod(floor(cell / 256), 2)
generate double pL0 = lBase + lU * U
generate double pA0 = invlogit(aBase + aL * L0 + aV * V)
generate double pL1 = lBase + lU * U + lA * A0
generate double pA1 = invlogit(aBase + aL * L1 + aA * A0 + aV * V)
generate double pL2 = lBase + lU * U + lA * A1
generate double pA2 = invlogit(aBase + aL * L2 + aA * A1 + aV * V)
generate double pL6 = lBase + lU * U + lA * A2
generate double p = cond(U == 1, pU, 1 - pU) * cond(V == 1, pV, 1 - pV)
replace p = p * cond(L0 == 1, pL0, 1 - pL0) * cond(A0 == 1, pA0, 1 - pA0)
replace p = p * cond(L1 == 1, pL1, 1 - pL1) * cond(A1 == 1, pA1, 1 - pA1)
replace p = p * cond(L2 == 1, pL2, 1 - pL2) * cond(A2 == 1, pA2, 1 - pA2)
replace p = p * cond(L6 == 1, pL6, 1 - pL6)
generate double risk = yBase + yU * U + yV * V + yL * (L0 + L1 + L2 + L6) + yA * (A0 + A1 + A2)
generate cumA = A0 + A1 + A2
* Naive 3 and 4 (one row per patient).
quietly regress risk cumA [aweight = p]
display "CANON w2.truth.naive_limit.cum_crude " strtrim(string(_b[cumA], "%12.4f"))
quietly regress risk cumA V L0 L1 L2 [aweight = p]
display "CANON w2.truth.naive_limit.cum_adj " strtrim(string(_b[cumA], "%12.4f"))
* Naive 1 and 2 (one row per patient-day): three rows per combination, one per treatment day.
expand 3
bysort cell: generate t = _n - 1
generate A = cond(t == 0, A0, cond(t == 1, A1, A2))
generate L = cond(t == 0, L0, cond(t == 1, L1, L2))
quietly regress risk A [aweight = p]
display "CANON w2.truth.naive_limit.unadj " strtrim(string(_b[A], "%12.4f"))
quietly regress risk A L V [aweight = p]
display "CANON w2.truth.naive_limit.adj " strtrim(string(_b[A], "%12.4f"))
display "These true values describe the simulation design (simulated data), not any real drug."

* ---- Follow-up length, the marginal MSM intercept, and how many days each patient was treated ----
* Death is counted up to day 28.
scalar deathDay = 28
display "CANON w2.truth.design.death_day " strtrim(string(deathDay, "%12.0f"))
* The intercept of the marginal MSM is the risk of death by day 28 if no patient were ever treated.
display "CANON w2.truth.msm.g0_marginal " strtrim(string(rNever, "%12.4f"))
* Read the simulated cohort file (not published) again: share who died, patients treated on one or two days.
import delimited using ../../datasets/W2/W2.csv, clear case(preserve) varnames(1)
sort id t
by id (t): generate cumA = sum(A)
summarize Y if t == 2
display "CANON w2.death_risk " strtrim(string(r(mean), "%12.4f"))
count if t == 2 & cumA == 1
display "CANON w2.treated_one_day " strtrim(string(r(N), "%12.0f"))
count if t == 2 & cumA == 2
display "CANON w2.treated_two_days " strtrim(string(r(N), "%12.0f"))
tabulate cumA if t == 2
display "Simulated data (ข้อมูลจำลอง): not evidence about any real drug or patient."

. * MSM: weighted regression of death on cumulative treatment and age, CI clust
> ered by patient.
. glm Y cumA V if t == 2 [pweight = sw], family(gaussian) link(identity) vce(cl
> uster id)

Iteration 0:  Log pseudolikelihood = -2096.0075

Generalized linear models                         Number of obs   =      3,000
Optimization     : ML                             Residual df     =      2,997
                                                  Scale parameter =   .2370902
Deviance         =  710.5593289                   (1/df) Deviance =   .2370902
Pearson          =  710.5593289                   (1/df) Pearson  =   .2370902

Variance function: V(u) = 1                       [Gaussian]
Link function    : g(u) = u                       [Identity]

                                                  AIC             =   1.399338
Log pseudolikelihood = -2096.007537               BIC             =  -23284.52

                                 (Std. err. adjusted for 3,000 clusters in id)
------------------------------------------------------------------------------
             |               Robust
           Y | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        cumA |  -.0572908   .0118918    -4.82   0.000    -.0805982   -.0339834
           V |   .0477372   .0242604     1.97   0.049     .0001876    .0952868
       _cons |   .4654255   .0266901    17.44   0.000     .4131138    .5177371
------------------------------------------------------------------------------

# Simulated ICU cohort, three treatment days: treatment-confounder feedback.
# Simulated data: not evidence about any real drug or patient.
# Question: what does each treated day of drug X do to the risk of death by day 28?
# Methods: naive regressions (with and without L_t) versus a marginal structural model (MSM)
# fitted with cumulative stabilised inverse probability of treatment weights.
# Variables: L_t = 1 when severe organ dysfunction is present on day t, A_t = 1 when the drug is given on day t,
# V = 1 for age 65 or older, Y = 1 for death by day 28.

# Check the packages this script uses; install into the user library when one is missing.
need <- c("sandwich", "ipw", "WeightIt")
for (p in need) {
  if (!requireNamespace(p, quietly = TRUE)) {
    try(install.packages(p, repos = "https://cloud.r-project.org"), silent = TRUE)
  }
  cat(sprintf("VERIFY %s %s\n", p, if (requireNamespace(p, quietly = TRUE)) "available" else "missing"))
}

# Print each result on its own labelled line: 4 decimals for estimates, integers for counts.
canon <- function(key, value, digits = 4) {
  cat(sprintf("CANON w2.%s %s\n", key, formatC(value, format = "f", digits = digits)))
}

# Read the simulated cohort file (not published): one row per patient-day (days 0, 2, 4).
d <- read.csv(file.path("..", "..", "datasets", "W2", "W2.csv"))
d <- d[order(d$id, d$t), ]
canon("n", length(unique(d$id)), 0)
canon("n_rows", nrow(d), 0)
canon("deaths", sum(d$Y[d$t == 2]), 0)
canon("treated_days", sum(d$A), 0)

# Cluster-robust (sandwich) 95% CI by patient for one coefficient; HC0 with the G/(G-1) factor, as Stata glm does.
ci_cl <- function(fit, term, cluster) {
  V <- sandwich::vcovCL(fit, cluster = cluster, type = "HC0", cadjust = TRUE)
  b <- unname(coef(fit)[term]); se <- sqrt(V[term, term])
  c(b = b, se = se, lo = b - qnorm(0.975) * se, hi = b + qnorm(0.975) * se)
}
print_est <- function(key, est) {
  canon(paste0(key, ".b"), est[["b"]]); canon(paste0(key, ".se"), est[["se"]])
  canon(paste0(key, ".lo"), est[["lo"]]); canon(paste0(key, ".hi"), est[["hi"]])
}

# Naive 1: pool all patient-days and regress death on that day's treatment (linear risk model).
f_unadj <- glm(Y ~ A, data = d, family = gaussian)
print_est("naive.unadj", ci_cl(f_unadj, "A", ~id))

# Naive 2: the same pooled regression adjusted for that day's organ dysfunction L_t and age V.
f_adj <- glm(Y ~ A + L + V, data = d, family = gaussian)
print_est("naive.adj", ci_cl(f_adj, "A", ~id))

# Cumulative treatment through each day, and organ dysfunction on each day, in wide form.
d$cumA <- ave(d$A, d$id, FUN = cumsum)
wide <- d[d$t == 2, c("id", "V", "Y", "cumA")]
Lw <- reshape(d[, c("id", "t", "L")], idvar = "id", timevar = "t", direction = "wide")
wide <- merge(wide, Lw, by = "id")
wide <- wide[order(wide$id), ]
canon("always_treated", sum(wide$cumA == 3), 0)
canon("never_treated", sum(wide$cumA == 0), 0)

# Naive 3: one row per patient, death on the number of treated days (no adjustment).
f_cc <- glm(Y ~ cumA, data = wide, family = gaussian)
print_est("cum.crude", ci_cl(f_cc, "cumA", ~id))

# Naive 4: the same, adjusted for age and organ dysfunction on days 0, 2 and 4.
f_ca <- glm(Y ~ cumA + V + L.0 + L.1 + L.2, data = wide, family = gaussian)
print_est("cum.adj", ci_cl(f_ca, "cumA", ~id))

# Treatment models, pooled over days: denominator uses history (L_t, A_{t-1}, V); numerators drop L_t.
den <- glm(A ~ L + A_prev + V + factor(t), data = d, family = binomial)
num <- glm(A ~ A_prev + V + factor(t), data = d, family = binomial)
numm <- glm(A ~ A_prev + factor(t), data = d, family = binomial)
pa <- function(fit) { p <- fitted(fit); ifelse(d$A == 1, p, 1 - p) }

# Cumulative product over days within each patient gives the weight at day 4.
d$w <- ave(1 / pa(den), d$id, FUN = cumprod)
d$sw <- ave(pa(num) / pa(den), d$id, FUN = cumprod)
d$swm <- ave(pa(numm) / pa(den), d$id, FUN = cumprod)
last <- d[d$t == 2, ]
last <- last[order(last$id), ]

# Weight summaries at day 4 (percentiles use the same averaging rule as Stata summarize, detail).
ess <- function(x) sum(x)^2 / sum(x^2)
q <- function(x, p) unname(quantile(x, p, type = 2))
canon("sw.mean", mean(last$sw)); canon("sw.sd", sd(last$sw))
canon("sw.min", min(last$sw)); canon("sw.p1", q(last$sw, 0.01))
canon("sw.p99", q(last$sw, 0.99)); canon("sw.max", max(last$sw))
canon("sw.ess", ess(last$sw))
canon("w.mean", mean(last$w)); canon("w.max", max(last$w)); canon("w.ess", ess(last$w))
canon("swm.mean", mean(last$swm)); canon("swm.max", max(last$swm))

# MSM: weighted regression of death on cumulative treatment and age, cluster-robust CI by patient.
msm <- glm(Y ~ cumA + V, data = last, weights = sw, family = gaussian)
e_msm <- ci_cl(msm, "cumA", ~id)
print_est("msm", e_msm)
canon("msm.rd_always_never", 3 * e_msm[["b"]])
canon("msm.rd_always_never.lo", 3 * e_msm[["lo"]])
canon("msm.rd_always_never.hi", 3 * e_msm[["hi"]])

# Marginal MSM (numerator without V, no V in the model): gives the risks under never and always treated.
msmm <- glm(Y ~ cumA, data = last, weights = swm, family = gaussian)
e_m <- ci_cl(msmm, "cumA", ~id)
print_est("msm_marg", e_m)
canon("msm_marg.risk_never", unname(coef(msmm)[1]))
canon("msm_marg.risk_always", unname(coef(msmm)[1] + 3 * coef(msmm)[2]))

# Cross-check: ipw::ipwtm builds the same pooled stabilised weights from its own code.
if (requireNamespace("ipw", quietly = TRUE)) {
  dd <- d
  dd$tf <- factor(dd$t)
  iw <- ipw::ipwtm(exposure = A, family = "binomial", link = "logit",
                   numerator = ~ A_prev + V + tf, denominator = ~ L + A_prev + V + tf,
                   id = id, timevar = t, type = "all", data = dd)
  cat(sprintf("CHECK ipwtm versus hand-built weights, largest absolute difference %.2e\n",
              max(abs(iw$ipw.weights - dd$sw))))
}

# Cross-check: WeightIt::weightitMSM fits one treatment model per day (not pooled), so it differs slightly.
if (requireNamespace("WeightIt", quietly = TRUE)) {
  wd <- reshape(d[, c("id", "t", "V", "L", "A", "Y")], idvar = c("id", "V", "Y"),
                timevar = "t", direction = "wide")
  wd <- wd[order(wd$id), ]
  wm <- WeightIt::weightitMSM(list(A.0 ~ L.0 + V, A.1 ~ A.0 + L.1 + V, A.2 ~ A.1 + L.2 + V),
                              data = wd, method = "glm", stabilize = TRUE)
  wd$cumA <- wd$A.0 + wd$A.1 + wd$A.2
  fw <- glm(Y ~ cumA + V, data = wd, weights = wm$weights, family = gaussian)
  ew <- ci_cl(fw, "cumA", ~id)
  canon("msm_weightit.b.r", ew[["b"]]); canon("msm_weightit.lo.r", ew[["lo"]]); canon("msm_weightit.hi.r", ew[["hi"]])
}

cat("Simulated data (ข้อมูลจำลอง): no number above is evidence about any real drug.\n")

# ---- The true values built into the simulated cohort (exact: computed from the design, no sampling) ----
# The cohort was simulated from these probabilities. U is an unmeasured patient trait; L_6 is the day-6
# value of L (after the last treatment day, never recorded).
#   P(U = 1) = pU and P(V = 1) = pV
#   P(L_t = 1) = lBase + lU*U + lA*A_{t-1}, with A_{t-1} = 0 on day 0
#   logit P(A_t = 1) = aBase + aL*L_t + aA*A_{t-1} + aV*V
#   P(death by day 28) = yBase + yU*U + yV*V + yL*(L_0 + L_1 + L_2 + L_6) + yA*(A_0 + A_1 + A_2)
pU <- 0.40; pV <- 0.45
lBase <- 0.32; lU <- 0.45; lA <- -0.30
aBase <- -1.2; aL <- 2.0; aA <- 1.5; aV <- -0.5
yBase <- 0.05; yU <- 0.25; yV <- 0.05; yL <- 0.15; yA <- -0.005
design <- c(p_u = pU, p_v = pV, l_base = lBase, l_u = lU, l_a = lA, a_base = aBase, a_l = aL, a_a = aA,
            a_v = aV, y_base = yBase, y_u = yU, y_v = yV, y_l = yL, y_a = yA)
for (k in names(design)) canon(paste0("truth.design.", k), design[[k]])

# Each treated day lowers the risk of death directly (yA) and by lowering the next day's L (yL x lA).
g1 <- yA + yL * lA
canon("truth.msm.g1_per_treated_day", g1)
canon("truth.msm.g1_direct_part", yA)
canon("truth.msm.g1_through_l_part", yL * lA)
# Risk of death by day 28 if every patient were treated on k of the three days (k = 0 never, k = 3 always).
# EL is P(L_t = 1) after an untreated day; L_0, L_1, L_2 and L_6 each add yL x EL, and each treated day adds g1.
EL <- lBase + lU * pU
r_never <- yBase + yU * pU + yV * pV + yL * 4 * EL
canon("truth.msm.risk_never_treated", r_never)
canon("truth.msm.risk_cum1", r_never + g1)
canon("truth.msm.risk_cum2", r_never + 2 * g1)
canon("truth.msm.risk_always_treated", r_never + 3 * g1)
canon("truth.msm.rd_always_vs_never", 3 * g1)
# The model with age V (the one fitted above): intercept at V = 0 and the coefficient of V.
canon("truth.msm.g0_conditional", yBase + yU * pU + yL * 4 * EL)
canon("truth.msm.g2_v", yV)

# Where each naive regression lands in an infinitely large cohort: list all 512 combinations of
# U, V, L_0, A_0, L_1, A_1, L_2, A_2 and L_6 with their probabilities under the observed treatment
# process, then fit the same regressions to the expected risk, weighting each combination by its probability.
g <- expand.grid(U = 0:1, V = 0:1, L0 = 0:1, A0 = 0:1, L1 = 0:1, A1 = 0:1, L2 = 0:1, A2 = 0:1, L6 = 0:1)
probL <- function(U, Aprev) lBase + lU * U + lA * Aprev
probA <- function(L, Aprev, V) plogis(aBase + aL * L + aA * Aprev + aV * V)
pick <- function(x, p) ifelse(x == 1, p, 1 - p)
g$p <- with(g, pick(U, pU) * pick(V, pV) * pick(L0, probL(U, 0)) * pick(A0, probA(L0, 0, V)) *
              pick(L1, probL(U, A0)) * pick(A1, probA(L1, A0, V)) *
              pick(L2, probL(U, A1)) * pick(A2, probA(L2, A1, V)) * pick(L6, probL(U, A2)))
g$risk <- with(g, yBase + yU * U + yV * V + yL * (L0 + L1 + L2 + L6) + yA * (A0 + A1 + A2))
g$cumA <- g$A0 + g$A1 + g$A2
# Naive 3 and 4 (one row per patient).
lim_cum <- coef(lm(risk ~ cumA, data = g, weights = p))[["cumA"]]
lim_cum_adj <- coef(lm(risk ~ cumA + V + L0 + L1 + L2, data = g, weights = p))[["cumA"]]
canon("truth.naive_limit.cum_crude", lim_cum)
canon("truth.naive_limit.cum_adj", lim_cum_adj)
# Naive 1 and 2 (one row per patient-day): three rows per combination, one per treatment day.
gd <- rbind(transform(g, A = A0, L = L0), transform(g, A = A1, L = L1), transform(g, A = A2, L = L2))
lim_unadj <- coef(lm(risk ~ A, data = gd, weights = p))[["A"]]
lim_adj <- coef(lm(risk ~ A + L + V, data = gd, weights = p))[["A"]]
canon("truth.naive_limit.unadj", lim_unadj)
canon("truth.naive_limit.adj", lim_adj)

# ---- summary tables (simulated data) ----
# Effect of one treated day on the risk of death by day 28: each estimate with its 95% CI (clustered by
# patient) beside the value it converges to in an infinitely large cohort (for the MSM: the true effect).
est3 <- function(e) unname(e[c("b", "lo", "hi")])
t_eff <- rbind(
  naive_pooled = c(est3(ci_cl(f_unadj, "A", ~id)), lim_unadj),
  naive_pooled_adj_L_V = c(est3(ci_cl(f_adj, "A", ~id)), lim_adj),
  naive_cumulative = c(est3(ci_cl(f_cc, "cumA", ~id)), lim_cum),
  naive_cumulative_adj = c(est3(ci_cl(f_ca, "cumA", ~id)), lim_cum_adj),
  msm_stabilised_weights = c(est3(e_msm), g1),
  msm_marginal = c(est3(e_m), g1))
colnames(t_eff) <- c("per_day", "lo", "hi", "large_sample_value")
cat("Effect of one treated day on death by day 28, simulated data\n")
print(round(t_eff, 4))
# Cumulative weights at day 4: stabilised versus unstabilised.
t_wt <- rbind(stabilised = c(mean(last$sw), max(last$sw), ess(last$sw)),
              unstabilised = c(mean(last$w), max(last$w), ess(last$w)))
colnames(t_wt) <- c("mean", "max", "ess")
cat("Cumulative weights at day 4 (3,000 patients), simulated data\n")
print(round(t_wt, 4))
# Risk of death by day 28 under each regime, from the marginal MSM, beside the true risk.
t_risk <- rbind(never_treated = c(coef(msmm)[[1]], r_never),
                always_treated = c(coef(msmm)[[1]] + 3 * coef(msmm)[[2]], r_never + 3 * g1))
colnames(t_risk) <- c("msm_estimate", "true_value")
cat("Risk of death by day 28 by treatment regime, simulated data\n")
print(round(t_risk, 4))
cat("These true values describe the simulation design (simulated data), not any real drug.\n")

# ---- Follow-up length, the marginal MSM intercept, and how many days each patient was treated ----
# Death is counted up to day 28.
canon("truth.design.death_day", 28, 0)
# The intercept of the marginal MSM is the risk of death by day 28 if no patient were ever treated.
canon("truth.msm.g0_marginal", r_never)
# Share who died, and patients treated on one or two of the three days.
canon("death_risk", mean(wide$Y))
canon("treated_one_day", sum(wide$cumA == 1), 0)
canon("treated_two_days", sum(wide$cumA == 2), 0)
print(table(treated_days = wide$cumA))
cat("Simulated data (ข้อมูลจำลอง): not evidence about any real drug or patient.\n")

> cat("Effect of one treated day on death by day 28, simulated data\n")
Effect of one treated day on death by day 28, simulated data

> print(round(t_eff, 4))
                       per_day      lo      hi large_sample_value
naive_pooled            0.0859  0.0631  0.1087             0.0860
naive_pooled_adj_L_V   -0.0128 -0.0363  0.0107            -0.0119
naive_cumulative        0.0676  0.0498  0.0854             0.0641
naive_cumulative_adj    0.0108 -0.0057  0.0273             0.0005
msm_stabilised_weights -0.0573 -0.0806 -0.0340            -0.0500
msm_marginal           -0.0566 -0.0803 -0.0329            -0.0500

> t_wt <- rbind(stabilised = c(mean(last$sw), max(last$sw),
+     ess(last$sw)), unstabilised = c(mean(last$w), max(last$w),
+     ess(last$w)))

> colnames(t_wt) <- c("mean", "max", "ess")

> cat("Cumulative weights at day 4 (3,000 patients), simulated data\n")
Cumulative weights at day 4 (3,000 patients), simulated data

> print(round(t_wt, 4))
               mean     max      ess
stabilised   0.9994 10.3018 1694.320
unstabilised 7.9648 94.9382 1456.026

> t_risk <- rbind(never_treated = c(coef(msmm)[[1]],
+     r_never), always_treated = c(coef(msmm)[[1]] + 3 * coef(msmm)[[2]],
+     r_never + 3 * g1))

> colnames(t_risk) <- c("msm_estimate", "true_value")

> cat("Risk of death by day 28 by treatment regime, simulated data\n")
Risk of death by day 28 by treatment regime, simulated data

> print(round(t_risk, 4))
               msm_estimate true_value
never_treated        0.4888     0.4725
always_treated       0.3189     0.3225