2.1 medfit::MediationData (the fitting contract)

fit_mediation() / extract_mediation() return a MediationData whose @estimates is a named numeric vector and @vcov a named matrix:

estimates: m_(Intercept), m_X, m_C, y_(Intercept), y_X, y_M, y_C, a, b, c_prime

The convenience aliases a (= m_X), b (= y_M), c_prime (= y_X) and the matching @vcov dimnames are the interface every downstream consumer resolves by name.

2.2 RMediation (the inference contract)

Namespace note: the package is RMediation (capital R-M), not rmediation — the source repo is data-wise/rmediation but DESCRIPTION: Package: RMediation. Code uses RMediation::.

• ci_mediation_data(mu, level, type, n.mc) — Monte-Carlo / distribution-of-the- product CI from a single named MediationData. Resolves a/b strictly by name (errors if unnamed). This is the type = "mc" path.
• mbco(h0, h1, ...) — likelihood-ratio MBCO, currently OpenMx-only; it has no per-imputation / MI entry point (see §4).

2.3 Dependency direction (no cycle)

missingmed  ->  medfit       (fit_mediation, MediationData)
missingmed  ->  RMediation   (ci_mediation_data, medci; ci for X:M models)
RMediation  ->  medfit       (consumes MediationData)
medsim          (Suggests, later phases)

missingmed imports both medfit and RMediation; neither depends back on missingmed, so the graph is acyclic.

Worked example: pooling three imputations

Time: about 3 minutes. Arithmetic only – no calculus, no software.

After this you can compute \(T\) by hand and say what fraction of your information the missing data cost you.

The punchline first. Three imputations whose estimates differ by only \(\pm 0.10\) turn a standard error of 0.100 into 0.153. That 53% inflation is the price of the missing data, and computing it is all pool() does.

Symbols in this example:

The three fits:

Step 1 – average the estimates. \(\bar Q = \tfrac{0.50+0.60+0.70}{3} = \mathbf{0.60}\)

Step 2 – average the within-imputation variances. \(\bar U = \tfrac{0.01+0.01+0.01}{3} = \mathbf{0.01}\)

Step 3 – measure the disagreement between imputations. \(B = \tfrac{1}{3-1}\big[(-0.10)^2 + 0^2 + (0.10)^2\big] = \tfrac{0.02}{2} = \mathbf{0.01}\)

Step 4 – combine the two sources. \(T = \bar U + \left(1+\tfrac13\right)B = 0.01 + \tfrac43(0.01) = \mathbf{0.0233\overline{3}}\), so \(SE = \sqrt{T} = \mathbf{0.1528}\).

Two more numbers fall straight out:

• Fraction of missing information \(\gamma = (1+\tfrac1m)B/T = 0.01\overline{3}/0.023\overline{3} = \mathbf{0.571}\). Over half the information about \(a\) is missing.
• Degrees of freedom \(\nu = (m-1)\big(1 + \bar U/[(1+\tfrac1m)B]\big)^2 = 2(1+0.75)^2 = \mathbf{6.125}\). Three imputations buy about six df, not three.

Explain this one to yourself before moving on. Step 4 multiplies \(B\) by \((1+\tfrac1m)\) instead of using \(B\) alone. Why should a larger \(m\) make that correction smaller?

Your turn – same four steps, one line blanked. Now \(Q = (0.20,\ 0.30,\ 0.40)\), still with \(U_i = 0.01\) and \(m = 3\). \(\bar Q = 0.30\); \(B = \underline{\hphantom{0.01}}\); \(T = \underline{\hphantom{0.0233}}\).

Why this cannot be done on \(\hat a\hat b\) directly. Every step above is linear in \(Q\). The indirect effect \(ab\) is not, so pool() deliberately stops here: it pools \((a, b)\) and their covariance and hands the nonlinear step to RMediation. Pooling \(\hat a_i \hat b_i\) across imputations and attaching a normal interval commits two errors at once.

You can now read res@tidy_table’s var_w, var_b and var_tot columns as steps 2, 3 and 4 of the arithmetic above.

The table also carries a Rubin-pooled Wald test for each coefficient: statistic, df, riv, fmi and p_value. df is the Barnard–Rubin (1999) small-sample df, using each model’s complete-data df. The mediator model’s df applies to the m_* rows and the outcome model’s to the y_* rows, and it is infinite for a binomial or poisson model. With a single fit (m = 1, as under IPW) the test is the ordinary one, with riv = fmi = 0. These p-values test one path at a time. A small p_value on a and on b is not a test of mediation: \(ab\) is nonlinear and its null is non-regular, which is the reason infer() exists (Section 4).

The pooled MediationData is built by copy-modifying a per-imputation one (S7 is copy-on-write), replacing only @estimates =\(\bar Q\) and @vcov =\(T\) — so all medfit metadata (roles, predictors, families) is inherited and the aliases stay name-addressable for ci_mediation_data().

fit <- run(md_mi)            # m named MediationData
res <- pool(fit)             # pooled named MediationData + within/between/total vcov
res@tidy_table               # + statistic, df, riv, fmi, p_value (per path)

3A.1 Estimand

With mediator model \(M = \beta_0 + aX + \sum_k \beta_k C_k + e_M\) and outcome model \(Y = \gamma_0 + c'X + bM + \theta_3 XM + \dots + e_Y\), the effect of a one-unit increase in \(X\) through \(M\), with \(X\) held at \(x\) in the outcome model, is

\[ a\,(b + \theta_3 x). \]

It depends on \(x\), so infer(type = "mc") requires treatment_level on such a model and refuses it on a model without X:M. For a 0/1 treatment, treatment_level = 1 gives the total natural indirect effect and 0 the pure natural indirect effect. The draws come from the pooled estimates of \((a, b, \theta_3)\) and the matching block of the pooled total covariance \(T\). That block is not block-diagonal: the between-imputation part couples the mediator and outcome models, because both are fitted to the same completed dataset.

3A.2 The four-way decomposition, pooled

medfit returns an InteractionMediationData per imputation, carrying the four-way decomposition (VanderWeele 2014):

with \(\text{nde} = \text{cde} + \text{int\_ref}\), \(\text{nie} = \text{int\_med} + \text{pie}\), and the total effect their sum. \(m^*\) is the mediator level at which the CDE is evaluated (medfit’s default is 0). \(m_{\text{ref}} = \beta_0 + \sum_k \beta_k \bar C_k\) is the mediator’s expected value at the reference treatment level with covariates at their means.

pool() does not average these components one by one. It recomputes all of them from one pooled reference profile:

1. \(a, b, c', \theta_3\) are the Rubin estimates \(\bar Q\).
2. \(m^*\) must be the same in every imputation; otherwise pool() stops, because the decompositions refer to different mediator levels.
3. \(m_{\text{ref}}\) changes across imputations whenever a covariate is imputed, since \(\bar C_k\) is a mean over the completed data. It is pooled as the mean of the per-imputation values. Each value is read back from medfit’s own int_ref as \(m^* + \text{int\_ref}/\theta_3\). When an imputation’s \(\theta_3\) is exactly 0, int_ref carries no information about \(m_{\text{ref}}\), and it is rebuilt from \(\beta_0\) and the covariate means medfit stores (medfit >= 0.4.0).
4. \(\text{int\_ref} = \bar\theta_3(\bar m_{\text{ref}} - m^*)\), and the other components follow from the formulas above.

The mean of the per-imputation int_ref values would be the mean of the products \(\theta_{3,i}\, m_{\text{ref},i}\). That is not \(\bar\theta_3 \bar m_{\text{ref}}\), and the pooled object would then fail medfit’s own identities (for example \(\text{nde} = \text{cde} + \text{int\_ref}\) built from pooled parts). Recomputing from one profile keeps every identity exact.

In the tidy table, the theta3 row takes the df of its source row (y_X:M), and b0 that of m_(Intercept), like the a, b and c_prime aliases in section 3.

The check below fits a model whose covariate is imputed, so \(m_{\text{ref}}\) differs across imputations, and fails the build if the pooled decomposition does not come from one profile:

library(missingmed)
set.seed(8)
n <- 300
d <- data.frame(X = rbinom(n, 1, 0.5), C = rnorm(n))
d$M <- 0.5 * d$X + 0.4 * d$C + rnorm(n)
d$Y <- 0.2 * d$X + 0.3 * d$M + 0.25 * d$X * d$M + 0.3 * d$C + rnorm(n)
d$M[sample(n, 40)] <- NA
d$C[sample(n, 50)] <- NA
imp_xm <- mice::mice(d, m = 5, method = "norm", printFlag = FALSE, seed = 8)
fit_xm <- run(set_md_mediation(imp_xm, Y ~ X * M + C, M ~ X + C,
  treatment = "X", mediator = "M"
))
p <- missingmed::pool(fit_xm)@pooled
fits <- fit_xm@per_imputation
m_ref <- vapply(fits, function(f) f@m_star + f@int_ref / f@interaction, numeric(1))
round(m_ref, 4) # differs across imputations: C is imputed
#>       1       2       3       4       5 
#> -0.0918 -0.1236  0.0431 -0.0605 -0.0239
stopifnot(
  all.equal(p@int_ref, p@interaction * (mean(m_ref) - p@m_star)),
  all.equal(p@nde, p@cde + p@int_ref),
  all.equal(p@nie, p@int_med + p@pie)
)
c(pooled_int_ref = p@int_ref,
  mean_of_int_ref = mean(vapply(fits, function(f) f@int_ref, numeric(1))))
#>  pooled_int_ref mean_of_int_ref 
#>     -0.01289910     -0.01314089

The last line shows the two candidates side by side. They differ, and only the first is consistent with the pooled \(\theta_3\).

3A.3 What stays refused

Only one product term is supported, and only in the outcome model. X:C, M:W, X:M:W, transforms of X or M and an X:M term in a non-Gaussian model are refused before fitting (the Supported models article lists every rule). For moderated models, mbco_d4() tests the indirect effect on the completed datasets directly.

4.1 Why MBCO does not commute with Rubin’s rules

MBCO tests \(H_0: a b = 0\). Since \(a b = 0 \iff a = 0 \lor b = 0\), the constrained log-likelihood is the branch union:

\[ T = 2\big[\ell_{\text{full}} - \max(\ell_{a=0},\ \ell_{b=0})\big]. \]

That max() is non-linear, so the MBCO statistic of the pooled estimate is not the pool of per-imputation MBCO statistics. You cannot pool first and test second — you need the per-imputation fits.

4.2 D4 combination

Hence MDMediationFit retains the per-imputation list (exposed by per_imputation_list()), and infer(type = "mbco") combines the per-imputation LRT statistics with the D4 rule (Chan & Meng 2022; Grund et al. 2023):

\[ d_S = \tfrac{\text{LRT(stacked data)}}{K},\quad r_4 = \max\!\Big(0, \tfrac{K+1}{k(K-1)}(\bar d - d_S)\Big),\quad D_4 = \tfrac{d_S}{k(1+r_4)} \sim F_{k,\nu},\quad \nu = k(K-1)\big(1 + \tfrac1{r_4}\big)^2 \]

The ariv argument sets which per-imputation statistics \(d_i\) enter \(\bar d\). With ariv = "fixed" (the default), every \(d_i\) is computed on the branch the stacked constrained fit chose, so all imputations share the stacked fit’s \(k\). With ariv = "own", each \(d_i\) uses that imputation’s own winning branch (the D4 and r4 of missingmed 0.4.0 on full-rank designs); it errors when the branches remove different numbers of parameters, and it can be liberal at high missingness. See vignette("mbco-mi") for both formulas.

4.3 Hosting decision

D4-MBCO under multiple imputation lives in missingmed (R/mbco_mi.R), and RMediation keeps complete-data MBCO (author decision, 2026-10-07; issue #19). The code was ported from a research prototype and reproduced its statistic, \(r_4\) and branch choice exactly (max abs diff \(\approx 5\times10^{-11}\) across design cells). The denominator degrees of freedom \(\nu\) were later changed to Chan and Meng’s eq. 2.15, as in mitml::testModels(method = "D4"), so \(\nu\) and \(p\) now differ slightly from the prototype’s (which used the Li et al. 1991 df). The exported mbco_d4() runs the same test on a plain list of completed data frames, so other packages can call it directly.

5.1 Pipeline

md <- set_md_mediation(df, Y ~ X + M + C, M ~ X + C, treatment = "X",
  mediator = "M", method = "ipw",
  weight_formula = NULL, weight_stabilize = TRUE, weight_trim = 0.99,
  se_type = "sandwich")
res <- pool(run(md))         # MDMediationFit(m = 1) -> MDMediationResult (B = 0)
infer(res, type = "mc")

Because there is a single weighted fit, MDMediationFit@m = 1 and pool() reduces to \(\bar Q = Q_1,\ B = 0,\ T = \bar U\) — the result is structurally identical to an MI result, so nothing downstream changes.

5.2 Weight estimation

Let \(R = 1\) mark a complete case over the model variables.

• Sequential (chain-rule) factorization (default, and a named list). With incomplete variables \(V_1, V_2, \dots\) in a given order, \(P(\text{complete}\mid Z) = P(V_1\text{ obs}\mid Z)\, P(V_2\text{ obs}\mid V_1\text{ obs}, Z)\cdots\), each factor a logistic model fitted on the rows where the earlier variables are observed. It is an identity, so it is correct whenever each factor is logistic, whether the variables go missing separately, together or one after another. A factor observed on every fitted row (the second variable under simultaneous missingness) is 1 and fits no model. weight_formula = NULL orders the variables by ascending share of missing values (ties in model-variable order) and predicts with the fully observed model variables (the MAR drivers); a list gives the order and the predictors. It is the standard factorization for monotone missingness (Robins, Rotnitzky and Zhao 1994, section 7.1; Seaman and White 2013, section 2.3, equations 4 and 5). For one fixed order it is an “observed past” missingness process, which Robins and Gill (1997, section 9) show is a randomized monotone missingness process. The order is an assumption: if the true order differs from subject to subject, the probability of being complete is a sum over orders, not one product (Robins and Gill 1997, section 4). The order matters when the factors are not all logistic. Under a monotone dropout the factor of the variable that drops out second is logistic given the first; in the reversed order the first factor is the probability that both are observed, a product of two logistic probabilities, which is not.
• Joint complete-case model (a single formula). One model, \(\text{logit }P(R=1\mid Z)\), with the formula’s predictors. Correct when the incomplete variables go missing together or only one is incomplete; biased at any \(n\) when they go missing separately, because \(P(R=1\mid Z)\) is then a product of logistics.
• The earlier marginal form (each variable on all rows, weights from the product) is not offered: it is biased at any \(n\) when the variables go missing together or one after another.

What is assumed. With the default predictors, and with any formula built from fully observed variables (the advice throughout), the probability of being complete is assumed to depend on the fully observed model variables only, not on the observed values of the other incomplete variables. This is the usual restriction (Seaman and White 2013, appendix 2: “commonly restricted to fully observed predictors”). For non-monotone data it is strictly stronger than missing at random, where the probability of a pattern may depend on any variable observed in that pattern (Sun and Tchetgen Tchetgen 2018, section 3.1 of the arXiv version). Modeling that general case is hard (Robins and Gill 1997; Sun and Tchetgen Tchetgen 2018) and is not attempted here.

Stabilized weights (weight_stabilize = TRUE, default) use a treatment-only numerator to reduce variance. The numerator needs no correct model of its own; the stabilized estimator is consistent when the outcome (mean) regression is correct (Seaman and White 2013, section 6): \[ w_i = \frac{\hat P(R=1\mid X_i)}{\hat P(R=1\mid Z_i)} . \]

5.3 Variance: sandwich vs model

IPW’s weighted-GLM model-based vcov is optimistic (it ignores that weights were estimated). se_type = "sandwich" (the IPW default) uses heteroskedasticity-consistent sandwich::vcovHC instead. This is implemented in medfit via an injectable vcov_fun threaded through extract_mediation() (so the named-vcov assembly stays in one place); missingmed just passes se_type through.

MBCO for IPW is not implemented (a weighted LRT is methodologically distinct); infer(type = "mbco") on an IPW fit errors with guidance.

5A.1 The scale caveat

For a non-identity link, a and b are link-scale coefficients, so their product – and the Monte-Carlo interval around it – is on the link scale too:

Two consequences worth stating plainly:

1. a*b under a non-identity link is not a risk difference, a risk ratio or an odds ratio, and exp(a*b) does not produce one. Reporting it as such is a category error, not a rounding issue.
2. Rubin’s rules are applied to the coefficients, which is standard. Pooling a transformed quantity (an odds ratio, say) is a different and much less well-behaved operation; pool() deliberately does not do it.

Where a response-scale contrast is wanted – a natural direct/indirect effect on the risk difference or risk-ratio scale – the mediation formula must be applied to the fitted models, which is outside missingmed’s orchestration remit. That is a fitting-layer concern (medfit), not a pooling one.

5A.2 IPW plus a binomial family: the non-integer-successes warning

stats::glm() emits non-integer #successes in a binomial glm! whenever a binomial fit receives non-integer prior weights. IPW weights are inverse sampling probabilities, not trial counts, so the warning fires on every method = "ipw" fit with a binomial family. It is a false alarm: the weighted score equations being solved are exactly the ones IPW specifies, and the point estimates and sandwich SEs are unaffected. It is left unsuppressed rather than silenced, so that a genuine non-integer-response mistake still surfaces.

5B.1 Method

Delta-adjusted imputation in the pattern-mixture sense (van Buuren, FIMD §9.2; Leacy et al. 2017). Per rung: re-impute with the shift applied, then run the existing pipeline unchanged.

For a continuous target the canonical procedure imputes under MAR and then adds the constant to the imputed values (Hayati Rezvan et al. 2018). For a binary target an additive shift on a drawn 0 or 1 is meaningless; the shift belongs on the imputation model’s linear predictor instead. Which construction runs depends on the target’s imputation method:

NARFCS is Tompsett et al. (2018); the mice methods are Moreno-Betancur, van Buuren & White (2020). For a constant delta on a norm target the first and third rows would give identical draws — the offset is a scalar added to a normal draw’s linear predictor, which is arithmetically the same as adding it to the draw — so a numeric delta stays on post, and mnar.norm is used only where it adds something: a covariate-varying ums. Only an exact method match is delegated: swapping norm.boot for mnar.norm, say, would change the imputation method itself, and delta = 0 would stop reproducing the MAR analysis. print() and tidy() state the mechanism and scale that ran.

A delegated target also accepts a covariate-varying delta through ums, passed verbatim to NARFCS: ums = c("0", "1 + 0.5*C") offsets each row by \(1 + 0.5\,C\). A string grid has no numeric ordering, so summary() does not compute a tipping point for it.

Re-imputation rebuilds from the mids object’s stored settings, never from its recorded call — a mids built inside a function stores a call referencing that function’s local variables, which does not resolve anywhere else.

sens <- sensitivity_mnar(md, delta = c(0, -0.5, -1, -1.5, -2))
tidy(sens)
summary(sens)   # adds the tipping point, if the grid contains one

5B.2 What delta actually is — read before choosing a value

delta is a conditional sensitivity parameter (CSP): a difference conditional on all remaining variables and their missingness indicators. The quantity an analyst can actually reason about — “non-respondents average half a point lower” — is a marginal one (MSP). They are different numbers.

Supplying an elicited MSP where a CSP is expected is the standard failure mode of this method. Tompsett et al. (2018) measure it: coverage fell from 95% to 49.3% when an elicited MSP was inserted directly as the CSP.

missingmed cannot make the specified parameter marginal, but it can show you the realized one. Every rung reports msp. The run below uses a data set where M and Y are both incomplete and M’s missingness depends on X:

library(missingmed)
set.seed(9)
n <- 400
C <- rnorm(n)
X <- rbinom(n, 1, 0.5)
M <- 0.5 * X + 0.4 * C + rnorm(n)
Y <- 0.45 * M + 0.2 * X + 0.3 * C + rnorm(n)
d <- data.frame(X, M, Y, C)
d$M[runif(n) < plogis(-1.6 + 1.6 * X + 0.4 * C)] <- NA
d$Y[runif(n) < plogis(-1.8 + 0.5 * X + 0.5 * C)] <- NA
meth <- mice::make.method(d)
meth["M"] <- "norm"
imp <- mice::mice(d, m = 5, maxit = 5, method = meth, seed = 9, printFlag = FALSE)
md <- set_md_mediation(imp, Y ~ X + M + C, M ~ X + C,
  treatment = "X", mediator = "M"
)
grid <- c(0, -0.5, -1, -1.5, -2)
sens_tab <- tidy(sensitivity_mnar(md, delta = grid, n.mc = 1e4, seed = 9))
knitr::kable(
  data.frame(`delta (CSP)` = grid, `msp (realized)` = round(sens_tab$msp, 3),
    estimate = round(sens_tab$estimate, 3), check.names = FALSE)
)

Two things to read off that table. First, msp is non-zero at delta = 0 (here 0.237): under MAR, imputed values legitimately differ from observed ones when missingness depends on covariates. Second, msp is not simply delta shifted: its increments differ from the grid’s by up to 0.047, a chained-equations effect. With a single incomplete variable the increments coincide exactly; with more, part of the shift circulates back through the other imputation models.

5B.3 Reading a tipping point

summary() reports the delta at which the interval first includes zero. Judge plausibility on the realized msp, not on the CSP: Tompsett et al. note a CSP-scale tipping point has no direct clinical interpretation. Resseguier et al. add the constraint that a tipping value must “correspond to reasonable hypotheses supported by epidemiologic evidences” — a tipping point that requires an implausible departure from MAR is reassurance, not a warning.

The figure below plots the run from 5B.2. It shows why the resolution warning matters: the interval’s lower bound crosses zero near -1.12, but the grid c(0, -0.5, -1, -1.5, -2) can only report -1.5.

5B.4 Binary targets: when is a log-odds delta of 1 plausible?

On the mnar.logreg route delta is a conditional log-odds difference: at delta = 1 the odds that a non-respondent’s value is 1 are \(e \approx 2.72\) times what MAR imputation implies, given every other variable. The CSP/MSP caution of 5B.2 applies with an extra twist — the two sit on different scales. delta is log-odds; the reported msp is a prevalence difference on the probability scale.

The same log-odds step moves prevalence by different amounts depending on where it starts. The chunk below imputes a binary mediator with mnar.logreg at three log-odds deltas and reports the mean imputed prevalence (it runs when the vignette is built, so the numbers always come from the code shown):

set.seed(7)
n <- 600
C <- rnorm(n)
X <- rbinom(n, 1, 0.5)
M <- rbinom(n, 1, plogis(-0.3 + 0.8 * X + 0.3 * C))
Y <- 0.5 * M + 0.2 * X + 0.3 * C + rnorm(n)
d <- data.frame(X, M, Y, C)
d$M[runif(n) < 0.35] <- NA # missing completely at random
meth <- mice::make.method(d)
meth["M"] <- "mnar.logreg"
prev <- vapply(c(0, 1, 2), function(delta) {
  imp <- suppressWarnings(mice::mice(d,
    m = 5, maxit = 5, method = meth, seed = 7, printFlag = FALSE,
    blots = list(M = list(ums = as.character(delta)))
  ))
  mean(unlist(imp$imp$M))
}, numeric(1))
knitr::kable(
  data.frame(`log-odds delta` = c("0 (MAR)", "1", "2"),
    `imputed prevalence` = round(prev, 3), check.names = FALSE)
)

The first step added 24.5 points, the second 16.2 — the logistic curve flattens toward 1. So:

• A delta of 1 is plausible when external evidence — a follow-up of non-respondents, a validation subsample, a substantive account of why people with the outcome stop responding — supports a conditional odds ratio of roughly 2.7 among the missing. Judge that against the realized msp, not against the delta.
• It is not plausible by default. An odds ratio of 2.7 is a large conditional effect. When the MAR prevalence is already near 0 or 1, even a modest delta pushes the imputed values toward all-0 or all-1, and the curve mostly reflects that the grid ran out of room.
• delta = 0 still reproduces the MAR analysis exactly, so the curve is anchored at the analysis you would otherwise report.

5B.5 Limitations

7A.1 Three common confusions

“delta = 0 should give msp = 0.” No. Under MAR, imputed values legitimately differ from observed ones whenever missingness depends on covariates, because the imputed subgroup has a different covariate mix. What \(\delta = 0\) guarantees is that the analysis reproduces the MAR result exactly – not that the marginal shift is zero. The table in section 5B.2 shows msp = 0.244 at delta = 0.

“se_type = "sandwich" corrects for having estimated the weights.” No. The sandwich is heteroskedasticity-consistent and treats \(\hat\pi\) as fixed. When the missingness model is correct, ignoring the estimation of the weights makes the variance estimate conservative for each coefficient: the asymptotic variance with estimated selection probabilities is no larger than with the true ones (Robins, Rotnitzky and Zhao 1994, section 6.1 and proposition 6.1; the same for propensity-score weights in Lunceford and Davidian 2004, section 3.1; stated for IPW generally in Seaman and White 2013, section 8). The intervals are then wider than necessary, not too narrow.

Three limits. The result is asymptotic and assumes a correct missingness model; with a misspecified or heavily weighted model the interval is centered in the wrong place and can under-cover. It is proved for a regression coefficient, not for a product of two, so for \(a b\) it is an extrapolation: in this package’s simulations the variance that accounts for estimated weights was 1% to 3% narrower, and in no setting did it repair an interval that under-covered. And Robins, Rotnitzky and Zhao (1994, after proposition 6.1) warn that a richly parameterized missingness model biases the estimated variance in finite samples and suggest the bootstrap.

“IPW and MI should agree; if they disagree, one is broken.” No. They target the same estimand under MAR but use different information: MI uses partially observed rows, IPW discards them and reweights what remains. Disagreement beyond Monte-Carlo noise is diagnostic – usually of a weight model missing a MAR driver, or an imputation model incompatible with the substantive model – not of a bug.