Canonical formal derivation for the JAIGP theory paper (setup, assumptions, proof)
Multiple imputation (MI) replaces each missing entry with several plausible draws, lets the analyst apply complete-data tools to each completed data set, and combines the results (Rubin 1987). Many of those complete-data tools — likelihood-ratio statistics, information criteria, profile likelihoods, deviance-based fit indices — are built, directly or indirectly, on the imputation-averaged complete-data log-likelihood: the Q-function of the EM and congeniality literatures (Schafer 1997; Meng 1994). Write \(\bar Q_\infty(\theta)\) for its idealized infinite-imputation limit. One would like \(\bar Q_\infty\) to behave like the complete-data log-likelihood \(\ell_{\text{com}}\) that would have been computed had nothing been missing. It does not. Evaluated at the observed-data estimator and measured against the complete-data gold standard \(\ell_{\text{com}}(\hat\theta_{\text{com}})\), the Q-function is systematically biased. This document derives that bias exactly, to leading order, and names the two mechanisms behind it.
Let \(\mathrm{RIV} = I_{\text{obs}}^{-1} I_{\text{mis}\mid\text{obs}}\) denote the matrix relative increase in variance due to missingness — Rubin’s scalar \(r\) in matrix form. Under congenial, proper MI, taking the complete-data MLE \(\hat\theta_{\text{com}}\) as the target, the Q-function deviance bias is \[ \mathbb{E}\big[\bar Q_\infty(\hat\theta_{\text{obs}}) - \ell_{\text{com}}(\hat\theta_{\text{com}})\big] = \underbrace{+\operatorname{tr}(\mathrm{RIV})}_{\text{imputation bias}} \;\;\underbrace{-\tfrac12\operatorname{tr}(\mathrm{RIV})}_{\text{estimation mismatch}} \;=\; +\tfrac12\operatorname{tr}(\mathrm{RIV}), \] to \(O(n^{-1})\) (on the deviance scale, double this). The first term is imputation bias: averaging \(\ell_{\text{com}}\) over completions drawn from the imputation model overstates the fit by \(\operatorname{tr}(\mathrm{RIV})\). The second is estimation mismatch: MI’s estimand is \(\hat\theta_{\text{com}}\), whereas \(\bar Q_\infty\) is centered at the observed-data maximizer \(\hat\theta_{\text{obs}}\), and the gap between them costs \(\tfrac12\operatorname{tr}(\mathrm{RIV})\). The two do not cancel — they leave a net upward bias of half the trace. This \(\tfrac12\operatorname{tr}(\mathrm{RIV})\) holds to leading order under MCAR (the design-imbalance term is \(O(n^{-1})\)) for any regular likelihood whose conditional missing-data entropy depends on the estimated parameters — concretely, an estimated scale or covariance, as in the multivariate normal we instantiate throughout. For canonical generalized linear models the same statement holds up to a computable base-measure offset \(\Delta_h\), zero for the Bernoulli and nonzero for the Poisson (Section 4.6). That dependence is essential: with a known scale (a location-only fit) the conditional entropy is parameter-free and the net bias collapses to \(\mathbb{E}[T]=0\) under FIML and \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\) under proper MI (scope remark after the theorem). Under MAR Term A carries one further \(O(1)\) piece — a design-imbalance term \((A)+(C)\), nonzero exactly when the missing and observed units have different conditioning-variable distributions and zero under MCAR — developed in Section 4.4.
Whatever is built on \(\bar Q_\infty\) inherits this bias. Two applications motivate the paper; both are developed below (Section 5, Section 6).
The result sits upstream of the MI test-calibration literature. Meng and Rubin’s combining rule for likelihood-ratio tests (Meng and Rubin 1992), its refinement by Chan and Meng (Chan and Meng 2022), and Chan’s general stacked-MI tests (Chan 2022) all concern how to calibrate the reference distribution of an MI test statistic, taking the statistic’s numerator as given. We address the bias of that numerator. The two are complementary: calibration makes the reference distribution right; the present correction makes the statistic referred to it centered. On the model-selection side, the complete-data-discrepancy criteria PDIO and AICcd (Cavanaugh and Shumway 1998) attach a \(2\operatorname{tr}(\mathrm{RIV})\) missing-data penalty to the complete-data deviance under EM-based estimation; (Shimodaira and Maeda 2018) showed this doubles the correct value and derived \(\mathrm{AIC}_{x;y}\) with the halved penalty \(\operatorname{tr}(\mathrm{RIV})\). Our net correction reproduces exactly that \(\operatorname{tr}(\mathrm{RIV})\) — so as a penalty the headline IC result is not new: it is \(\mathrm{AIC}_{x;y}\). What is new is threefold: (i) the decomposition of that penalty into an imputation-bias term \(+\operatorname{tr}(\mathrm{RIV})\) and an estimation-mismatch term \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\) (MI’s target is \(\hat\theta_{\text{com}}\), not the EM pivot \(\hat\theta_{\text{obs}}\)); (ii) the MAR design-imbalance correction \((A)+(C)\), which lies beyond the i.i.d. fixed-pattern setting of (Shimodaira and Maeda 2018) — named there (§8) as future work in its concluding remarks (“combine a missing mechanism with other sampling mechanisms, such as the covariate-shift … problem”), the covariate-shift machinery being (Shimodaira 2000); and (iii) the extension from the deterministic EM Q-function to proper multiple imputation. Claeskens and Consentino’s missing-covariate AIC (Claeskens and Consentino 2008) targets a different (EM-adjusted-likelihood) discrepancy with a Takeuchi-type penalty, and so sits outside this line. A third route — Hens, Aerts and Molenberghs’s weighted AIC (Hens, Aerts, and Molenberghs 2006) — reaches the same complete-data-equivalent target by inverse-probability (Horvitz–Thompson) reweighting of the complete cases rather than a deviance correction: it needs estimated selection probabilities, keeps the penalty at \(2q\), and explicitly leaves an imputation-based criterion — which the present \((A)+(C)\) term bias-corrects — as open. The design-imbalance bias it diagnoses (its complete-case discrepancy shifts by \(\operatorname{Cov}(y,\pi)/\mathbb{E}[\pi]\) under MAR) is the same phenomenon \((A)+(C)\) captures intrinsically within MI.
Section 2 fixes notation: the model and missingness, the two likelihoods and their maximizers, the three Fisher informations and the identity binding them, RIV, the Q-function, and the bias decomposition \(T = A + B\). Section 3 collects the regularity, congeniality, and properness conditions together with the asymptotic lemmas the proof rests on. Section 4 states the theorem, sketches the proof strategy, and gives the proof: Term B (estimation mismatch) is an elementary Taylor argument, while Term A (imputation bias) is a second-order expansion turning on a cancellation we isolate as a lemma; a closing subsection carries that cancellation beyond the multivariate normal, to canonical generalized linear models, where a base-measure offset survives it (Section 4.6). Section 5 develops the likelihood-ratio application: the closed-form leading-order differential bias of the MI deviance numerator at the null, why the naive difference of the two models’ \(\operatorname{tr}(\mathrm{RIV})\)s overstates it, and the noise collapse that makes the comparison-relevant correction precisely estimable. Section 6 develops the model-selection application: the corrected criterion and its place in the missing-data information-criterion literature, the pre-registered selection evidence, and a graded contrast — mean, two moments, all moments — of how much of the complete-data criterion distribution a correction can match, with the null/noncentral split derived in advance and then measured. Section 7 reproduces the cited source passages verbatim, and Section 8 verifies the trace identities symbolically and the normal instantiation numerically, executed inline at render time.
A note on provenance. Because this is a human-prompted AI derivation, every external claim is sourced. The main column carries the mathematics; the right margin carries the provenance — the source’s own equation or theorem number, the reproduced statement, and a line reference into the local literature/ copy (with a short verbatim snippet, so the pointer survives re-pagination). Steps proved here rather than cited are flagged [derived]; those additionally checked symbolically in Section 8 are flagged [mc]. Independently, the pre-registered Monte Carlo study in verification/ estimates \(\mathbb{E}[A],\mathbb{E}[B],\mathbb{E}[T]\) across configurations and compares them with the analytic \(\tfrac12\operatorname{tr}(\mathrm{RIV})\). One honest caveat on that comparison: the study was pre-registered against \(\tfrac12\operatorname{tr}(\mathrm{RIV})\) alone, before the \(O(1)\) MAR term \((A)+(C)\) had been derived, and its per-replicate statistic is heavy-tailed with standard deviation \(\propto\sqrt n\) — so at fixed replication the per-cell estimates at larger \(n\) cannot resolve an \(O(1)\) offset and one cell (\(N{=}1000\)) misses the uncorrected target at \(2.6\) standard errors. Pooling the four sample sizes by inverse variance gives \(\hat T = 2.43 \pm 0.26\) against the pre-registered target \(2.77\) and the \((A)+(C)\)-augmented mean \(2.42\)–\(2.55\): the pooled pre-registered data corroborate the completed theorem (verification/w1-h1-pooled-reconciliation.R; the criterion amendment is documented in the pre-registration’s version history). A second candor note at the component level: the individual terms \(\mathbb{E}[A]\) and \(\mathbb{E}[B]\) miss their Monte-Carlo targets by roughly \(9\)–\(20\%\) at the sample sizes studied (\(B\) overshoots its \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\) target by \(-0.3\) to \(-0.55\) with tight Monte-Carlo error; verification/w1-h1-pooled-reconciliation.R, part 4) — only the total \(T\) and the paired contrasts are Monte-Carlo-tight at finite \(n\), which is why every quantitative claim in this document is stated for totals and paired differentials, never for components.
This section fixes the objects the theorem is about: the data and model, the two likelihoods and their maximizers, the three Fisher informations and the identity binding them, the relative increase in variance, the Q-function, and the bias decomposition \(T = A + B\). The conditions and asymptotic results the proof rests on — regularity, congeniality, properness, and the expansions — are collected separately in Section 3. Definitions carry a label (D1, D2, …) so that the margin provenance and the later proof can point back to them.
D1 (data and model). Let \(Y \in \mathbb{R}^{N\times p}\) collect \(N\) complete observations, the rows drawn independently from a \(p\)-variate normal, \(Y_i \sim \mathcal{N}_p(\mu,\Sigma)\) with \(\Sigma \succ 0\). The estimand is the mean–covariance parameter \(\theta = (\mu, \operatorname{vech}\Sigma)\), of dimension \(k = p + \tfrac12 p(p+1)\). Missingness, introduced next, will keep us from observing \(Y\) in full.
Convention. Complete data \(Y=(Y_{\text{obs}},Y_{\text{mis}})\) follows Schafer (1997), :1477, where \(\ell(\theta\mid Y)\), \(\ell(\theta\mid Y_{\text{obs}})\), \(P(Y_{\text{mis}}\mid Y_{\text{obs}},\theta)\) appear and Little and Rubin (2019), :438. Meng’s congeniality framework writes \(Z_c=\{X,Y_{\text{inc}},I\}\), \(Z_o=\{X,Y_{\text{obs}},I,R_{\text{inc}}\}\) with \(X\) = always-observed covariates (Meng 1994), :95 — a different \(X\), which is why we avoid the symbol. The verification code (verification/00-setup.R) labels the realized complete matrix X; that is our \(Y\).
D2 (missingness). A response indicator \(R\) (with \(1\) marking a missing entry) splits each row into its observed and missing parts, \(Y=(Y_{\text{obs}},Y_{\text{mis}})\). Throughout we assume the missingness is ignorable — missing at random, with the missingness parameters distinct from \(\theta\) — so the mechanism itself need not be modeled. The precise conditions are stated as assumptions in Section 3.
Missingness forces a distinction the whole argument turns on: the likelihood we would have maximized with complete data versus the one we can maximize from what is observed, and their two different maximizers.
D3 (complete-data log-likelihood and MLE). Had \(Y\) been observed in full, we would form \(\ell_{\text{com}}(\theta; Y) = \sum_i \log f_{\mathcal{N}}(Y_i; \mu,\Sigma)\) and maximize it at the complete-data MLE \(\hat\theta_{\text{com}} = \arg\max_\theta \ell_{\text{com}}(\theta; Y)\). This is the gold standard MI aims to recover.
D4 (observed-data log-likelihood and MLE). From the observed entries alone, each row contributes only its marginal density, giving \(\ell_{\text{obs}}(\theta; Y_{\text{obs}}) = \sum_i \log f_{\mathcal{N}}(Y_{i,\text{obs}}; \mu_{O_i},\Sigma_{O_iO_i})\), maximized at the observed-data (FIML) MLE \(\hat\theta_{\text{obs}} = \arg\max_\theta \ell_{\text{obs}}\). The gap between \(\hat\theta_{\text{obs}}\) and \(\hat\theta_{\text{com}}\) is the source of Term B.
Each estimator’s precision is governed by a Fisher information, and the two informations differ by exactly the information the missing data would have supplied.
D5 (complete information). The complete-data Fisher information is \[ I_{\text{com}}(\theta) = -\,\mathbb{E}_{Y\sim\theta}\!\big[\partial^2\ell_{\text{com}}(\theta;Y)/\partial\theta\,\partial\theta^\top\big] = \mathbb{E}\big[S_{\text{com}}(\theta)\,S_{\text{com}}(\theta)^\top\big]. \]
D6 (observed information). The observed-data (expected) Fisher information — the asymptotic precision of \(\hat\theta_{\text{obs}}\) (Lemma L1) — is \[ I_{\text{obs}}(\theta) = -\,\mathbb{E}_{Y_{\text{obs}}\sim\theta}\!\big[\partial^2\ell_{\text{obs}}(\theta;Y_{\text{obs}})/\partial\theta\,\partial\theta^\top\big]. \]
Source. Schafer (1997), :1489: observed information \(I_o=-\partial^2\ell(\theta\mid Y_{\text{obs}})/\partial\theta^2\); we use its expectation, the Fisher form (Vaart (1998), :2320), taken unconditionally over \((Y,R)\) — see Section 2.4. For MVN the naive (MCAR) form is the pattern-mixture sum (fisher_info_obs_mvn); under MAR \(I_{\text{obs}}\) carries an extra mean–covariance cross term (Appendix C) and is estimated by the observed information.
D7 (missing information and the information identity). Their difference is the missing information, \(I_{\text{mis}\mid\text{obs}}(\theta) := I_{\text{com}}(\theta) - I_{\text{obs}}(\theta)\), equivalently the additive identity \[ I_{\text{com}}(\theta) = I_{\text{obs}}(\theta) + I_{\text{mis}\mid\text{obs}}(\theta). \tag{1}\] It is positive semidefinite, \(I_{\text{mis}\mid\text{obs}}\succeq 0\): missingness can only cost information, never create it.
Source. Schafer (1997), missing-information principle: Schafer establishes \(I_c=I_o+I_m\) conditional on \(Y_{\text{obs}}\) via the EM decomposition \(Q(\theta\mid\phi)=\ell_{\text{obs}}+H(\theta\mid\phi)\), with \(I_c=-\partial^2 Q(\theta\mid\theta)/\partial\theta^2\) [:1483], \(I_m=-\partial^2 H(\theta\mid\theta)/\partial\theta^2\) [:1495], eq (3.19) [:1501]. Taking \(\mathbb{E}_{Y_{\text{obs}}}\) gives the identity for the Fisher informations D5-D6 (his conditional \(I_c,I_o,I_m\) average to our \(I_{\text{com}},I_{\text{obs}},I_{\text{mis}\mid\text{obs}}\)). \(I_{\text{mis}\mid\text{obs}}\succeq 0\) [mc] (min eigenvalue \(\approx 0\) on the fixture, Listing 6).
The ratio of missing to observed information is the quantity in which the whole result is written.
D8 (relative increase in variance). The matrix relative increase in variance, \[ \mathrm{RIV}(\theta) := I_{\text{obs}}^{-1}\, I_{\text{mis}\mid\text{obs}}, \tag{2}\] is Rubin’s scalar \(r\) in matrix form — a measure of how much missingness inflates the variance of \(\hat\theta_{\text{obs}}\) relative to the complete-data optimum. Its trace is the scalar that drives the theorem.
Source (name, scalar form). Rubin (1987), :2194: the relative increase in variance due to nonresponse \(r_m = (1+m^{-1})B_m/\bar U_m\); with within \(\bar U_m=\sum_l U_{*l}/m\) [:2167], between \(B_m=\sum_l(\hat Q_{*l}-\bar Q_m)^{t}(\hat Q_{*l}-\bar Q_m)/(m-1)\) [:2173], total \(T_m=\bar U_m+(1+m^{-1})B_m\) [:2179].
An equivalent form, used throughout: substituting \(I_{\text{mis}\mid\text{obs}}=I_{\text{com}}-I_{\text{obs}}\) (D7) and using \(I_{\text{obs}}^{-1}I_{\text{obs}}=I_k\), \[ \mathrm{RIV} = I_{\text{obs}}^{-1}I_{\text{com}} - I_k, \qquad \operatorname{tr}(\mathrm{RIV}) = \operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{com}}) - k, \tag{3}\] so the trace appearing in the theorem is \(\operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{com}}) - k\).
Under congenial proper MI as \(M\to\infty\): \(\bar U_\infty = I_{\text{com}}^{-1}\) (complete-data MLE variance) and \(T_\infty = I_{\text{obs}}^{-1}\) (variance recovery; Lemma L4), so \(B_\infty = I_{\text{obs}}^{-1}-I_{\text{com}}^{-1}\) and \[ r_\infty = B_\infty \bar U_\infty^{-1} = (I_{\text{obs}}^{-1}-I_{\text{com}}^{-1})I_{\text{com}} = I_{\text{obs}}^{-1}I_{\text{com}} - I_k = \mathrm{RIV}. \] Confirmed numerically in Listing 5.
Schafer’s fraction of missing information / EM rate matrix is \(D = I_{\text{com}}^{-1}I_{\text{mis}\mid\text{obs}}\) (Schafer 1997 eq (3.21), \(D = I_c^{-1}(\hat\theta)I_m(\hat\theta)\)), :1517. It shares \(I_{\text{mis}\mid\text{obs}}\) with RIV but uses \(I_{\text{com}}^{-1}\) on the left, not \(I_{\text{obs}}^{-1}\). They relate by \(D = (I_k+\mathrm{RIV})^{-1}\mathrm{RIV}\) (scalar: \(\gamma = r/(1+r)\)). Our headline \(\operatorname{tr}(\mathrm{RIV})=\operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{com}})-k\) lives specifically on the \(I_{\text{obs}}^{-1}\) side; a reader who substitutes \(D\) gets the wrong quantity.
The RIV’s denominator is the observed-data information \(I_{\text{obs}}=-\mathbb{E}[\partial^2\ell_{\text{obs}}]\) (D6), and one subtlety decides what it equals under MAR. The expectation must be taken unconditionally — over the joint law of \((Y,R)\), with the missingness mechanism included — because under MAR the response pattern \(R\) is not ancillary for \(\theta\) (Kenward and Molenberghs 1998, sec. 2), :89. Taking the expectation as if each pattern’s observed block were a fresh marginal sample — the naive expectation, valid only under MCAR — gives a different, biased information: “the expected information … is wrong because the expectation is taken under MCAR” (Kenward and Molenberghs 1998), :60.
For the MVN the difference is explicit. The naive expectation makes \(I_{\text{obs}}\) block-diagonal in \((\mu,\Sigma)\) — the familiar mean/covariance orthogonality — and is exactly what fisher_info_obs_mvn computes. Under MAR, selection on the observed values shifts each pattern’s observed mean off the marginal, and the unconditional information picks up a non-zero mean–covariance cross term, so “the orthogonality of mean and variance-covariance parameters is lost under the MAR mechanism,” vanishing iff the mechanism is MCAR (Kenward and Molenberghs 1998, sec. 3.2), :374 (Appendix C gives the cross term in closed form). Because the unconditional information needs the usually-unknown mechanism, the standing recommendation — which we adopt — is to estimate \(I_{\text{obs}}\) by the observed information \(-\partial^2\ell_{\text{obs}}/\partial\theta\partial\theta^\top\), which is consistent and needs no mechanism model: “provided the observed information matrix is used, conventional likelihood based frequentist inference is applicable in the MAR setting” (Kenward and Molenberghs 1998, sec. 5), :544. So there is a single RIV, \(\operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{mis}\mid\text{obs}})\), evaluated with the observed-data information; both deviance-bias terms calibrate to it (note after Equation 9).
A separate axis — score variance vs curvature — does not bite here. The expected (curvature) information \(-\mathbb{E}[\partial^2\ell]\) and the score variance \(\mathbb{E}[SS^\top]\) coincide under correct specification — the second Bartlett identity, or information equality — which is why D5 writes \(I_{\text{com}}=-\mathbb{E}[\partial^2\ell]=\mathbb{E}[S_{\text{com}}S_{\text{com}}^\top]\) as one object (Pawitan 2013, sec. 13.6: “assuming the model is correct we get \(\mathcal{J}=\mathcal{I}\), and the formula reduces to … the ‘naive’ variance”), :14699. Their inequality under misspecification is what turns AIC’s penalty \(k\) into Takeuchi’s \(\operatorname{tr}(\mathcal{J}\mathcal{I}^{-1})\) (TIC/GIC) (Pawitan 2013; Konishi and Kitagawa 1996), :14707, :90 — a regime that is not ours: the analysis model is correctly specified, so information equality holds, confirmed numerically with \(\operatorname{tr}(\hat{\mathcal{J}}\,I_{\text{obs}}^{-1})\approx k\) across cells (run_all.R phase 6, verification/cache/run_all-prod/phase6-info/). That axis is about misspecification; the MAR distinction above is about which sampling distribution the expectation is taken over, and it is the one that matters here.
Symbol caution. The proof below writes \(\mathcal{J}\) for the observed information (the negative Hessian \(-\partial^2\ell/\partial\theta\partial\theta^\top\) at the estimate), the EM/Louis convention. This is not Pawitan’s \(\mathcal{J}\) (the score variance \(\mathbb{E}[SS^\top]\)); under correct specification both reduce to \(I\), so the body’s \(\mathcal{J}\to_p I\) steps are unaffected.
We can now name the object the theorem is about, and the bias attributed to it.
D9 (MI Q-function). Multiple imputation does not maximize \(\ell_{\text{obs}}\) directly; it averages the complete-data log-likelihood over completions drawn from the analyst’s fitted imputation model. Write \(\tilde\phi\) for the imputation parameter — under congenial proper MI a draw from the imputation posterior fit to \(Y_{\text{obs}}\), centered at \(\hat\theta_{\text{obs}}\) (the symbol \(\phi\) rather than \(\theta\) marks it as the imputer’s; see Section 4). Each completion is \(\tilde Y_{\text{mis}}^{(l)} \sim P(\,\cdot \mid Y_{\text{obs}}, \tilde\phi^{(l)})\), and in the infinite-imputation limit \[ \bar Q_\infty(\theta) := \lim_{M\to\infty}\frac1M\sum_{l=1}^M \ell_{\text{com}}\big(\theta; Y_{\text{obs}}, \tilde Y_{\text{mis}}^{(l)}\big) = \mathbb{E}_{\tilde\phi}\,\mathbb{E}_{Y_{\text{mis}}\mid Y_{\text{obs}},\,\tilde\phi}\!\big[\ell_{\text{com}}(\theta)\big]. \] This is the EM/Meng Q-function \(Q(\theta\mid\hat\theta_{\text{obs}})\) in the proper-MI limit: the inner expectation conditions on the fitted \(\tilde\phi\approx\hat\theta_{\text{obs}}\), exactly as the EM E-step conditions \(Q(\theta\mid\theta^{(t)})\) on the current iterate — not on the unknown truth \(\theta_0\). Whether this mean-level approximation is enough to leave \(\bar Q_\infty\)’s bias unchanged depends on the imputation engine, not only on the mean of \(\tilde\phi\); the callout below makes this precise.
Source. Deterministic-FIML evaluation barQ_fiml_analytic\((\cdot,\hat\theta_{\text{obs}})\) (verification/00-setup.R) conditions the completion at the fitted \(\hat\theta_{\text{obs}}\); [mc] equal to this construction on the fixture (Listing 7). The fitted-vs-oracle distinction (callout below) is checked in verification/rederivation-L0-fitted-vs-oracle.R.
It is tempting to write \(\bar Q_\infty(\theta)=\mathbb{E}[\ell_{\text{com}}(\theta)\mid Y_{\text{obs}}]\), the genuine conditional expectation under the data-generating law. That is wrong, and the slip is not cosmetic: it flips the sign of the headline. The genuine conditional expectation averages \(Y_{\text{mis}}\) over the true law \(P(\cdot\mid Y_{\text{obs}},\theta_0)\); \(\bar Q_\infty\) averages over the fitted imputation law \(P(\cdot\mid Y_{\text{obs}},\tilde\phi)\) with \(\tilde\phi\approx\hat\theta_{\text{obs}}\). The two agree only if \(\tilde\phi=\theta_0\), which it is not — the analyst imputes from a model fit to the data, not from an oracle. The gap between them is the imputation bias Term A (Section 4.4), \(+\operatorname{tr}(\mathrm{RIV})\).
The distinction is decisive because \(\hat\theta_{\text{obs}}\) is a function of \(Y_{\text{obs}}\) alone. Were \(\bar Q_\infty\) the oracle conditional expectation, iterated expectation (the tower property) would give \[ \mathbb{E}\big[\,\mathbb{E}[\ell_{\text{com}}(\hat\theta_{\text{obs}})\mid Y_{\text{obs}}] - \ell_{\text{com}}(\hat\theta_{\text{obs}})\,\big] = 0, \] so Term A would vanish and the deviance bias would collapse to \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\). Term A is instead \(+\operatorname{tr}(\mathrm{RIV})\) precisely because the completions are drawn at the fitted \(\tilde\phi\), not at \(\theta_0\). This is confirmed directly: imputing from the true \(\theta_0\) drives \(\mathbb{E}[A]\to 0\) and the total to \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\), while imputing from the fitted \(\hat\theta_{\text{obs}}\) gives \(\mathbb{E}[A]=+\operatorname{tr}(\mathrm{RIV})\) and the total \(+\tfrac12\operatorname{tr}(\mathrm{RIV})\) for a second-moment-matched imputer — deterministic FIML or bootstrap resampling with an EM plug-in, the engine this paper’s own simulations use — (verification/rederivation-L0-fitted-vs-oracle.R); a genuine Bayesian posterior draw is not second-moment-matched and carries the correction stated in the callout below.
Write \(\Xi=b_\Sigma+\operatorname{Var}(\mu^\ast)\) for an imputer’s centering bias in \(\Sigma\) plus its draw covariance in \(\mu\), both evaluated against the fitted \(\hat\theta_{\text{obs}}\). For an imputer whose draw is regular — \(\hat\theta_{\text{obs}}+O_p(n^{-1/2})\) with controlled third moments, the condition every engine below satisfies, and which excludes only pathological constructions such as a draw that rarely jumps to a fixed distant value — the \(+\tfrac12\operatorname{tr}(\mathrm{RIV})\) headline holds, to leading order, whenever \(\Xi=o(1/n)\) per completed cell and the design-imbalance term below also vanishes or is controlled; call an imputer meeting both conditions second-moment-matched. Deterministic FIML meets both trivially, since it does not draw a parameter at all. Bootstrap resampling with an EM plug-in, the engine this paper’s own simulations use, meets the first: for complete data \(\mathbb{E}^\ast[\tilde\Sigma+(\tilde\mu-\hat\mu)(\tilde\mu-\hat\mu)^\top]=\hat\Sigma\) exactly, so the bootstrap’s covariance deflation cancels its mean’s draw variance, and \(\Xi\) is zero up to a residual of the observed-data estimator’s own bias under missingness; like every drawing engine, though, the bootstrap still carries the small design-imbalance term under MAR.
A genuine Bayesian posterior draw is not second-moment-matched, by construction: a proper draw is over-dispersed relative to the plug-in (Rubin’s between-imputation variance rule is built to capture that extra spread), so \(\Xi\succ0\). Such a draw carries a persistent, computable correction to the headline, \[ D=-\tfrac12\operatorname{tr}(\mathrm{RIV}_{\mu\mu})-\tfrac12\sum_i\operatorname{tr}\!\big(\hat V_i^{-1}C_i\,b_\Sigma\,C_i^\top\big)-\mathrm{IMB}+o(1). \] \(\mathrm{RIV}_{\mu\mu}\) is the mean block of \(\mathrm{RIV}\). For a unit with missing block \(M\) and observed block \(O\), \(\hat B=\hat\Sigma_{MO}\hat\Sigma_{OO}^{-1}\) and \(\hat V=\hat\Sigma_{MM}-\hat B\hat\Sigma_{OM}\) are the fitted regression map and conditional covariance of Section 3.2’s Schur complement, and \(C=[I,-\hat B]\) carries a full covariance perturbation to its effect on \(\hat V\). \(\mathrm{IMB}\) shares its design-imbalance origin with \((A)+(C)\) (Section 4.4) but arises from the imputer’s own \(\Sigma\)-draw rather than from the fitted model. It is present under MAR for every drawing engine (deterministic FIML excepted, since it never draws \(\Sigma\)) and exactly zero under MCAR; measured at \(\approx-0.006\operatorname{tr}(\mathrm{RIV})\) in this paper’s own W1 cell, not yet derived in closed form. Every other piece is computable from quantities the derivation already defines; \(b_\Sigma\) is prior-dependent, so no engine-free correction constant exists — a Jeffreys-prior draw measured here gives \(\approx-0.78\operatorname{tr}(\mathrm{RIV})\) at \(N=200\), while recentring the same draws at their own measured \(b_\Sigma\) collapses the correction to \(\approx-0.15\operatorname{tr}(\mathrm{RIV})\), matching a mean-centered Gaussian draw.
[mc] the four algebraic identities behind this reduction — the Schur-complement derivative and its concavity, the resulting cancellation, and the block-inverse identity that isolates \(\mathrm{RIV}_{\mu\mu}\) — reduce to zero symbolically (verification/cas-wolfram/verify_engine_gap_sympy.py); the numerical values above are measured directly from each engine’s parameter draws, with no imputation or log-likelihood evaluation needed (verification/_modules/engine-gap-audit.R).
\(\bar Q_\infty\) here is the averaged log-likelihood. It is not Rubin’s pooled point estimate \(\bar Q_m = \sum_l \hat Q_{*l}/m\) (Rubin 1987), :2161, which averages complete-data estimates. We reserve \(\bar Q_\infty\) for the averaged log-likelihood throughout.
D10–D12 (deviance-bias terms). The theorem compares \(\bar Q_\infty\) at \(\hat\theta_{\text{obs}}\) against the complete-data gold standard \(\ell_{\text{com}}(\hat\theta_{\text{com}})\), split at the pivot \(\hat\theta_{\text{obs}}\) into an imputation-bias term \(A\) and an estimation-mismatch term \(B\) (on the log-likelihood scale; the deviance scale doubles each): \[ A := \bar Q_\infty(\hat\theta_{\text{obs}}) - \ell_{\text{com}}(\hat\theta_{\text{obs}}) \quad(\text{imputation bias}), \] \[ B := \ell_{\text{com}}(\hat\theta_{\text{obs}}) - \ell_{\text{com}}(\hat\theta_{\text{com}}) \quad(\text{estimation mismatch}), \] \[ T := A + B = \bar Q_\infty(\hat\theta_{\text{obs}}) - \ell_{\text{com}}(\hat\theta_{\text{com}}) \quad(\text{total deviance bias}). \tag{4}\]
Source. This work; operationalized in the Monte Carlo study (verification/). The theorem (Section 4): \(\mathbb{E}[T]=+\tfrac12\operatorname{tr}(\mathrm{RIV})\), \(\mathbb{E}[A]=+\operatorname{tr}(\mathrm{RIV})\), \(\mathbb{E}[B]=-\tfrac12\operatorname{tr}(\mathrm{RIV})\).
D13 (congeniality). We assume the imputation and analysis models are congenial in Meng’s sense — loosely, the imputer’s posterior moments agree with the analyst’s procedure. This is what licenses the variance-recovery and Q-function identities the proof relies on (Lemma L4); it is stated precisely as an assumption in Section 3.
Source. Meng (1994), Definitions 1-2: posterior moments match the analyst’s procedure — \([\hat\theta(Z_o),U(Z_o)]\simeq[E_f(\theta\mid Z_o),V_f(\theta\mid Z_o)]\) [eq (2.3.1), :103] and \([\hat\theta(Z_c),U(Z_c)]\simeq[E_f(\theta\mid Z_c),V_f(\theta\mid Z_c)]\) [eq (2.3.2), :109]; Definition 2 adds \(f(Y_{\text{mis}}\mid Z_o)=g(Y_{\text{mis}}\mid Z_o,A)\) [eq (2.3.3), :121].
Because \(\theta\) carries \(\Sigma\) through \(\operatorname{vech}\Sigma\), the normal information matrices below are written with the duplication matrix relating the half- and full-vectorizations of a symmetric matrix.
D14 (half-vectorization and duplication matrix). For symmetric \(A\) (\(n\times n\)), \(\operatorname{vech}(A)\) stacks the \(\tfrac12 n(n+1)\) entries of its lower triangle, and the duplication matrix \(D_n\) satisfies \[ D_n \operatorname{vech}(A) = \operatorname{vec} A, \qquad \operatorname{vech}(A) = (D_n'D_n)^{-1}D_n'\operatorname{vec} A. \]
For the normal model the complete information has a closed form — the one the verification in Section 8 computes with.
D15 (closed-form MVN information). On \(\theta=(\mu,\operatorname{vech}\Sigma)\) the complete information is block-diagonal, \[ I_{\text{com}}(\theta) = \begin{pmatrix} N\Sigma^{-1} & 0 \\ 0 & \tfrac{N}{2}D_p'(\Sigma^{-1}\otimes\Sigma^{-1})D_p \end{pmatrix}. \]
Source. Magnus and Neudecker (2019), MVN ML section: information for \(\mu\) and \(\operatorname{vech}(\Omega)\) as \(n\Omega^{-1}\) and \(\tfrac{n}{2}D_m'(\Omega^{-1}\otimes\Omega^{-1})D_m\) [:1800-:1813] (\(\Omega\!\leftrightarrow\!\Sigma\), \(m\!\leftrightarrow\!p\), \(n\!\leftrightarrow\!N\)). \(I_{\text{obs}}\) is the pattern-mixture sum of per-pattern marginal informations (fisher_info_obs_mvn, verification/00-setup.R); [mc] against this closed form in Listing 5.
This section collects what the proof rests on but does not itself establish: the regularity, ignorability, congeniality, and properness conditions (R1–R7); the standard matrix and expectation identities; and four asymptotic lemmas (L1–L4). Each labeled item is cited by the proof in Section 4. As in Section 2, provenance — source, equation or theorem number, and line reference, with any [mc]/[derived] status — sits in the margin.
Seven conditions are in force throughout, and each is simply the hypothesis under which some later step is exact to leading order. R1–R2 are the standard likelihood-asymptotics regularity that makes the Fisher informations and the MLE expansions of Section 3.3 well-defined. R3 makes the missingness ignorable, so the mechanism drops out of the likelihood (D2). R4–R7 are the multiple-imputation conditions — congeniality, properness, self-efficiency, and the infinite-imputation idealization — that let the imputation limit recover the observed-data variance (Lemma L4) and that fix \(\operatorname{Cov}(\tilde\phi,\hat\theta_{\text{obs}})\) in Term A. None is exotic.
R1 (regular model). The complete-data model is differentiable in quadratic mean, with Lipschitz log-density and a consistent MLE; differentiation and integration interchange, so \(P_\theta\dot\ell_\theta=0\) (scores are mean-zero) and \(P_\theta\ddot\ell_\theta=-I_\theta\).
R2 (nonsingular information). Both informations are positive definite, \(I_{\text{com}},I_{\text{obs}}\succ 0\), so the inverses in RIV and in the variance formulas below exist.
Source. Vaart (1998), §5.5: nonsingular Fisher information \(I_{\text{com}},I_{\text{obs}}\succ 0\) (see Section 7.3).
R3 (ignorability). The data are missing at random, with the missingness parameters distinct from \(\theta\), so the missingness mechanism need not be modeled (D2).
R4 (congeniality). The imputation and analysis models are congenial in Meng’s sense (D13) — the imputer’s posterior moments agree with the analyst’s procedure.
R5 (properness). The imputation is second-moment proper in Rubin’s sense: its posterior draws carry the correct first two moments, so the between-imputation variance behaves as the combining rules assume.
R6 (self-efficiency). The analyst’s estimator is self-efficient: no convex combination of the observed- and complete-data estimators beats the complete-data one, i.e. no \(\lambda\) makes \(\operatorname{MSE}(\lambda\hat\theta(W_o)+(1-\lambda)\hat\theta(W_c))<\operatorname{MSE}(\hat\theta(W_c))\).
R7 (\(M\to\infty\)). Infinite-imputation idealization (\(\bar Q_\infty,T_\infty\)); finite-\(M\) corrections are Rubin’s.
Source. Rubin (1987), §4: infinite-imputation idealization (see Section 7.5).
The proof uses only standard identities, named at the point of use and reproduced verbatim with line references in Section 7:
Four facts about the estimators and the infinite-imputation limit supply the asymptotic inputs the proof draws on.
L1 (MLE asymptotic linearity and normality). Under R1-R2, \[ \sqrt{n}(\hat\theta_n-\theta_0)=\tfrac{1}{\sqrt n}\textstyle\sum_i I_{\theta_0}^{-1}\dot\ell_{\theta_0}(X_i)+o_P(1)\;\rightsquigarrow\;N(0,I_{\theta_0}^{-1}), \] hence \(\operatorname{Var}(\hat\theta_{\text{obs}})=I_{\text{obs}}^{-1}\), \(\operatorname{Var}(\hat\theta_{\text{com}})=I_{\text{com}}^{-1}\).
L3 (Fisher/Louis score identity). Differentiating the loglik decomposition \(\ell(\theta\mid Y)=\ell(\theta\mid Y_{\text{obs}})+\log P(Y_{\text{mis}}\mid Y_{\text{obs}},\theta)+c\) once and taking \(\mathbb{E}[\cdot\mid Y_{\text{obs}}]\) (conditional score mean-zero, R1) gives \(U_{\text{obs}}=\mathbb{E}[U_{\text{com}}\mid Y_{\text{obs}}]\); hence by the score-covariance identity (Section 3.2), \(\operatorname{Cov}(U_{\text{obs}},U_{\text{com}})=\operatorname{Var}(U_{\text{obs}})=I_{\text{obs}}\).
L2 (variance of the estimator gap). \[ \operatorname{Var}(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}}) = I_{\text{obs}}^{-1}-I_{\text{com}}^{-1}. \tag{5}\] Proof. By L1, \(\hat\theta_\bullet-\theta_0\approx I_\bullet^{-1}U_\bullet\); by L3, \(\operatorname{Cov}(U_{\text{obs}},U_{\text{com}})=I_{\text{obs}}\), so \(\operatorname{Cov}(\hat\theta_{\text{obs}},\hat\theta_{\text{com}})=I_{\text{obs}}^{-1}I_{\text{obs}}I_{\text{com}}^{-1}=I_{\text{com}}^{-1}\); then \(\operatorname{Var}(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}})=I_{\text{obs}}^{-1}-I_{\text{com}}^{-1}-I_{\text{com}}^{-1}+I_{\text{com}}^{-1}\). The variance is of the gap, not of \(\hat\theta_{\text{obs}}\) alone — the distinction that keeps Term B at \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\) rather than a spurious \(-\tfrac12 k\). This identity is exact only to leading order: it rests on \(\operatorname{Cov}(\hat\theta_{\text{obs}},\hat\theta_{\text{com}})=I_{\text{com}}^{-1}\) (via L3), which holds asymptotically with an \(O(n^{-1})\) remainder. Evaluated with the observed-data information \(I_{\text{obs}}\) (Section 2.4), it gives the single RIV that Term B calibrates to.
Source. [derived] from L1, L3, and the score-covariance identity (Section 3.2); verified symbolically (Tier B, Listing 11).
L4 (MI variance recovery). Under R4-R7, \(\bar\theta_\infty\) is consistent and \(T_\infty=\operatorname{Var}(\hat\theta_{\text{obs}})=I_{\text{obs}}^{-1}\). (Anchors the D8 RIV identification.)
Source. The equality \(T_\infty=\operatorname{Var}(\hat\theta_{\text{obs}})\) is the variance-recovery property of congenial, proper MI (R4): under congeniality the analyst’s procedure is the observed-data MLE, and Rubin’s proper-imputation variance is randomization-valid for the repeated-sampling variance of the complete-data statistic (Rubin 1987, sec. 4.2, conditions (4.2.1)-(4.2.2) and Result 4.1, :3489/:3563); \(\operatorname{Var}(\hat\theta_{\text{obs}})=I_{\text{obs}}^{-1}\) then holds by regularity (L1). Meng (1994) (Lemma 3 [:486]; Main result [:503/:511]) gives the efficiency dominance — MI is at least as efficient as the analyst’s incomplete-data estimator — with equality under congeniality.
Each step names the assumptions and lemmas of Section 3 that it invokes; a margin [mc] marks a step checked symbolically in Section 8. Let \(\theta_0\) denote the data-generating parameter and \(\tilde\phi\) the imputation parameter — the symbol \(\phi\) rather than \(\theta\) marks it as the imputer’s, distinct from the analysis parameter even where the two spaces coincide under congeniality, and the tilde marks a posterior draw (which collapses to the point value \(\hat\theta_{\text{obs}}\) under FIML). We write \(S\) for the score and \(\mathcal{J}\) for the observed information (the negative Hessian at the estimate) — not the score-variance \(\mathcal{J}\) of Pawitan (2013) (:14618); under correct specification both reduce to the expected information \(I\) (see Section 2.4), so the \(\mathcal{J}\to_p I\) steps below are unaffected.
Theorem (Q-function deviance bias). Under R1-R7, for a model that estimates a scale or covariance (so the conditional missing-data entropy depends on the estimated parameters — the multivariate normal with estimated \(\Sigma\) throughout; scope remark below), to leading order \[ \mathbb{E}[T] = \mathbb{E}\big[\bar Q_\infty(\hat\theta_{\text{obs}}) - \ell_{\text{com}}(\hat\theta_{\text{com}};Y)\big] = +\tfrac12\operatorname{tr}(\mathrm{RIV}) \;+\; [(A)+(C)] \;+\; O(n^{-1}), \tag{6}\] decomposing (D10-D12) as \(\mathbb{E}[A]=\operatorname{tr}(\mathrm{RIV})+[(A)+(C)]+O(n^{-1})\) (imputation bias) and \(\mathbb{E}[B]=-\tfrac12\operatorname{tr}(\mathrm{RIV})+O(n^{-1})\) (estimation mismatch). The deviance-scale statement multiplies by 2.
The term \((A)+(C)\) is a design-imbalance correction carried by Term A (Section 4.4): it vanishes under MCAR (missing and observed units share a conditioning-variable distribution), recovering the clean \(\mathbb{E}[T]=+\tfrac12\operatorname{tr}(\mathrm{RIV})\) that holds for any regular likelihood with an estimated scale or covariance (the scope remark below; the headline is not universal — a known-scale, location-only fit gives \(\mathbb{E}[T]=0\) under FIML, and a canonical-GLM family whose base measure is not constant adds the computable offset \(\Delta_h\) of Section 4.6 even under MCAR), and is a nonzero \(O(1)\) under MAR. Its general-MVN structure is derived in Section 4.4 — \((A)+(C)=\alpha^\top\mathbb{E}[\delta]+\tfrac12\operatorname{tr}(\mathcal H\,I_{\text{obs}}^{-1})-\operatorname{tr}(\mathrm{RIV})\), with the entropy gradient \(\alpha\) supported only on the covariance parameters — and for the bivariate monotone MVN it has the closed form \[ (A)+(C) = \tfrac{n_{\text{mis}}}{n_{\text{obs}}}\Big[1-\tfrac12\operatorname{tr}\!\big(Q_{\text{mis}}Q_{\text{obs}}^{-1}\big)\Big], \] \(Q_{\text{mis}},Q_{\text{obs}}\) the conditioning-variable second-moment matrices of the missing vs observed units, so \((A)+(C)=0\) iff \(Q_{\text{mis}}=Q_{\text{obs}}\) (MCAR). For the general-MVN monotone case a leading-order analytic approximation (expected information) gives the \(n\)-independent constant \(\alpha^\top\mathbb{E}[\delta]+\tfrac12\operatorname{tr}(\mathcal H\,I_{\text{obs}}^{-1})\approx-0.22\) (verification/term-a-mar-closedform.R, with the MCAR control returning \(\approx0\)). For non-monotone patterns the observed-data likelihood does not factor; its covariance-MLE bias \(b_\Sigma\) is instead derived from the general second-order (Cox–Snell) bias of the MLE — confirmed for the bivariate design by two computer-algebra systems and direct Monte-Carlo, and now carried to general \(p\) including the four-variate phase-8 ampute design (\(k{=}14\); verification/cas-wolfram/verify_term_ac_nonmonotone{,_genp}.py), where it reduces exactly to the bivariate value, meets the score/information identities to machine precision, and recovers the complete-data bias \(-\sigma\) on the always-observed block. The information convention for \((A)+(C)\) is settled; its magnitude is an \(O(1)\) that direct Monte Carlo does not pin sharply. The convention is the realized (observed) information (Kenward and Molenberghs 1998) — which under MAR carries the mean–covariance cross term, is biased otherwise, and differs from the naive expected information by the \(O(1)\) realized-information correction \(\Delta_{\mathrm{KM}}\) (\(\approx-0.55\) non-monotone, the observed-minus-naive RIV gap; derived in closed form in Section 4.4, its block structure and \(O(1)\) order independently cross-checked off-lineage, exactly \(0\) under MCAR, and itself an ultra-precise paired contrast — Monte-Carlo standard error \(<10^{-3}\)). The leading-order (\(\delta\)-method) analytic value is \(\approx-0.22\) for the four-variate design (right sign and \(O(1)\) order; MCAR control \(\approx0\)). The magnitude, by contrast, is the mean of a heavy-tailed statistic — the per-replicate remainder \(A_{\mathrm{rb}}-\operatorname{tr}(\mathrm{RIV})\) has standard deviation \(\propto\sqrt n\) and a long right tail (from over-fitted draws) — so its direct Monte-Carlo estimate is imprecise and untrended: \(-0.34,-0.46,-0.29,+0.17\) across \(n=800,1500,3000,6000\) (standard errors \(0.09,0.13,0.25,0.35\); inverse-variance pooled \(\approx-0.35\pm0.07\)), broadly consistent with the \(-0.22\) leading order and showing no real growth with \(n\) (verification/term-ac-rivgap-lavaan.R, verification/term-ac-rem-distribution.R). A variance-reduced (control-variate) estimator probes the asymptote: pairing each replicate’s \(A_{\mathrm{rb}}\) against its decoupled second-order Taylor (population coefficients) cancels the \(\sqrt n\) noise and isolates the higher-order remainder on top of that anchor, which is a positive \(O(1/n)\) (coefficient \(\approx100\)–\(130/n\); standard error \(<0.013\); a separate paired diagnostic gives the third-order piece as \(O(1/n)\) likewise), with a blind off-lineage derivation confirming the same \(O(1/n)\) remainder order. This establishes that the finite-\(n\) correction beyond the \(-0.22\) leading order is \(O(1/n)\) — so \(-0.22\) is consistent with the \(n\to\infty\) asymptote — but the result is conditional on the \(-0.22\) analytic anchor: the estimator is anchored to that value and does not independently validate it (and a finite-difference cross-check of the second-order term does not cleanly reproduce \(-0.22\), plausibly Monte-Carlo noise but undiagnosed). The absolute \((A)+(C)\) value is thus not independently pinned; what is firm is its sign, \(O(1)\) order, MCAR vanishing, and the \(O(1/n)\) correction order (verification/term-ac-controlvariate-pin.R, verification/term-ac-higherorder-paired.R). And because \((A)+(C)\) enters model comparison and likelihood-ratio numerators only through a difference across candidates — a paired contrast that cancels the common heavy-tailed noise — the comparison-relevant correction is structurally better-conditioned than the standalone value. That clause is now derived and measured, not just asserted: at the null the paired differential’s \((A)+(C)\) content cancels exactly and its noise is \(O(1)\) against the level’s \(O(\sqrt n)\) (Section 5, Propositions L3–L4; sd \(\approx 1.5\) flat in \(n\) against \(18\to37\) for the level). What is firm: \((A)+(C)\) is a nonzero \(O(1)\) under MAR and vanishes under MCAR (controls \(\approx0\) throughout).
Here \(\mathrm{RIV}=\operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{mis}\mid\text{obs}})\) (D8) is evaluated with the observed-data information \(I_{\text{obs}}\) — the unconditional expectation that, under MAR, carries the mean–covariance cross term and is estimated by the observed information (Section 2.4). Both Term A and Term B calibrate to this single RIV; the design-imbalance \((A)+(C)\) is additional to it (it is the part of Term A beyond \(\operatorname{tr}(\mathrm{RIV})\), not a re-grading of the information).
The \(+\tfrac12\operatorname{tr}(\mathrm{RIV})\) result is not universal over regular likelihoods. It requires the conditional missing-data entropy \(\mathcal{E}(\theta)=-\mathbb{E}_\theta[\log P(Y_{\text{mis}}\mid Y_{\text{obs}},\theta)]\) to depend on the estimated parameters — equivalently, a nonzero entropy gradient \(\alpha=-\mathcal{E}'\). For the multivariate normal this dependence lives entirely in the covariance (\(\alpha_\mu=0\), \(\alpha_\Sigma\neq0\); the warning in Section 4.4), so the headline is a statement about models that estimate a scale or covariance — the regime the paper instantiates and applies.
When the scale is known and only a location is estimated the result changes. Then \(\alpha=0\), so piece \((A)\) vanishes; but the cross-entropy curvature does not — for a fixed conditional covariance \(V\) a fourth-moment (Isserlis) calculation gives \(H_{\phi\phi}=-I_{\text{mis}\mid\text{obs}}\), so piece \((C)=\tfrac12\operatorname{tr}(H_{\phi\phi}I_{\text{obs}}^{-1})=-\tfrac12\operatorname{tr}(\mathrm{RIV})\) survives uncancelled, even under MCAR. The imputation bias halves, \(\mathbb{E}[A]=\operatorname{tr}(\mathrm{RIV})-\tfrac12\operatorname{tr}(\mathrm{RIV})=\tfrac12\operatorname{tr}(\mathrm{RIV})\), and with \(\mathbb{E}[B]=-\tfrac12\operatorname{tr}(\mathrm{RIV})\) the net bias collapses to \[ \mathbb{E}[T]_{\text{known scale}}=\begin{cases}0 & \text{deterministic FIML},\\[2pt]-\tfrac12\operatorname{tr}(\mathrm{RIV}) & \text{proper MI},\end{cases} \] the two arms differing by the posterior-draw imputation variance \(C_n=\tfrac12\operatorname{tr}(\mathrm{RIV})\) — so deterministic FIML and proper MI are not interchangeable for a known-scale fit. This \(C_n\) does not vanish once the scale is estimated: for any imputer whose draw is not second-moment-matched to the plug-in, it survives as the \(\mu\)-block \(\tfrac12\operatorname{tr}(\mathrm{RIV}_{\mu\mu})\) of the correction in Section 2.5, and only a second-moment-matched imputer such as bootstrap resampling with an EM plug-in cancels it.
A known-variance experiment (estimate the mean, \(\Sigma\) fixed) confirms every value to Monte-Carlo error (verification/rederivation-knownvar.R, \(R=2\times10^5\): \(\mathbb{E}[A_{\text{fiml}}]/\tfrac12\operatorname{tr}(\mathrm{RIV})=1.04\), \(\mathbb{E}[T_{\text{fiml}}]\approx0\), \(\mathbb{E}[T_{\text{proper}}]/(-\tfrac12\operatorname{tr}(\mathrm{RIV}))=0.96\), \(C_n/\tfrac12\operatorname{tr}(\mathrm{RIV})=1.000\)), against the estimated-covariance baseline \(\mathbb{E}[T]=+\tfrac12\operatorname{tr}(\mathrm{RIV})\) (verification/rederivation-L0-fitted-vs-oracle.R). The information-criteria application (Section 1.3) compares multivariate-normal models with an estimated covariance and so sits squarely in the estimated-scale regime; the known-scale collapse affects only location-only fits with a fixed dispersion, which the applications do not use. Both regimes are instantiated beyond the multivariate normal in Section 4.6: for canonical generalized linear models the estimated-parameter dependence delivers the headline up to a computable base-measure offset (\(\Delta_h\) — zero for the Bernoulli, nonzero for the Poisson), and a fixed-scale location family confirms the collapse outside the normal family.
The proof pivots on the observed-data estimator \(\hat\theta_{\text{obs}}\), splitting the total bias into the two terms of Equation 4, \[ T = \underbrace{\big[\bar Q_\infty(\hat\theta_{\text{obs}}) - \ell_{\text{com}}(\hat\theta_{\text{obs}})\big]}_{A:\ \text{imputation bias}} \;+\; \underbrace{\big[\ell_{\text{com}}(\hat\theta_{\text{obs}}) - \ell_{\text{com}}(\hat\theta_{\text{com}})\big]}_{B:\ \text{estimation mismatch}}, \] and evaluating each in expectation to leading order.
Term B is elementary. Both quantities are complete-data log-likelihoods, so their difference is a second-order Taylor expansion of \(\ell_{\text{com}}\) about its maximizer \(\hat\theta_{\text{com}}\); its expectation is \(-\tfrac12\operatorname{tr}\!\big(I_{\text{com}}\operatorname{Var}(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}})\big)\). The variance of the estimator gap is \(I_{\text{obs}}^{-1}-I_{\text{com}}^{-1}\) (the variance-of-the-gap lemma below), and the trace collapses to \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\). The single trap — substituting \(\operatorname{Var}(\hat\theta_{\text{obs}})\) for the variance of the gap — would produce a spurious \(-\tfrac12 k\); the lemma is stated precisely to forestall it.
Term A is the delicate half. The missing-data log-density is now averaged over completions drawn from the imputation model rather than the truth, and a second-order expansion in the imputation parameter \(\tilde\phi\) yields three pieces. The leading cross-term gives \(+\operatorname{tr}(\mathrm{RIV})\); the remaining two — an entropy-gradient term \((A)\) and a cross-entropy-curvature term \((C)\), each \(O(1)\) — cancel only when the missing and observed units share a conditioning-variable distribution (MCAR). Under MAR their sum is a nonzero \(O(1)\) design-imbalance term \((A)+(C)\), which we isolate as a lemma inside the Term A argument. (An earlier shortcut wrongly set one piece to zero; a second wrongly asserted the cancellation in general — both corrected below.)
Summing, \(\mathbb{E}[T] = \operatorname{tr}(\mathrm{RIV}) + [(A)+(C)] - \tfrac12\operatorname{tr}(\mathrm{RIV}) = +\tfrac12\operatorname{tr}(\mathrm{RIV}) + [(A)+(C)]\), the design-imbalance term vanishing under MCAR.
This term is elementary and exact to leading order.
B1. \(\hat\theta_{\text{com}}=\arg\max_\theta\ell_{\text{com}}(\theta;Y)\) [D3], so \(\ell_{\text{com}}(\hat\theta_{\text{com}})\ge\ell_{\text{com}}(\hat\theta_{\text{obs}})\) and \(B\le 0\).
B2. A second-order Taylor expansion of \(\ell_{\text{com}}\) about its maximizer \(\hat\theta_{\text{com}}\), whose first-order term vanishes since \(S_{\text{com}}(\hat\theta_{\text{com}})=0\) [D3]: \[ B = -\tfrac12(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}})^\top \mathcal{J}_{\text{com}}(\hat\theta_{\text{com}})(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}}) + O_p(n^{-3/2}). \]
B3. Take \(\mathbb{E}\); replace \(\mathcal{J}_{\text{com}}(\hat\theta_{\text{com}})\to_p I_{\text{com}}\) [R1, R2]; apply the quadratic-form identity \(a^\top Ba=\operatorname{tr}(Baa^\top)\); both estimators are asymptotically unbiased for \(\theta_0\), so \(\mathbb{E}[\cdot\,(\cdot)^\top]\) is the variance of the gap: \[ \mathbb{E}[B] = -\tfrac12\operatorname{tr}\!\big(I_{\text{com}}\operatorname{Var}(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}})\big)+O(n^{-3/2}). \]
B4. Substitute L2, then reduce the trace: \[ \mathbb{E}[B] = -\tfrac12\operatorname{tr}\!\big(I_{\text{com}}(I_{\text{obs}}^{-1}-I_{\text{com}}^{-1})\big)\ \text{[L2]}\ = -\tfrac12\big(\operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{com}})-k\big)\ \text{(trace algebra)}\ = -\tfrac12\operatorname{tr}(\mathrm{RIV})\ \text{[D8]}. \tag{7}\]
[mc] Listing 3 (B + tr(RIV)/2 == 0).
Which information. Term B is governed by the variance of the estimator gap (L2), \(\operatorname{Var}(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}})=I_{\text{obs}}^{-1}-I_{\text{com}}^{-1}\), where \(I_{\text{obs}}\) is the observed-data information (Section 2.4). Graded against the naive block-diagonal information the gap variance looked \(\approx 10\%\) inflated (\(\operatorname{tr}(I_{\text{com}}\widehat{\operatorname{Var}}(\text{gap}))\approx 6.0\) vs \(5.5\) at \(n=2000\), phase6-info/); graded against the observed-data information it is the single \(\operatorname{tr}(\mathrm{RIV})\), with no residual gap (note after Equation 9).
A1. Factor the complete-data log-likelihood \(\ell_{\text{com}}(\theta;Y)=\ell_{\text{obs}}(\theta;Y_{\text{obs}})+\ell_{\text{mis}\mid\text{obs}}(\theta)\), \(\ell_{\text{mis}\mid\text{obs}}(\theta):=\log P(Y_{\text{mis}}\mid Y_{\text{obs}},\theta)\). The \(\ell_{\text{obs}}\) parts cancel in \(A\) (D9, D10), leaving
\[ A = \mathbb{E}\big[\,\mathbb{E}_{Y_{\text{mis}}\mid Y_{\text{obs}},\tilde\phi}[\ell_{\text{mis}\mid\text{obs}}(\hat\theta_{\text{obs}})] - \ell_{\text{mis}\mid\text{obs}}(\hat\theta_{\text{obs}})\,\big]. \]
A2. By iterated expectation over the true completion \(Y_{\text{mis}}\sim\theta_0\), \[ A = \mathbb{E}_{Y_{\text{obs}}}\big[\,\mathbb{E}_{\tilde\phi}[\ell_{\text{mis}\mid\text{obs}}(\hat\theta_{\text{obs}})] - \mathbb{E}_{\theta_0}[\ell_{\text{mis}\mid\text{obs}}(\hat\theta_{\text{obs}})]\,\big], \] the discrepancy between averaging the missing-data log-density under the fitted \(\tilde\phi\) versus the truth \(\theta_0\). These two inner expectations would cancel were \(\tilde\phi=\theta_0\) (the oracle case, whence \(\mathbb{E}[A]=0\) by the tower property — Section 2.5); \(A\) is nonzero precisely because the imputations use \(\tilde\phi\approx\hat\theta_{\text{obs}}\neq\theta_0\).
A3. Second-order Taylor in \(\tilde\phi\) about \(\theta_0\), using the score-function identity \(\partial_\phi\mathbb{E}_\phi[g]\big|_{\theta_0}=\mathbb{E}_{\theta_0}[g\,S_{\text{mis}\mid\text{obs}}]\) [R1, differentiate under the integral] and Bartlett \(\mathbb{E}_{\theta_0}[S_{\text{mis}\mid\text{obs}}S_{\text{mis}\mid\text{obs}}^\top]=I_{\text{mis}\mid\text{obs}}\) [R1, D7], then expanding \(\ell_{\text{mis}\mid\text{obs}}(\hat\theta_{\text{obs}})\) about \(\theta_0\), gives three pieces: \[ A = \underbrace{\mathbb{E}[(\tilde\phi-\theta_0)^\top\alpha]}_{(A)} + \underbrace{\mathbb{E}[(\tilde\phi-\theta_0)^\top I_{\text{mis}\mid\text{obs}}(\hat\theta_{\text{obs}}-\theta_0)]}_{(B)\ \text{main}} + \underbrace{\tfrac12\mathbb{E}[(\tilde\phi-\theta_0)^\top H_{\phi\phi}(\tilde\phi-\theta_0)]}_{(C)}. \]
\(\alpha = \mathbb{E}_{\theta_0}[\ell_{\text{mis}\mid\text{obs}}(\theta_0)\,S_{\text{mis}\mid\text{obs}}(\theta_0)] = \operatorname{Cov}_{\theta_0}(\ell_{\text{mis}\mid\text{obs}},S_{\text{mis}\mid\text{obs}}) = -\mathcal{E}'(\theta_0)\), the gradient of the conditional entropy — generally nonzero (e.g. \(N(\mu,\sigma^2)\): \(\alpha_\mu=0\), \(\alpha_{\sigma^2}=-n_{\text{mis}}/2\sigma^2\neq 0\)). Gibbs’ inequality concerns varying the evaluation point of \(\mathbb{E}_{\theta_0}[\log p_\theta]\); \(\alpha\) varies the sampling distribution. The two are different functions — an earlier shortcut that set \(\alpha=0\) was wrong (v4.1 correction).
A4 (the \((A)+(C)\) term). Pieces \((A)\) and \((C)\) are each \(O(1)\) — \((A)\) pairs the bias of \(\tilde\phi\) with the entropy gradient, \((C)\) the curvature of the conditional cross-entropy with \(\operatorname{Var}(\tilde\phi)\), linked through the information geometry of \(P(Y_{\text{mis}}\mid Y_{\text{obs}},\theta)\). Their leading parts cancel iff the missing and observed units share a conditioning-variable distribution: \[ (A)+(C) = \begin{cases} O(n^{-1}) & \text{MCAR } (Q_{\text{mis}}=Q_{\text{obs}}),\\[2pt] \tfrac{n_{\text{mis}}}{n_{\text{obs}}}\big[1-\tfrac12\operatorname{tr}(Q_{\text{mis}}Q_{\text{obs}}^{-1})\big]+O(n^{-1}) & \text{MAR (bivariate monotone).} \end{cases} \] The MAR case is a nonzero \(O(1)\) design-imbalance correction; the mechanism, the closed-form MCAR cases, and the MAR formula are in the lemma below. So the main piece \((B)\) gives the leading \(\operatorname{tr}(\mathrm{RIV})\), plus this design-imbalance term: \(A = \operatorname{tr}\!\big(I_{\text{mis}\mid\text{obs}}\operatorname{Cov}(\tilde\phi,\hat\theta_{\text{obs}})\big)+[(A)+(C)]+O(n^{-1})\).
[derived] Mechanism + two closed-form families in the lemma below; closed forms computed in the companion MI-IC derivation (v4, Steps 6-7) and reproduced here so this document is self-contained; the bivariate-monotone \((A)+(C)\) closed form is independently re-derived by a second computer-algebra system (Wolfram) in verification/cas-wolfram/verify_term_ac.py. \(A\) and \(B\) both calibrate to the single observed-data RIV (note after Equation 9), now augmented by the \((A)+(C)\) design-imbalance term under MAR (Section 4.4 lemma).
Let \(f_0(\phi):=\mathbb{E}_{Y_{\text{mis}}\sim\phi}\big[\ell_{\text{mis}\mid\text{obs}}(\theta_0)\big]\) be the conditional cross-entropy — the expected missing-data log-density evaluated at the truth \(\theta_0\), over completions drawn from \(\phi\). With \(\delta:=\tilde\phi-\theta_0\), pieces \((A)\) and \((C)\) are exactly the first- and second-order terms of the Taylor expansion of \(f_0\) about \(\theta_0\): \[ (A)+(C) = \mathbb{E}[\delta^\top\nabla f_0(\theta_0)] + \tfrac12\mathbb{E}[\delta^\top\nabla^2 f_0(\theta_0)\delta] + O(n^{-3/2}), \] with \(\nabla f_0(\theta_0)=\alpha\) (the entropy gradient of A3) and \(\nabla^2 f_0(\theta_0)=H_{\phi\phi}\).
Curvature identity [R1]. Differentiating \(\nabla f_0(\phi)=\mathbb{E}_\phi[\ell_{\text{mis}\mid\text{obs}}(\theta_0)\,S_{\text{mis}\mid\text{obs}}(\phi)]\) once more under the integral (R1) and using \(\nabla_\phi S_{\text{mis}\mid\text{obs}}=-\mathcal{J}_{\text{mis}\mid\text{obs}}\): \[ H_{\phi\phi} = \mathbb{E}_{\theta_0}\!\big[\ell_{\text{mis}\mid\text{obs}}(\theta_0)\,(S_{\text{mis}\mid\text{obs}}S_{\text{mis}\mid\text{obs}}^\top - \mathcal{J}_{\text{mis}\mid\text{obs}})\big],\qquad \mathcal{J}_{\text{mis}\mid\text{obs}}=-\nabla^2_\phi\log P(Y_{\text{mis}}\mid Y_{\text{obs}},\phi). \]
Order counting. Under proper MI centered at \(\hat\theta_{\text{obs}}\) [R5, L1], \(\mathbb{E}[\delta]=\mathbb{E}[\hat\theta_{\text{obs}}-\theta_0]=O(n^{-1})\) (MLE bias) and \(\operatorname{Var}(\delta)=I_{\text{obs}}^{-1}=O(n^{-1})\), while \(\alpha=O(n)\) and \(H_{\phi\phi}=O(n)\) (both extensive in \(n_{\text{mis}}\)). Hence \[ (A)=\mathbb{E}[\delta]^\top\alpha = O(1),\qquad (C)=\tfrac12\operatorname{tr}\!\big(H_{\phi\phi}\,I_{\text{obs}}^{-1}\big)+O(n^{-1}) = O(1) \] are two \(O(1)\) pieces of opposite sign. \((A)\) pairs the entropy gradient \(\alpha=-\mathcal{E}'(\theta_0)\) with the MLE bias of \(\tilde\phi\) (scaling with \(1/n_{\text{obs}}\)); \((C)\) pairs the cross-entropy curvature \(H_{\phi\phi}\) with \(\operatorname{Var}(\tilde\phi)\) (scaling with the inverse design moment \(Q_{\text{obs}}^{-1}\)), over \(P(Y_{\text{mis}}\mid Y_{\text{obs}},\theta)\). Writing \(Q_{\text{obs}},Q_{\text{mis}}\) for the conditioning-variable second-moment matrices over the observed vs missing units, their leading parts cancel iff the per-unit design moments match, \(Q_{\text{mis}}/n_{\text{mis}}=Q_{\text{obs}}/n_{\text{obs}}\) — i.e. under MCAR. When they differ (MAR) the residual is the nonzero \(O(1)\) design-imbalance term \[ (A)+(C) = \tfrac{n_{\text{mis}}}{n_{\text{obs}}}-\tfrac12\operatorname{tr}\!\big(M_{\text{mis}}M_{\text{obs}}^{-1}\big) = \tfrac{n_{\text{mis}}}{n_{\text{obs}}}\big[1-\tfrac12\operatorname{tr}(Q_{\text{mis}}Q_{\text{obs}}^{-1})\big]+O(n^{-1}), \] with \(M_\bullet=\sum_{\text{units}}[1,X_2][1,X_2]^\top=n_\bullet Q_\bullet\), derived in closed form for the bivariate monotone MVN.
Closed-form checks. With no conditioning structure the family is MCAR-equivalent and \((A)+(C)\) cancels to leading order (MLE imputation \(\tilde\phi=\hat\theta_{\text{obs}}\); \(n_{\text{obs}}\) observed, \(n_{\text{mis}}\) missing):
| Family (no conditioning ⇒ MCAR) | \((A)\) | \((C)=\tfrac12\operatorname{tr}(H_{\phi\phi}I_{\text{obs}}^{-1})\) | \((A)+(C)\) |
|---|---|---|---|
| \(N(\mu,\sigma^2)\) | \(+\,n_{\text{mis}}/(2n_{\text{obs}})\) | \(-\,n_{\text{mis}}/(2n_{\text{obs}})\) | \(0\) (exact) |
| \(\text{Exp}(\lambda)\) | \(+\,n_{\text{mis}}/(n_{\text{obs}}-1)\) | \(-\,n_{\text{mis}}/n_{\text{obs}}\) | \(n_{\text{mis}}/[n_{\text{obs}}(n_{\text{obs}}-1)]=O(n^{-1})\) |
A univariate \(N(\mu,\sigma^2)\) cannot be MAR — there is no second variable to select on — so these are necessarily the cancelling cases. Under MAR with conditioning the cancellation fails by the design-imbalance amount above: the bivariate-monotone formula matches the simulation across selection strength \(b\) (verification/term-a-mar-diagnostic.R, study S1), and the 4-variate ampute MAR confirms a material \(O(1)\) residual (clean phase-8 sweep, \(R{=}40\text{k}\), _modules/term-a-mar-correction.R: non-monotone sample means \(-0.34\,(n{=}800),-0.46\,(n{=}1500)\), monotone \(-0.18,-0.27\) — \(O(1)\) in size, but heavy-tailed and imprecise, see below). A leading-order analytic approximation gives \(\approx-0.22\) (right sign and order); the reported convention is the realized (observed) information of Kenward and Molenberghs (1998) — the principled MAR choice. The realized RIV exceeds the naive one by the \(O(1)\) realized-information correction \(\Delta_{\mathrm{KM}}\) (the observed-minus-naive information gap), derived in closed form below.
General-MVN closed form (monotone). The exact Term-A integrand is, per unit \(i\) with missing block \(M_i\) and conditional moments \(m_{\theta,i},S_{\theta,i}\) under \(\theta\), \[ A = \sum_{i:\,M_i\neq\varnothing}\tfrac12\Big[(m_{0,i}-m_{\hat\theta,i})^\top S_{\hat\theta,i}^{-1}(m_{0,i}-m_{\hat\theta,i}) + \operatorname{tr}(S_{\hat\theta,i}^{-1}S_{0,i}) - |M_i|\Big], \] verified line-by-line against barQ_fiml_analytic (it is exactly the Rao-Blackwellized \(A_{\text{rb}}\)). Let \(f_0(\phi)=\sum_{i:\text{mis}}\mathbb{E}_{X_{M_i}\sim\phi}[\log p(X;\theta_0)]\) be the conditional cross-entropy — the imputation slot of \(A\) with the evaluation point held at \(\theta_0\). Expanding \(f_0\) to second order in \(\delta=\hat\theta_{\text{obs}}-\theta_0\) gives the clean two-piece form \[ (A)+(C)=\underbrace{\alpha^\top\mathbb{E}[\delta]}_{(A)}+\underbrace{\tfrac12\operatorname{tr}\!\big(\mathcal H_\phi\,I_{\text{obs}}^{-1}\big)}_{(C)},\qquad \alpha=\nabla_\phi f_0|_{\theta_0},\ \ \mathcal H_\phi=\nabla^2_\phi f_0|_{\theta_0}, \] where \(\mathcal H_\phi\) is the cross-entropy curvature (the imputation-slot Hessian); the evaluation-slot curvature of \(A\) supplies the leading \(\operatorname{tr}(\mathrm{RIV})\), so to this order no \(-\operatorname{tr}(\mathrm{RIV})\) correction enters. That eval–impute cross-curvature is, by Louis, the conditional missing information \(\mathcal C=\mathbb{E}[U_{\text{mis}\mid\text{obs}}U_{\text{mis}\mid\text{obs}}^\top]\), which under MAR is the selected missing information (the incomplete units are a selected subsample). Because the paper’s RIV uses the realized observed information — for which \(I_{\text{com}}-I_{\text{obs}}^{O}\) is exactly that selected missing information — the cancellation \(\operatorname{tr}(\mathcal C\,I_{\text{obs}}^{O,-1})=\operatorname{tr}(\mathrm{RIV})\) is exact under the realized (observed-information) convention. The naive expected-information RIV is the smaller one (Appendix C), and differs from it by the \(O(1)\) correction \(\Delta_{\mathrm{KM}}\) below; that convention difference separates the (noisy) measured naive value \(\approx+0.1\to+0.3\) from the realized \(\approx-0.3\), a shift of \(\Delta_{\mathrm{KM}}\approx-0.55\). The entropy gradient is covariance-only (\(\partial S_i/\partial\mu=0\Rightarrow\alpha_\mu=0\)), with closed form \[ \alpha_a = -\tfrac12\!\!\sum_{i:\text{mis}}\!\operatorname{tr}\!\big(S_{0,i}^{-1}\,\partial_a S_i\big),\quad \partial_a S_i=\Delta_{MM}-\Delta_{MO}B^\top-B\Delta_{OM}+B\Delta_{OO}B^\top,\ \ B=\Sigma_{MO}\Sigma_{OO}^{-1}, \] (\(\partial_a\) the symmetric perturbation of \(\Sigma\)-coordinate \(a\), \(\Delta=\partial_a\Sigma\)). It matches the finite-difference gradient to \(3\times10^{-5}\), and \(\nabla_\phi f_0=\alpha\) confirms \(\mathcal H_\phi\) is the correct curvature. At \(\theta_0\) the block-inverse identity \((\Sigma_0^{-1}\bar e)_M=0\) (with \(\bar e=(e_O;\,m_\phi-\mu_{0,M})\), \(e_O=y_O-\mu_{0,O}\)) collapses the mean part of \(\mathcal H_\phi\) to a Gauss–Newton form, giving the explicit per-unit Hessian \[ \mathcal H_\phi=\sum_{i:\text{mis}}\Big[-(\partial_\phi m_i)^\top S_{0,i}^{-1}(\partial_\phi m_i)+\tfrac12\,\partial^2_\Sigma\operatorname{tr}\!\big(S_{0,i}^{-1}S_{\cdot,i}\big)\Big], \] with the conditional-mean Jacobian \(\partial_\phi m_i\) (\(\mu\)-block \(I_M,-B\) deterministic; \(\Sigma\)-block \(G_a\Sigma_{OO}^{-1}e_O\), \(G_a=\Delta_{MO}-B\Delta_{OO}\)) and \(\partial^2_\Sigma\operatorname{tr}\) the Schur second differential. Taking expectations, every block reduces to the per-pattern selected moments of \(e_O\) (\(\mathbb E[e_O\mid P]\), \(\mathbb E[e_Oe_O^\top\mid P]\)), so \(\mathcal H_\phi\) — like \(\alpha\) and \(b_\Sigma\) — is analytic (no finite differences). It is validated against a finite-difference Hessian (deterministic \(\mu\)-block to \(0.2\%\); \((C)\) to \(3\times10^{-3}\)).
The one nonstandard ingredient is \(\mathbb{E}[\delta]\), the second-order FIML MLE bias of \(\hat\Sigma\) under MAR; only its \(\Sigma\)-block enters, since \(\alpha\) is covariance-only. For a monotone pattern the observed-data likelihood factors into a sequence of complete-data regressions of block \(j\) on blocks \(1{:}j{-}1\), each fit on the \(r_j=f_j n\) units observed through step \(j\) (Anderson 2003; Little and Rubin 2019; Schafer 1997). Because each step’s regression slope is conditionally unbiased given its predictors (\(\operatorname{Cov}(\hat\beta_j,\hat\Sigma_{<j})=0\)), the per-step biases combine without cross-coupling into a closed recursion for \(b_\Sigma:=\lim_n n\,\mathbb{E}[\hat\Sigma_{\text{obs}}-\Sigma_0]\): \[ b_{11}=-\sigma_{11},\quad b_{j,P}=\beta_j\,b_{P,P},\quad b_{jj}=-\frac{j}{f_j}\sigma_{jj\cdot<j}+\beta_j\,b_{P,P}\,\beta_j^\top+\frac{\sigma_{jj\cdot<j}}{f_j}\operatorname{tr}\!\big(\Sigma_{PP}\,(\Sigma^{(j)}_{PP})^{-1}\big), \] \(P=1{:}j{-}1\), \(\beta_j=\Sigma_{j,P}\Sigma_{PP}^{-1}\), \(\sigma_{jj\cdot<j}\) the step-\(j\) residual variance, and \(\Sigma^{(j)}_{PP}=\operatorname{Cov}(X_P\mid\text{observed through }j)\) the selected (design-imbalanced) predictor covariance. The first two terms are the classical regression-residual MLE bias and its recursive propagation through the sweep; the third is the design imbalance, which collapses to \((j{-}1)\sigma_{jj\cdot<j}/f_j\) — and the whole correction to zero — under MCAR, where \(\Sigma^{(j)}_{PP}=\Sigma_{PP}\). (The exact small-sample bias of \(\hat\Sigma\) under monotone missingness is classical (Hyodo et al. 2013; Kanda and Fujikoshi 1998); this recursion is its \(n\to\infty\) coefficient for the FIML estimator.)
Assembling the three ingredients, the leading-order analytic approximation to \((A)+(C)\) is the \(n\)-independent constant \(\alpha^\top b_\Sigma+\tfrac12\operatorname{tr}(\mathcal H_\phi I_{\text{obs}}^{-1})\). For the four-variate monotone-MAR design (Section 3.2) it evaluates to \[ (A)+(C)_{\text{approx}}=\underbrace{+0.54}_{\alpha^\top b_\Sigma}\ \underbrace{-0.76}_{\frac12\operatorname{tr}(\mathcal H_\phi I_{\text{obs}}^{-1})}\ \approx\ -0.22, \] with the MCAR control returning \(\approx0\) (verification/term-a-mar-closedform.R) — the correct sign and \(O(1)\) order, and (per the convention discussion below) consistent with the heavy-tailed direct estimate to within its Monte-Carlo error.
Value: settled to the observed-information convention. Which information enters the RIV is not a free choice. Under MAR the principled object is the observed (realized) information: the naive expected information requires correctly specifying the missingness mechanism and is otherwise biased, the classical Kenward–Molenberghs recommendation — “provided the observed information matrix is used, conventional likelihood based frequentist inference is applicable in the MAR setting” (Kenward and Molenberghs 1998, sec. 5), :544. The paper’s RIV (Section 2.4, Appendix C) and the empirically measured \((A)+(C)\) both use it, so the reported value is the realized one — an \(O(1)\) whose direct Monte-Carlo estimate is heavy-tailed and imprecise: the per-replicate remainder has standard deviation \(\propto\sqrt n\) and a long right tail, so its sample mean is untrended and only loosely determined (non-monotone \(-0.34,-0.46,-0.29,+0.17\) at \(n=800,1500,3000,6000\), inverse-variance pooled \(\approx-0.35\pm0.07\); monotone \(-0.18,-0.27\)), broadly consistent with the \(-0.22\) leading order. What is sharply pinned — because it is a paired contrast that cancels that sampling noise — is the convention gap \(\Delta_{\mathrm{KM}}\) below.
How much the convention matters is itself an \(O(1)\) quantity — the realized-information correction \[ \Delta_{\mathrm{KM}}=\operatorname{tr}\!\big(I_{\text{obs}}^{N,-1}I_{\text{com}}\big)-\operatorname{tr}\!\big(\mathbb{E}[I_{\text{obs}}^{O}]^{-1}I_{\text{com}}\big)=\operatorname{tr}\!\big(I_{\text{obs}}^{N,-1}\,\big(\mathbb{E}[I_{\text{obs}}^{O}]-I_{\text{obs}}^{N}\big)\,\mathbb{E}[I_{\text{obs}}^{O}]^{-1}I_{\text{com}}\big), \] the RIV-trace gap between the naive information \(I_{\text{obs}}^{N}\) (block-diagonal, population \(\Sigma_{OO}\)) and the realized \(I_{\text{obs}}^{O}\). By Appendix C, \(\mathbb{E}[I_{\text{obs}}^{O}]-I_{\text{obs}}^{N}\) has a zero \(\mu\mu\) block, a \(\mu\times\operatorname{vech}\Sigma\) cross block carried by the per-pattern selection shift \(\mathbb{E}[e_O\mid P]\) (the K&M term that the block-diagonal \(I_{\text{obs}}^{N}\) sets to zero), and a \(\operatorname{vech}\Sigma\) block carried by the selected-minus-population second moment \(\mathbb{E}[e_Oe_O^\top\mid P]-\Sigma_{OO}\) — both nonzero under MAR and exactly zero under MCAR. So \(\Delta_{\mathrm{KM}}\) is \(O(1)\) under MAR and \(0\) under MCAR. Evaluated in closed form from the selected moments (no fitting), it is \(n\)-stable at \(\approx-0.55\) (non-monotone) and \(\approx-0.44\) (monotone), with the MCAR control \(\to0\) (verification/term-ac-realized-info.R); a paired lavaan-FIML measurement at the fit returns the same gap, approaching \(\approx-0.55\) as \(n\) grows (\(-0.627,-0.582,-0.569\) at \(n=800,1500,3000\); the residual at \(-0.569\) being the \(\hat\theta\)-vs-\(\theta_0\) offset; verification/term-ac-rivgap-lavaan.R).
Switching from the (biased) naive to the (principled) observed information therefore shifts \((A)+(C)\) by exactly \(\Delta_{\mathrm{KM}}\), \[ (A)+(C)_{\text{realized}}=(A)+(C)_{\text{naive}}+\Delta_{\mathrm{KM}}, \] where \((A)+(C)_{\text{naive}}=\mathbb{E}[A_{\text{rb}}]-\operatorname{tr}(\mathrm{RIV}^{N})\approx+0.1\to+0.3\) is the expected-information-convention value — the anchored, higher-variance leg (the low-variance carriers are the paired gap and the analytic \(\Delta_{\mathrm{KM}}\)). The realized value is an \(O(1)\approx-0.3\) (heavy-tailed and imprecise — see above); naive and realized are tied by \(\Delta_{\mathrm{KM}}\) as a near-identity (they share \(A_{\text{rb}}\), so \(\mathrm{rem}_{\text{real}}\equiv\mathrm{rem}_{\text{exp}}+\mathrm{gap}\)), so the genuine precise check is the analytic \(\Delta_{\mathrm{KM}}\approx-0.55\) against the measured paired gap \(\approx-0.56\) (se \(<10^{-3}\)). What is firm: \((A)+(C)\) is a nonzero \(O(1)\) under MAR, vanishes under MCAR, and its realized−naive convention gap \(\Delta_{\mathrm{KM}}\approx-0.55\) is pinned in closed form under the principled observed-information convention.
Scope. The closed-form recursion above is for monotone MVN missingness, where the observed-data likelihood factors into complete-data regressions. For non-monotone patterns it does not factor, so \(b_\Sigma\) has no such recursion; it instead follows from the general second-order (Cox–Snell) bias of the maximum-likelihood estimate (Cox and Snell 1968; Cordeiro and Klein 1994), :269, evaluated from the first three derivatives of the (non-factoring) observed-data log-likelihood — the multivariate-normal recipe of Patriota and Lemonte (2009), adapted to the pattern-mixed observed-data score. For the minimal bivariate non-monotone design — complete / \(X_1\)-missing / \(X_2\)-missing units under an additive half-sigmoid MAR mechanism — this \(b_\Sigma\) is now derived and independently confirmed by two computer-algebra systems that agree to \(\sim10^{-9}\) (Wolfram adaptive quadrature and SymPy Gauss–Hermite, verification/cas-wolfram/verify_term_ac_nonmonotone{,_sympy}.py). It satisfies the score (\(\mathbb{E}[U]=0\)) and information (\(\mathbb{E}[UU^\top]=\mathcal I\)) identities to machine precision, reduces exactly to the monotone closed form when the \(X_2\)-missing group is removed (\(b_{\sigma_{22}}=-\sigma_{22}\), \(b_{\sigma_{12}}=-\sigma_{12}\), and the regression-residual recursion for \(b_{\sigma_{11}}\)), and matches a direct Monte-Carlo measurement of \(n\,\mathbb{E}[\hat\Sigma-\Sigma]\) in every component within Monte-Carlo error (e.g. \(b_{\sigma_{11}}=-1.19\) analytic vs \(-1.21\pm0.16\) at \(n{=}800\)). Assembling \((A)+(C)=\alpha^\top b_\Sigma+\tfrac12\operatorname{tr}(\mathcal H_\phi
I_{\text{obs}}^{-1})\) with this \(b_\Sigma\) then gives a nonzero \(O(1)\) design-imbalance that grows from \(\approx0\) under MCAR to \(\approx-0.10\) as selection strengthens, with the direct \(A_{\text{rb}}-\operatorname{tr}(\mathrm{RIV})\) remainder corroborating its sign and order (\(\approx-0.07\), the leading-order-vs-realized-information gap matching the monotone case) and the MCAR control returning \(\approx0\) (verification/verify_term_ac_nonmonotone.R). The same Cox–Snell \(b_\Sigma\) is now carried to general \(p\), including the four-variate phase-8 ampute design (\(k{=}14\); verification/cas-wolfram/verify_term_ac_nonmonotone_genp.py), where it reduces exactly to the bivariate value and passes the score/information identities. The four-variate assembly’s leading-order approximation gives \(\approx-0.22\) (verification/verify_term_ac_nonmonotone_4v.R; right sign and order); the reported convention is the realized (observed) information of Kenward and Molenberghs (1998), and the value is an \(O(1)\approx-0.3\) whose direct estimate is heavy-tailed and imprecise (Section 4.4). The realized RIV exceeds the naive one by the \(O(1)\) realized-information correction \(\Delta_{\mathrm{KM}}\approx-0.55\) (non-monotone, the observed-minus-naive information gap), the selected-vs-population information effect under MAR, derived in closed form (Section 4.4; its block structure and \(O(1)\) order independently cross-checked off-lineage) and itself sharply pinned as a paired contrast; it is not a \(b_\Sigma\) shortfall (\(b_\Sigma\) is derived and nearly selection-invariant). An earlier version asserted \((A)+(C)=O(n^{-1})\) “for any regular conditional family”; that holds only under MCAR and is retracted for MAR. Independent leading-order support that the main piece gives \(\operatorname{tr}(\mathrm{RIV})\): the FIML \(A\) matches the R code to machine precision (Listing 7, Listing 8), and the high-\(R\) phase-6 diagnostic graded against the observed-data RIV (Section 2.4) shows \(A\) and \(B\) calibrate to the single \(\operatorname{tr}(\mathrm{RIV})\) once the \((A)+(C)\) residual is accounted for (note after Equation 9).
[derived] Realized-information correction \(\Delta_{\mathrm{KM}}\) in closed form: the analytic derivation from the per-pattern selected moments (verification/term-ac-realized-info.R; the analytic \(\mathbb{E}[I_{\text{obs}}^{O}]\) reproduces the realized Hessian to \(10^{-11}\)). A blind, off-lineage cross-check by a second model family (GPT-5.5) independently reproduces the \(\Delta_{\mathrm{KM}}\) block structure and the \(O(1)\)-under-MAR / \(0\)-under-MCAR order — the latter correcting a prior cross-model blind spot (an erroneous \(O(n^{-1})\) call). The trace functional and the choice of reference were supplied in the prompt, not produced blind (todo/019). Transparent provenance per the authorship model.
A5. The cross-term needs the imputer’s draws centered at \(\hat\theta_{\text{obs}}\) with the observed-data covariance. Under congeniality [R4] the imputer’s posterior is the Bayesian posterior for \(\theta\mid Y_{\text{obs}}\), which by Bernstein–von Mises under regularity [R1; Vaart (1998), §10.2] is asymptotically \(N(\hat\theta_{\text{obs}},I_{\text{obs}}^{-1})\). This supplies the posterior alignment \(\mathbb{E}[\tilde\phi\mid Y_{\text{obs}}]=\hat\theta_{\text{obs}}+o_P(n^{-1/2})\) and, writing \(\tilde\phi=\hat\theta_{\text{obs}}+\Delta\tilde\phi\), the posterior-draw orthogonality \(\operatorname{Cov}(\Delta\tilde\phi,\hat\theta_{\text{obs}})=0\), so \(\operatorname{Cov}(\tilde\phi,\hat\theta_{\text{obs}})=\operatorname{Var}(\hat\theta_{\text{obs}})=I_{\text{obs}}^{-1}\) [L1]. Properness [R5] alone does not fix this centering — it is the congeniality/BvM alignment that does; likewise the matching \(T_\infty=I_{\text{obs}}^{-1}\) (L4) rests on self-efficiency [R6], not R5. Hence, carrying the design-imbalance term of A4, \[ \mathbb{E}[A] = \operatorname{tr}\!\big(I_{\text{mis}\mid\text{obs}}\,I_{\text{obs}}^{-1}\big) + [(A)+(C)] + O(n^{-1}) = \operatorname{tr}(\mathrm{RIV}) + [(A)+(C)] + O(n^{-1})\ \text{[D7, D8]}, \tag{8}\] with \((A)+(C)=0\) under MCAR (so \(\mathbb{E}[A]=\operatorname{tr}(\mathrm{RIV})\)) and the nonzero \(O(1)\) design-imbalance term of A4 under MAR. (The MCAR cancellation here is between two nonzero pieces and so uses \(\alpha\neq0\) — an estimated scale; for a known-scale, location-only fit piece \((C)\) instead survives at \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\) even under MCAR, halving \(\mathbb{E}[A]\). See the scope remark after the theorem.)
[mc] Listing 2 (A - tr(RIV) == 0) for the leading-order trace identity; the \((A)+(C)\) design-imbalance term is separate. The non-monotone \(b_\Sigma\) is derived via Cox–Snell, confirmed by two computer-algebra systems plus direct Monte-Carlo (verification/cas-wolfram/verify_term_ac_nonmonotone{,_sympy}.py, verification/verify_term_ac_nonmonotone.R) and extended to general \(p\) (..._genp.py, ..._4v.R). The \((A)+(C)\) information convention is settled to the observed information (Kenward and Molenberghs 1998): a leading-order approximation gives \(\approx-0.22\) (right sign and order); the magnitude is a heavy-tailed, imprecise \(O(1)\approx-0.3\) (Section 4.4). What is sharply pinned is the realized−naive information gap \(\Delta_{\mathrm{KM}}\approx-0.55\) (non-monotone), derived two ways (Section 4.4). Firm: \((A)+(C)\) is \(O(1)\) under MAR, \(\approx0\) under MCAR. The MVN FIML instantiation of \(\bar Q_\infty\) matches the R code to machine precision (Listing 7).
Which information. Term A uses \(\operatorname{Var}(\hat\theta_{\text{obs}})=I_{\text{obs}}^{-1}\) (L1), with \(I_{\text{obs}}\) the observed-data information; the high-\(R\) diagnostic confirms \(\mathbb{E}[A]\) tracks \(\operatorname{tr}(\mathrm{RIV})\) on this same observed-data information up to the \((A)+(C)\) design-imbalance term (zero under MCAR, \(O(1)\) under MAR) — the single RIV being the one Term B also calibrates to (Section 2.4; note after Equation 9).
\[ \mathbb{E}[T] = \mathbb{E}[A]+\mathbb{E}[B] = \operatorname{tr}(\mathrm{RIV}) + [(A)+(C)] - \tfrac12\operatorname{tr}(\mathrm{RIV}) = +\tfrac12\operatorname{tr}(\mathrm{RIV}) + [(A)+(C)] + O(n^{-1}) \tag{9}\] [D10-D12], combining the Term A (Equation 8) and Term B (Equation 7) results. Under MCAR \((A)+(C)=0\) and \(\mathbb{E}[T]=+\tfrac12\operatorname{tr}(\mathrm{RIV})\); under MAR the design-imbalance term remains — a nonzero \(O(1)\) whose structure is derived (entropy gradient, cross-entropy curvature, and the now general-\(p\) Cox–Snell \(b_\Sigma\)) and whose information convention is settled to the realized (observed) information (Kenward and Molenberghs 1998). Its leading-order value is \(\approx-0.22\); the magnitude is the mean of a heavy-tailed statistic, not sharply pinnable by direct Monte Carlo (\(\approx-0.3\), untrended; Section 4.4). What is precise is the realized-vs-naive information gap \(\Delta_{\mathrm{KM}}\approx-0.55\), derived in closed form (Section 4.4), its block structure and \(O(1)\) order independently cross-checked by a blind off-lineage model. \(\qquad\blacksquare\)
[mc] Listing 4 (T - tr(RIV)/2 == 0), the leading-order identity (MCAR, or before adding the \((A)+(C)\) term).
The same total has a second, more compact reading that ties the pieces together. Let \(C_n(\theta)=\sum_i\mathbb{E}_\theta[\log P(Y_{\text{mis},i}\mid Y_{\text{obs},i},\theta)\mid Y_{\text{obs},i}]\) be the (negative) conditional entropy of the missing block given the observed block; for the multivariate normal \(C_n(\theta)=-\tfrac12\sum_i[\,|M_i|\log 2\pi e+\log|\Sigma_{M\mid O,i}(\theta)|\,]\), a function of \(\Sigma\) alone. Evaluating this entropy at the estimate rather than the truth incurs the plug-in bias \(\Delta_n=\mathbb{E}[C_n(\hat\theta_{\text{obs}})-C_n(\theta_0)]\), and to leading order \[ \mathbb{E}[T]=\Delta_n=+\tfrac12\operatorname{tr}(\mathrm{RIV})+[(A)+(C)]+O(n^{-1}), \] so the conditional-entropy story and the \((A)/(B\text{-main})/(C)\) split of Section 4.4 are the same bias, regrouped. The bridge is a single identity for the entropy curvature at the truth, \[ \nabla^2 C_n(\theta_0)=\underbrace{H_{\phi\phi}}_{\text{cross-entropy curvature}}+\underbrace{I_{\text{mis}\mid\text{obs}}}_{\text{conditional information}}, \] read off by the chain rule from \(C_n(\theta)=F(\theta,\theta)\), where \(F(\phi,\psi)=\mathbb{E}_{Y_{\text{mis}}\sim\phi}[\log P(Y_{\text{mis}}\mid Y_{\text{obs}},\psi)]\) has sampling-slot Hessian \(H_{\phi\phi}\), evaluation-slot Hessian \(-I_{\text{mis}\mid\text{obs}}\), and mixed block \(+I_{\text{mis}\mid\text{obs}}\) (Bartlett). Taylor-expanding \(\Delta_n\) in \(\hat\theta_{\text{obs}}-\theta_0\) then sends its first-order part to \((A)=\alpha^\top\mathbb{E}[\delta]\) (since \(\nabla C_n(\theta_0)=\alpha\)) and its second-order part to \(\tfrac12\operatorname{tr}(H_{\phi\phi}I_{\text{obs}}^{-1})+\tfrac12\operatorname{tr}(\mathrm{RIV})=(C)+\tfrac12\operatorname{tr}(\mathrm{RIV})\). This unifies the scope cases of Section 4.1: with an estimated covariance \(\Delta_n=\tfrac12\operatorname{tr}(\mathrm{RIV})+[(A)+(C)]\), whereas for a known scale \(C_n\) is parameter-free, \(\Delta_n=0\), and the \(\mu\)-block of the very same split is the Isserlis identity \(H_{\phi\phi}=-I_{\text{mis}\mid\text{obs}}\) used there. The reading is also unambiguous: the full normalized normal density carries no parameter-free, \(Y_{\text{mis}}\)-dependent additive constant (every \(Y\)-quadratic enters through \(\Sigma^{-1}\)), so \(T\) — and with it \(\Delta_n\) — is convention-free. We keep \(+\tfrac12\operatorname{tr}(\mathrm{RIV})\) as the headline; \(\mathbb{E}[T]=\Delta_n\) is recorded as a consistency check, confirmed symbolically by two computer-algebra systems and by direct Monte Carlo (verification/cas-wolfram/verify_term_ac{,_sympy}.py, verification/rederivation-L0-fitted-vs-oracle.R).
This callout concerns which information defines the RIV — a separate matter from the \((A)+(C)\) term (Section 4.4), which is an additional Term-A correction, not a re-grading of the RIV. Setting \((A)+(C)\) aside (it is \(0\) under MCAR), the leading constants \(\mathbb{E}[A]=+\operatorname{tr}(\mathrm{RIV})\) and \(\mathbb{E}[B]=-\tfrac12\operatorname{tr}(\mathrm{RIV})\) hold to \(O(n^{-1})\) against one RIV — the observed-data RIV of Section 2.4, with \(I_{\text{obs}}\) the unconditional information (estimated by the observed information). Term A is governed by \(\operatorname{Var}(\hat\theta_{\text{obs}})=I_{\text{obs}}^{-1}\) (L1) and Term B by \(\operatorname{Var}(\hat\theta_{\text{obs}}-\hat\theta_{\text{com}})=I_{\text{obs}}^{-1}-I_{\text{com}}^{-1}\) (L2); both reduce to the same \(\operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{mis}\mid\text{obs}})\).
A high-\(R\) diagnostic confirms this directly: graded against the observed-data RIV, both \(\mathbb{E}[A]-\operatorname{tr}(\mathrm{RIV})\) and \(\mathbb{E}[B]+\tfrac12\operatorname{tr}(\mathrm{RIV})\) are small and consistent with \(O(n^{-1})\), and the asymptotic inputs L1-L3 hold against that same information (verification/lemma-revalidation.R; also run_all.R phase 6; todo/008). The earlier reading — that the two terms attach to different (expected vs sample) RIVs with a slowly closing gap — was an artifact of grading the bias against the naive, block-diagonal information (fisher_info_obs_mvn), which is the MCAR form and omits the MAR cross term (Section 2.4, Appendix C; (Kenward and Molenberghs 1998)). Against the correct observed-data information there is a single RIV and no residual gap.
Three layers of evidence: this proof (the algebra), the inline symbolic verification (Section 8), and the pre-registered + diagnostic Monte Carlo study (verification/, including the phase-6 information diagnostic).
Robustness to the imputation engine. The sample RIV in that diagnostic is built from the imputations actually drawn — here by the EMB sampler, expectation–maximization run on bootstrap resamples (Honaker, King, and Blackwell 2011), whose imputation posterior is only an asymptotic approximation. The exact alternative is data augmentation, a Markov chain Monte Carlo scheme whose draws follow the observed-data posterior at any sample size (Tanner and Wong 1987; Schafer 1997), but at much higher computational cost. On the most demanding cells — the non-monotone MAR design with as few as \(N=200\) observations, where the bootstrap approximation is weakest — the sample RIV reproduces the observed-data RIV to within about \(2\)–\(3\%\), with no widening of the gap as \(N\) grows, and the information-criterion model selection agrees between full-information maximum likelihood and the bootstrap-imputed analysis to within \(0.011\). The proper-imputation results reported here are therefore insensitive to this choice of sampler for the designs studied (run_all.R phases 5–6, verification/cache/run_all-prod/).
This robustness is empirical, and we scope it as such. The bootstrap sampler is only approximately Bayesian: it propagates parameter uncertainty by resampling rather than by drawing from the exact observed-data posterior, and such resampling rests on large-sample arguments whose small-sample behaviour is, in Little and Rubin’s assessment, less assured than the Bayesian alternative (Little and Rubin 2019, sec. 5.5). Proper data augmentation (Tanner and Wong 1987; Schafer 1997) removes that approximation at the cost of a Markov chain. We observe no resulting discrepancy on the designs examined — including the small-\(N\), non-monotone cells where the approximation is most exposed — but this does not establish engine-independence beyond them, and a fully Bayesian sampler remains the conservative choice when the fraction of missing information is large.
The one step of the proof whose general form is established above only for the multivariate normal is the A4 cancellation. For the MVN the entropy gradient \(\alpha\) is supported on the covariance coordinates alone, and the lemma’s closed-form checks cover scalar families with no covariates; for a regression family \(\alpha\) has support on every coordinate of the estimated parameter, so the MCAR cancellation needed its own proof. This subsection supplies that proof for one-parameter exponential-family responses with a canonical link, identifies the term that survives when the family’s base measure is not constant, and records the corollaries that a set of pre-registered scans then tested blind. Everything here concerns response-only missingness under MCAR; the MAR design-imbalance term of A4 is a separate object and is not re-derived for this family.
Setup. The response follows a one-parameter exponential family with canonical link, \[ f(y\mid x;\beta)=h(y)\exp\{y\eta-A(\eta)\},\qquad \eta=x^\top\beta,\qquad \beta\in\mathbb{R}^q, \] with covariates \(X\) fully observed and the response missing completely at random at rate \(r\). The estimated parameter is \(\beta\) alone (any dispersion is fixed by the family), responses are conditionally independent given \(X\), and under response-only missingness \(\ell_{\text{mis}\mid\text{obs}}(\theta)=\sum_{i\in\text{mis}}\log f(y_i\mid x_i;\theta)\) and \(\hat\theta_{\text{obs}}=\hat\beta_{\text{obs}}\) is the complete-case MLE. Write \(w=A''(\eta)\) (the information weight), \(A'''\) (the third cumulant), \(g(\eta)=\mathbb{E}_\eta[\log h(Y)]\) (the base-measure mean), \(J_{\text{obs}}=\sum_{j\in\text{obs}}w_j x_j x_j^\top\) (the realized observed-block information), and \(\Omega=\mathbb{E}[wXX^\top]\) (its per-unit population version), so that \(I_{\text{com}}=n\Omega\), \(I_{\text{obs}}=(1-r)\,I_{\text{com}}\), and \(\operatorname{tr}(\mathrm{RIV})=q\,r/(1-r)\) exactly.
Two facts transfer from the main proof unchanged. First, the deterministic-FIML evaluation is exact for any family, not only the normal: a missing row’s contribution to \(\bar Q_\infty(\hat\theta_{\text{obs}})\) is \(\mathbb{E}_{\tilde Y\sim f(\cdot;\hat\theta_{\text{obs}})}[\log f(\tilde Y;\hat\theta_{\text{obs}})]=-H_i(\hat\theta_{\text{obs}})\), the response entropy at the fit, so the analytic anchor the scans evaluate is the \(M=\infty\) FITTED limit with no Monte-Carlo gap. Second, the A/B split and the A3 expansion are family-generic under R1–R2: Term B is \(-\tfrac12\operatorname{tr}(\mathrm{RIV})+O(n^{-1})\) with no normal input, and Term A is \(\operatorname{tr}(\mathrm{RIV})+[(A)+(C)]+O(n^{-1})\) with the A4 objects instantiated for this family. The conditional cross-entropy per missing unit is \(f_{0,i}(\phi)=g(\eta_i(\phi))+A'(\eta_i(\phi))\,\eta_{0,i}-A(\eta_{0,i})\) with \(\eta_i(\phi)=x_i^\top\phi\); differentiating twice in \(\phi\) at \(\theta_0\) gives \[ \alpha=\sum_{i\in\text{mis}}\big(g'_i+\eta_i w_i\big)x_i,\qquad \mathcal H_\phi=\sum_{i\in\text{mis}}\big(g''_i+\eta_i A'''_i\big)x_i x_i^\top,\qquad (A)=\alpha^\top\mathbb{E}[\delta],\quad (C)=\tfrac12\operatorname{tr}\!\big(\mathcal H_\phi J_{\text{obs}}^{-1}\big), \] with \(J_{\text{obs}}^{-1}\) and \(I_{\text{obs}}^{-1}\) interchangeable at this order under MCAR (\(J_{\text{obs}}=(1-r)n\,\Omega\,(1+o_P(1))\)). One bookkeeping convention is stated explicitly: \(\alpha\) and \(\mathcal H_\phi\) are built from the realized design points and the realized missing set, so the \((A)\) and \((C)\) displays here — and the population-design forms below — carry an outer expectation over the covariate design and the missingness pattern as well as over the response.
The MLE-bias input. Piece \((A)\) pairs \(\alpha\) with \(b:=\mathbb{E}[\hat\beta_{\text{obs}}-\beta_0]\), the second-order bias of the complete-case MLE (for the FITTED arm \(\tilde\phi=\hat\theta_{\text{obs}}\), so \(\mathbb{E}[\delta]=b\); a proper draw centered there adds no further mean). For a canonical link the observed information carries no \(Y\), so \(J_{\text{obs}}\) is deterministic given the design, and expanding the score equation \(0=U(\hat\beta)\) to second order and taking expectations gives \[ b=-\tfrac12\,J_{\text{obs}}^{-1}\sum_{j\in\text{obs}}A'''_j\,\big(x_j^\top J_{\text{obs}}^{-1}x_j\big)\,x_j+O(n^{-2}), \] the canonical-link case of the general second-order MLE bias of Cox and Snell (1968, sec. 3), :269 — the same classical input the non-monotone \(b_\Sigma\) of the A4 lemma already draws on, here in its original regression form (their §§8–9 work it for a binomial logit design, :748). On the intercept-only Bernoulli design the formula reduces to the delta-method value \(\mathbb{E}[\operatorname{logit}\bar y]-\eta=\sinh(\eta)/n+O(n^{-2})\) — a bias directed away from zero, the direction Firth (1993) records for the logistic MLE (:211).
Source. Cox and Snell (1968), §3: the second-order score equation [:180], scalar bias \(b=(K+2J)/(2I^2)\) [:229], vector case, their eq. (20) [:269]. Direction: Firth (1993), §3.3 — the logistic-regression MLE is biased away from \(\beta=0\) [:211]. [mc] the canonical-link formula reproduces the intercept-only \(\sinh(\eta)/n\) delta-method value symbolically and the exact enumerated \(\mathbb{E}[\hat\beta]\) on a two-point design (verification/cas-wolfram/verify_bernoulli_fitted_gap.py, checks 3 and 5).
Lemma P (projection cancellation). [derived] Substituting \(b\) into \((A)\) and passing to the population design (\(J_{\text{obs}}\to(1-r)n\,\Omega\) and \(\sum_{i\in\text{mis}}\to rn\,\mathbb{E}[\cdot]\) — the step where MCAR enters), the two \(O(1)\) pieces built from the canonical-structure parts of \(\alpha\) and \(\mathcal H_\phi\) are \[ (A)_{\eta w}=-\tfrac{r}{2(1-r)}\,\mathbb{E}[\eta wX^\top]\,\Omega^{-1}\,\mathbb{E}\big[A'''\,(X^\top\Omega^{-1}X)\,X\big],\qquad
(C)_{\eta A'''}=+\tfrac{r}{2(1-r)}\,\mathbb{E}\big[\eta A'''\,(X^\top\Omega^{-1}X)\big], \] and these cancel exactly, for every design distribution and every canonical family, because \(\eta=X^\top\beta_0\) is linear in \(X\) and \(\Omega=\mathbb{E}[wXX^\top]\) uses the same weight \(w\): \(\mathbb{E}[\eta wX^\top]=\beta_0^\top\Omega\), so \[ \mathbb{E}[\eta wX^\top]\,\Omega^{-1}\,\mathbb{E}\big[A'''(X^\top\Omega^{-1}X)X\big]=\beta_0^\top\,\mathbb{E}\big[A'''(X^\top\Omega^{-1}X)X\big]=\mathbb{E}\big[\eta A'''(X^\top\Omega^{-1}X)\big]. \tag{10}\] What survives is the base-measure part alone: \[ (A)+(C)=\Delta_h+O(n^{-1}),\qquad
\Delta_h=\tfrac{r}{1-r}\Big[-\tfrac12\,\mathbb{E}[g'X^\top]\,\Omega^{-1}\,\mathbb{E}\big[A'''(X^\top\Omega^{-1}X)X\big]+\tfrac12\,\mathbb{E}\big[g''\,X^\top\Omega^{-1}X\big]\Big], \tag{11}\] which the same projection does not cancel, because \(g'(\eta)\) is not linear in \(X\).
[derived] + [mc] Lemma P (Equation 10) reduces to zero symbolically on a free-symbol three-point design and on 12 exact-rational four-point designs, and numerically on two quadrature designs (verification/cas-wolfram/verify_bernoulli_fitted_gap.py, 15 checks, all passing). An independent blind re-derivation from the setup alone (GPT-5.5 Pro, no access to this derivation or to any measured value) produced the same cancellation, the same base-measure survivor, and the same family constants (verification/cache/layer3-gpt-5.5-pro-bern-glm-modeA.md); the adversarial red-team pass returned zero refutations (...-modeB.md). Transparent provenance per the authorship model.
Result (canonical-GLM FITTED gap, response-only MCAR). Under R1–R7, \[ \mathbb{E}[T_{\text{fitted}}]=+\tfrac12\operatorname{tr}(\mathrm{RIV})+\Delta_h+O(n^{-1}). \tag{12}\] The base-measure offset \(\Delta_h\) is the refinement this family class forces on the theorem’s MCAR statement: the A4 cancellation to zero is a property of families whose base measure contributes no entropy of its own, not a consequence of MCAR alone.
The scans, registered blind. Every family constant above was derived and committed before the production run that tested it (registration commit 4aaeaae, production run committed at afa4438); the finite-\(n\) remainder constants are the stated exception, reported as measured. One disclosure per the authorship model: the derivation postdates a single \(N{=}200\) Bernoulli measurement (\(+0.846\pm0.161\) on the ratio scale at \(R{=}1000\)); its out-of-sample content is everything else — the \(N\)-scan, the Poisson family (no simulation of which existed in this repository when its target was registered), and the Student-\(t\) exemplar. On the ratio scale \(\mathbb{E}[T]/\operatorname{tr}(\mathrm{RIV})\) with the analytic denominator, the production results are: Bernoulli FITTED \(N\)-scan asymptote \(\hat a=+0.512\pm0.037\) (weighted fit of \(a+C/N\) over \(N=200/500/1000/2000\)) against the derived \(+0.500\), the \(R{=}200{,}000\) point at \(N{=}200\) being \(+0.5176\pm0.0116\) — which identifies the earlier \(+0.846\) as Monte-Carlo noise at \(R{=}1000\) on top of a real \(\approx+3.5/n\) finite-sample term of the same order as the exact-enumeration constants; the blind Poisson arm \(+0.113\pm0.027\) at \(N{=}200\) and \(+0.284\pm0.112\) at \(N{=}1000\), both within two standard errors of the registered \(+0.104\); the paired draw discount \(-0.313\pm0.070\) at \(N{=}2000\) against the derived \(-0.2606\), with a real decaying finite-\(n\) excess at \(N{=}200\) (\(-0.383\pm0.023\)); and the Student-\(t\) exemplar FITTED \(+0.025\pm0.049\) against \(0\) and proper draw \(-0.483\pm0.049\) against \(-0.500\) (\(R{=}200{,}000\)). Two registered hygiene checks failed marginally, with recorded causes that touch no claim: a continuity check landed at \(2.03\) against a \(2.00\) combined-standard-error threshold because the earlier \(R{=}1000\) cache value was itself the outlier, and a \(3\%\) sample-RIV tolerance was registered against the \(M{=}\infty\) analytic value, omitting Rubin’s \((1+1/M)\) factor that the sample estimator deliberately carries.
Source. Scan and scoring scripts verification/spikes/phase3-family-scan.R, verification/spikes/phase3-family-scan-score.R; caches verification/cache/phase3-family-scan-{bern_fitted,bern_draw,pois_fitted,t_exemplar}-prod.rds (registered criteria P1–P6 in todo/045 §3.2–§3.3, scored in §3.4). [mc] the quadrature constants rest on three independent numerical routes: Gauss–Hermite (verification/cas-wolfram/verify_bernoulli_fitted_gap.py), the blind re-derivation’s own quadrature (\(-0.260\)), and a composite-Simpson third route reproducing \(K=-1.042393\) and the Poisson parts \(+0.687997/-1.084067\) to six decimals (verification/cas-wolfram/verify_045b_escalations.py). Cross-model records: verification/cache/layer3-gpt-5.5-pro-bern-glm-mode{A,B}.md; repo-grounded compliance validation of the staged texts: tkal-reviews/gemini-045c-compliance-20260704/REPORT.md.
Fences. Four limits bound the result. (i) Canonical link only: the derivation uses \(\eta=x^\top\beta\) linear in \(X\) and \(\Omega\) sharing the weight \(A''\) at every occurrence; a non-canonical link breaks both, and the estimated-shape Gamma arm run alongside the scans (log link — not the Gamma family’s canonical link) therefore remains a directional check outside this lemma. (ii) Response-only missingness; general patterns are not derived. (iii) MCAR: the population-design substitution is where MCAR enters, and the MAR design-imbalance term of A4 is a separate \(O(1)\) not instantiated for this family. (iv) Asymptotic: the \(O(n^{-1})\) remainder’s constant is measured, not derived — exact enumeration gives \(\approx0.6/n\) (intercept-only) and \(\approx1.6/n\) (two-point) on the ratio scale, and the Gaussian-covariate scan design’s \(\approx+3.5/n\) comes from the \(N\)-scan itself. The draw corollary additionally assumes a second-moment-proper, regular draw in the sense of Section 2.5 — the \((C)\) discount is this family class’s analog of the posterior-draw correction stated there, and engines that are not second-moment proper are not covered by it.
The theorem prices the bias of one model’s \(\bar Q_\infty\). A likelihood-ratio comparison consumes two such quantities at once — the same imputed data, two nested fits — and everything in this section follows from asking what survives the subtraction. Three things are derived here: the closed-form leading-order differential bias of the MI deviance numerator at the null (it is not the difference of the two models’ \(\operatorname{tr}(\mathrm{RIV})\)s); the collapse of the per-dataset noise from \(O_p(\sqrt n)\) to \(O_p(1)\) under pairing; and the exact cancellation of the design-imbalance term \((A)+(C)\) at the null.
Let the null model be a smooth submodel \(\theta = g(\gamma)\), \(\gamma\in\mathbb{R}^{k_0}\), with Jacobian \(G = \partial g/\partial\gamma^\top\) (full column rank, evaluated at \(\gamma_0\)) and \(q_{\mathrm{d}} = k - k_0\) tested constraints; under the null the truth is \(\theta_0 = g(\gamma_0)\). Both models are fit to the same imputations: with \(\bar Q_\infty\) the averaged log-likelihood of Section 2.5 (imputation parameter the full model’s \(\hat\theta_{\text{obs}}\)), the MI likelihood-ratio numerator is \[ \hat d_L = 2\big[\bar Q_\infty(\hat\psi^*) - \bar Q_\infty(\hat\psi_0^*)\big], \qquad \hat\psi^* = \arg\max_\theta \bar Q_\infty(\theta), \quad \hat\psi_0^* = \arg\max_{\theta = g(\gamma)} \bar Q_\infty(\theta), \] the infinite-\(M\) limit of the averaged-likelihood maximizers of Chan and Meng — note \(\hat\psi_0^*\) maximizes the shared \(\bar Q_\infty\) over the null manifold; it is not the null model’s own observed-data MLE. The complete-data counterpart is \(\hat d_{\text{com}} = 2[\ell_{\text{com}}(\hat\theta_{\text{com}}) - \ell_{\text{com}}(g(\hat\gamma_{\text{com}}))]\), with \(\mathbb{E}[\hat d_{\text{com}}] = q_{\mathrm{d}} + O(n^{-1})\) at the null by Wilks’ theorem. The object of interest is the differential \(\hat d_L - \hat d_{\text{com}}\): the bias the missing data inject into the deviance comparison.
Source. Averaged-likelihood maximizers \(\hat\psi^*, \hat\psi_0^*\) and \(\hat d_L = 2\{\bar L(\hat\psi^*) - \bar L(\hat\psi_0^*)\}\): Chan and Meng (2022, :191–97); the original average-the-LRT combining rule is Meng and Rubin (1992, (2.9), :168). Wilks’ \(\chi^2_{q_{\mathrm{d}}}\) limit: Wilks (1938, Theorem, :99).
Proposition L1 (differential bias). Under R1–R7 and the construction above, at the null, \[ \mathbb{E}\big[\hat d_L - \hat d_{\text{com}}\big] = \operatorname{tr}(\mathrm{RIV}_\perp) + O(n^{-1}), \qquad \operatorname{tr}(\mathrm{RIV}_\perp) := \operatorname{tr}\big(I_{\text{obs}}^{-1} I_{\text{mis}\mid\text{obs}}\big) - \operatorname{tr}\big[(G^\top I_{\text{com}} G)^{-1}\, G^\top I_{\text{com}} I_{\text{obs}}^{-1} I_{\text{mis}\mid\text{obs}}\, G\big]. \tag{13}\] The null-model term is the missing information projected onto the null tangent space in the \(I_{\text{com}}\) metric — not the null model’s own self-contained \(\operatorname{tr}(\mathrm{RIV}_0) = \operatorname{tr}[(G^\top I_{\text{obs}} G)^{-1} G^\top I_{\text{mis}\mid\text{obs}} G]\).
Why \(I_{\text{com}}\) and not \(I_{\text{obs}}\). Write \(Z = \nabla\bar Q_\infty(\theta_0)\) and \(J = -\nabla^2\bar Q_\infty(\theta_0)\). Because the conditional density integrates to one, the evaluation-slot gradient of the imputed part vanishes identically at \((\theta_0,\theta_0)\), so \(Z = S + I_{\text{mis}\mid\text{obs}}\,\delta_{\text{obs}} + O_p(1)\) with \(\delta_{\text{obs}} = \hat\theta_{\text{obs}}-\theta_0\) — and hence \(Z = I_{\text{com}} I_{\text{obs}}^{-1} S + O_p(1)\), \(\operatorname{Var}(Z) = I_{\text{com}} I_{\text{obs}}^{-1} I_{\text{com}} + O(\sqrt n)\). The curvature is \(J \to_p I_{\text{obs}} + I_{\text{mis}\mid\text{obs}} = I_{\text{com}}\) — the EM identity that the \(Q\)-function’s curvature is the complete-data information. The two are consistent only together: the unconstrained maximizer must reproduce \(\hat\theta_{\text{obs}}\) (EM self-consistency, the fixed point of Section 2.5), and indeed \(J^{-1}Z = I_{\text{com}}^{-1} I_{\text{com}} I_{\text{obs}}^{-1} S = \delta_{\text{obs}}\) only for \(J = I_{\text{com}}\). The constrained maximizer therefore projects \(Z\) onto \(\operatorname{col}(G)\) in the \(I_{\text{com}}\) metric, the deviance is the standard difference of quadratic forms \(\hat d_L = Z^\top[I_{\text{com}}^{-1} - G(G^\top I_{\text{com}}G)^{-1}G^\top]Z + O_p(n^{-1/2})\), and taking expectations against \(\operatorname{Var}(Z)\), subtracting Wilks’ \(\mathbb{E}[\hat d_{\text{com}}] = q_{\mathrm{d}} + O(n^{-1})\), gives Equation 13. Substituting \(I_{\text{obs}}\) for the curvature — conflating \(\hat\psi_0^*\) with the null model’s own observed-data MLE — is the natural error here; it reproduces exactly the naive \(\operatorname{tr}(\mathrm{RIV}) - \operatorname{tr}(\mathrm{RIV}_0)\) of Proposition L2.
Source. Derived this session, three independent routes (difference-of-levels, direct quadratic expansion, realized pairing) with independent referee recomputation (verification/cache/w2diff-stage1-routes.json); the realized-pairing route made precisely the \(I_{\text{obs}}\)-metric error and its referee caught and corrected it. Route-form equivalences and the projection algebra verified symbolically in verification/cas-wolfram/verify_w2_differential.py (checks A1–A4). MLE expansion and quadratic-form geometry: Vaart (1998) [Thm 5.39, :7388; Thm 16.7, :7457].
Numerical confirmation. On a four-variate design built to separate the two candidate formulas — \(H_0\!: \mu_1 = 0\) with \(X_1\) missing for most units under MAR, where Equation 13 gives \(2.64\) while the naive difference gives \(8.47\) — the Monte-Carlo mean of \(\hat d_L - \hat d_{\text{com}}\) lands on Equation 13 to within sampling error and \(52\) standard errors away from the naive value (\(2.639 \pm 0.113\) at \(N{=}500\), \(2.701 \pm 0.112\) at \(N{=}1000\), \(R = 2000\); verification/w2-paired-differential-mu1-cell.R). On the \(\sigma_{12}\)-test design of the pre-registered likelihood-ratio study the two formulas differ by only \(\approx 5\%\) (\(1.31\) vs \(1.37\)) and the Monte Carlo is consistent with both (verification/w2-paired-differential-spotcheck.R).
Proposition L2 (overstatement). Always \(\operatorname{tr}(\mathrm{RIV}_\perp) \le \operatorname{tr}(\mathrm{RIV}) - \operatorname{tr}(\mathrm{RIV}_0)\), with equality if and only if \(\operatorname{col}(I_{\text{obs}}^{1/2}G)\) is an invariant subspace of \(I_{\text{obs}}^{-1/2} I_{\text{mis}\mid\text{obs}} I_{\text{obs}}^{-1/2}\) — in particular under equal fractions of missing information, where both sides equal \(r\,q_{\mathrm{d}}\) for the common odds \(r\). The proof is one matrix inequality: the Gram matrix of \((I_{\text{obs}}^{1/2}G,\; I_{\text{obs}}^{-1/2}I_{\text{com}}G)\) is positive semidefinite, so its Schur complement \(G^\top I_{\text{com}} I_{\text{obs}}^{-1} I_{\text{com}} G \succeq (G^\top I_{\text{com}}G)(G^\top I_{\text{obs}}G)^{-1}(G^\top I_{\text{com}}G)\), and tracing against \((G^\top I_{\text{com}}G)^{-1}\) gives the claim.
Consequently, correcting an MI deviance comparison by the difference of the two models’ self-contained \(\operatorname{tr}(\mathrm{RIV})\)s over-corrects at the null whenever the tested directions mix unequally-missing information — mildly in the W2 design (\(\approx 5\%\)), by a factor above three in the \(\mu_1\) design above. The preregistered W2 study’s correction differential used the naive difference; its conclusions are unaffected there (the design sits near the invariant case), but Equation 13 is the correct general form.
The overstatement has an exact closed form, which makes the equality condition transparent. Write \(H = I_{\text{obs}}^{-1/2} I_{\text{mis}\mid\text{obs}} I_{\text{obs}}^{-1/2}\) for the standardized missing information, and partition it by the retained model directions (\(\operatorname{col}(I_{\text{obs}}^{1/2}G)\), block \(H_{11}\)) versus the tested directions (block \(H_{22}\), coupling \(H_{12}\)). Then \[\big[\operatorname{tr}(\mathrm{RIV}) - \operatorname{tr}(\mathrm{RIV}_0)\big] - \operatorname{tr}(\mathrm{RIV}_\perp) = \operatorname{tr}\!\big[(I + H_{11})^{-1} H_{12} H_{12}^\top\big] \;\ge\; 0,\] positive semidefinite and zero exactly when \(H_{12} = 0\) — when the tested and retained directions carry independent missing information. The overstatement is therefore precisely the missing-information coupling between the two subspaces, screened by the retained block’s own information \(I + H_{11}\); it vanishes in the equal-fractions and block-diagonal cases because there the coupling does.
Source. Inequality + equality condition: verification/cas-wolfram/verify_w2_differential.py (B1–B3: symbolic Gram/Schur check, strictness over 400 random instances, exact equality under constructed invariance; B4: the closed-form gap \(\operatorname{tr}[(I+H_{11})^{-1}H_{12}H_{12}^\top]\) reproduces the direct difference to \(6\times10^{-14}\)); independently derived in operator-variance form by the quadratic route (w2diff-stage1-routes.json) and as a red-team byproduct of the cross-model audit (Section 5.5 source note).
Proposition L3 (paired noise). Define the per-dataset paired differential \(D = (\hat d_L - \hat d_{\text{com}})/2\). At the null and under local alternatives \(\theta_n = g(\gamma_0) + h/\sqrt n\), \[ \operatorname{sd}(D) = O(1), \] in contrast to the single-model statistic \(T\) of Equation 4, whose standard deviation grows as \(c\sqrt n\) (Section 4.4). At a fixed alternative — the two models’ pseudo-true values separated at \(O(1)\), as in non-nested comparison (Vuong 1989) — the cancellation fails and \(\operatorname{sd}(D)\) reverts to \(O(\sqrt n)\).
The mechanism is worth stating precisely, because it is not a projection of the linear noise onto the tested subspace. The \(O_p(\sqrt n)\) noise of each model’s level lives in a fit-independent realized constant: the imputed conditional term \(F(\theta_0, \hat\theta_{\text{obs}})\)-type pieces (the entropy gradient contracted with the shared \(\delta_{\text{obs}}\)) plus the realized \(\ell_{\text{com}} - \ell_{\text{obs}}\) at \(\theta_0\). (Argument-order note: in this section \(F\) takes the evaluation parameter first and the sampling/imputation parameter second, \(F(\psi, \phi)\) — the reverse positional order of Section 4.5’s \(F(\phi, \psi)\); the roles are identical, with \(\phi\) the imputer’s parameter in both.) Both fits maximize the same realized \(\bar Q_\infty\) built from the same imputation parameter, so these pieces are literally the same realization in both levels and cancel exactly, dataset by dataset — not merely in expectation. What survives is a difference of \(O_p(1)\) quadratic forms. Empirically, on the W2 design at the null, \(\operatorname{sd}(T) = 18.4 \to 26.6 \to 36.5\) across \(N = 500 \to 1000 \to 2000\) (the \(\sqrt n\) growth) while \(\operatorname{sd}(D) = 1.58 \to 1.47 \to 1.57\) (flat); at the fixed alternative \(\sigma_{12} = 0.6\), \(\operatorname{sd}(D) = 12.4 \to 16.8 \to 24.9\) — the \(\sqrt n\) tail returns (verification/w2-paired-differential-spotcheck.R). This is the precise content of the “structurally better-conditioned” clause of Section 4.1, and the design basis for estimating the comparison-relevant correction by paired contrasts rather than differences of separately-estimated levels.
Proposition L4 (\((A)+(C)\) differential). The \((A)+(C)\) contributions of the two levels are built from the entropy gradient \(\alpha\), the curvature \(\mathcal{H}_{\phi\phi}\), and the shared imputation error \(\delta_{\text{obs}}\) — all properties of the imputation, not of which model is fit. At the null they are identical realizations in both levels and cancel exactly to leading order; the first surviving model-specific contribution is \(O(n^{-1})\). Under local alternatives the differential is \(O(n^{-1/2}) = o(1)\): the trace differential Equation 13 dominates whatever the \((A)+(C)\) term does locally. Only at a fixed alternative — pseudo-true values \(O(1)\) apart, the non-nested regime of Vuong (1989) — do the two models’ \((A)+(C)\) terms decouple into a genuine \(O(1)\) differential, alongside the return of the \(O(\sqrt n)\) noise of Proposition L3.
This sharpens the standalone-value discussion of Section 4.4: the \((A)+(C)\) level is an \(O(1)\) whose magnitude is hard to pin (heavy-tailed, \(\operatorname{sd}\propto\sqrt n\)), but the LRT consumes only its differential, which is exactly zero at the null, asymptotically negligible under contiguous alternatives, and re-emerges only for non-nested comparison.
A direct measurement sharpens this, with one quantitative refinement. Fitting compound-symmetry and AR(1) covariance candidates to the same congenial imputations under a strong \(X_1\)-heavy MAR design — a genuine non-nested comparison (Vuong 1989) with the two pseudo-true correlations well separated (\(0.35\) vs \(0.50\); \(0.40\) vs \(0.47\)) — and isolating the design-imbalance contribution by differencing against an MCAR twin (\(20{,}000\) replicates per cell, \(N \in \{500, 1000, 2000\}\)), the run reproduces the mechanism of Proposition L4. The single-model \((A)+(C)\) level is a large \(O(1)\) for each curved candidate — \(\approx +2.5\) deviance units, tens of standard errors from zero — that collapses to zero under MCAR: the two candidates’ levels decouple (\(+2.33\) vs \(+2.85\) in the AR(1)-truth cell), no longer the identical realizations that cancel at the null, and the paired noise returns to \(O(\sqrt n)\) (the same run shows this for any pseudo-true-separated pair — a nested-but-false restriction triggers it just as a non-nested one does, so it is separation, not nesting, that matters). What is not large is the differential of the two levels: compound symmetry and AR(1) respond to the imbalance so similarly that their \((A)+(C)\) contributions nearly cancel, leaving an isolated differential consistent with zero (below \(\approx 0.5\) at \(N = 2000\); the pilot’s apparent \(-2\) was Monte-Carlo noise the full run dissolved). The un-absorbed \(O(1)\) that an information criterion would consume is therefore real per candidate but mild for similar candidates; its magnitude scales with how differently the competitors handle the missing data. A second, maximally dissimilar pair confirms this directly: a structure-blind diagonal model against AR(1) (same truth, same isolation, same \(R\)) carries an isolated design-imbalance differential that is a resolved \(O(1)\) — \(-1.2\) at \(-3.5\) standard errors, three-to-four times the similar pair’s and un-trending across \(N\), where the similar pair’s sat below resolution. The mechanism is legible: the correct, richer AR(1) is more entangled with the \(X_1\)-heavy missingness than the variance-only diagonal, so it carries the larger \((A)+(C)\) (\(+2.4\) vs \(+1.2\)). One honest caveat travels with the dissimilar case: a badly misspecified candidate also carries a mechanism-independent misspecification \(O(1)\) (visible as a non-vanishing MCAR residual), which the MCAR contrast separates from the design-imbalance term; the two can partially offset in the net bias an information criterion sees, so “the \((A)+(C)\) differential is a resolved \(O(1)\)” is a statement about the design-imbalance term itself, not a guarantee that it dominates the net ranking bias for every pair.
Source. Local-alternative machinery: Vaart (1998, secs. 16.4, :7553–:7570). Pseudo-true values and the nested/non-nested LR framework: Vuong (1989, :81, :40). Monte-Carlo check of the null cancellation: w2-paired-differential-spotcheck.R (mean of \(D\) at the null matches Equation 13\(/2\) with no \(O(1)\) offset, \(\pm 0.02\) pooled). Non-nested measurement (decoupled levels, \(\sqrt n\) noise return, small differential): pre-registered todo/029 (§6 grades M1–M4); script verification/w4-nonnested-vuong.R; caches verification/cache/w4-nonnested-summary.csv and the 12 per-(cell\(\times\)N) files (e.g. verification/cache/w4-nonnested-A-N500.rds); CAS verification/cas-wolfram/verify_w4_nonnested.py (general-Jacobian spectrum identity to \(2\times10^{-15}\) on the curved \(G_k\)). Dissimilar-pair anchor (Diagonal vs AR(1); resolved \(-1.2\) differential at \(-3.5\) se): todo/029 Amendment 2; script verification/w4-dissimilar-vuong.R; caches verification/cache/w4-dissimilar-summary.csv + verification/cache/w4-dissimilar-A-N500.rds.
At the null, Equation 13 is precisely the mean inflation that a correctly calibrated reference distribution already absorbs: the limiting law of \(\hat d_L\) is the weighted sum \(\sum_{j=1}^{q_{\mathrm{d}}} \lambda_j \chi^2_1\) whose weights \(\lambda_j = 1 + r_j\) exceed one by the odds-of-missing-information of the tested directions (Chan 2022, eq. (1.10)), and \(\sum_j \lambda_j = q_{\mathrm{d}} + \operatorname{tr}(\mathrm{RIV}_\perp)\). A test that refers the uncorrected numerator to such a reference is therefore approximately calibrated at the null — consistent with the preregistered W2 null rates (uncorrected against Chan’s simulated reference: \(0.042\); numerator additionally corrected: \(0.034\), slightly conservative, as double-counting predicts). The differential bias matters instead where no reference distribution absorbs it: information-criterion-style ranking, which compares penalized \(-2\bar Q_\infty\) values directly (Section 1.3); numerator corrections, which must use Equation 13 rather than the naive difference (Proposition L2); and non-nested comparison, where both the \((A)+(C)\) differential and the heavy-tailed noise return (Proposition L4). The first of these — information-criterion ranking — is taken up in Section 6, together with the question of how much more than the mean a correction can match.
Scope. The propositions are stated for the deterministic-FIML \(\bar Q_\infty\) (imputation at the point \(\hat\theta_{\text{obs}}\)). Under proper MI the imputation-side quantities acquire posterior smearing, but they remain common to both fits — the same posterior, the same draws — so the cancellation structure of Propositions L3–L4 and the leading-order form of Equation 13 carry over; the Monte-Carlo confirmations above are deterministic-FIML. A blind off-lineage re-derivation (a different model family, given only the setup and the tasks) independently reproduced every load-bearing conclusion — the \(I_{\text{com}}\) curvature and projection metric, Equation 13 term for term, the definite-sign overstatement of Proposition L2 with a third independent proof, and the noise orders and exact \((A)+(C)\) cancellation of Propositions L3–L4 — and a companion adversarial pass, instructed to break each claim (including the natural objection that the two fits’ \(O_p(n^{-1/2})\)-separated maximizers should leave \(O_p(\sqrt n)\) uncancelled noise), sustained all four, locating the objection’s error precisely: it attaches the imputation-slot gradient to an evaluation-slot displacement that the conditional-score identity annihilates (verification/cache/layer3-gpt-5.5-w2diff-modeA.md and verification/cache/layer3-gpt-5.5-w2diff-modeB.md; protocol and grading record in todo/024-w2-differential-analytic-first.md).
The likelihood-ratio application referred a bias-corrected statistic to a reference distribution. Model selection by information criterion has no reference distribution to hide behind: it compares penalized deviances as numbers, across candidates, and takes the smallest. Whatever bias each candidate’s \(-2\bar Q_\infty\) carries lands directly in the ranking. This section does two things: it states the corrected criterion and places it in the missing-data information-criterion literature (the penalty itself turns out to be known; the decomposition, the MAR term, and the proper-MI extension are what is new), and it then asks a sharper question the literature has not: since selection depends on the distribution of pairwise criterion differences, not just their means, how much of that distribution can a correction match? The answer — derived from the Section 5 null law, pre-registered, then measured — splits cleanly along the null/noncentral axis.
The MI information criterion built on the averaged deviance is \(\mathrm{AIC}_{\mathrm{MI}}(k) = -2\bar Q_\infty(\hat\psi^*_k) + 2q_k\), with \(q_k\) the candidate’s parameter count. By the theorem, each candidate’s \(\bar Q_\infty\) overstates its complete-data-target log-likelihood by \(\tfrac12\operatorname{tr}(\mathrm{RIV}_k)\) — so the deviance is too small by \(\operatorname{tr}(\mathrm{RIV}_k)\), and the bias is model-specific: \(\mathrm{RIV}_k\) is computed on the candidate’s own parameter space. The corrected criterion is \[ \mathrm{AIC}_{\mathrm{MI}}^{c}(k) = -2\bar Q_\infty(\hat\psi^*_k) + 2q_k + \operatorname{tr}(\mathrm{RIV}_k). \tag{14}\] As a penalty this is not new: it reproduces the missing-data information criterion \(\mathrm{AIC}_{x;y}\) of Shimodaira and Maeda (2018), whose missing-data surcharge \(\operatorname{tr}(I_{z\mid y} I_y^{-1}) = \operatorname{tr}(\mathrm{RIV})\) is exactly half that of the earlier complete-data-discrepancy criteria — Cavanaugh and Shumway’s \(\mathrm{AIC}_{cd}\) (Cavanaugh and Shumway 1998) and the related PDIO attach \(2\operatorname{tr}(\mathrm{RIV})\) to the EM deviance, which Shimodaira and Maeda (2018) showed doubles the correct value. What the present derivation adds is the anatomy and the scope of that penalty: the decomposition into the imputation-bias term \(+\operatorname{tr}(\mathrm{RIV})\) and the estimation-mismatch term \(-\tfrac12\operatorname{tr}(\mathrm{RIV})\); the MAR design-imbalance term \((A)+(C)\), which lies beyond the fixed-pattern setting of the earlier criteria; and the extension from the deterministic EM pivot to proper multiple imputation (Section 1.4; the alternative EM-likelihood and reweighting routes of Claeskens and Consentino (2008) and Hens, Aerts, and Molenberghs (2006) are positioned there as well).
Source. \(\mathrm{AIC}_{x;y} = -2\ell_y(\hat\theta_y) + d +
\operatorname{tr}\!\big(I_x(\hat\theta_y)\, I_y(\hat\theta_y)^{-1}\big)\) and the halved surcharge: Shimodaira and Maeda (2018) [eq. (8), :80; eq. (9), :86; the trace decomposition at :59]. \(\mathrm{AIC}_{cd}\) with the doubled trace: Cavanaugh and Shumway (1998) [eq. (3.20), :302; the \(2d + 2\operatorname{tr}\) split, eq. (4.4), :353], reproduced as eq. (4) of Shimodaira and Maeda (2018, :68).
The ranking consequence is immediate. For two candidates the uncorrected criterion difference carries a differential bias whose leading term is \(-[\operatorname{tr}(\mathrm{RIV}_j) - \operatorname{tr}(\mathrm{RIV}_k)]\): the candidate with more missing information about its own parameters looks artificially better, by exactly the amount of extra missing information. Two caveats from Section 5 frame what the correction can and cannot promise. First, candidate sets are generally not nested chains, and for non-nested pairs the \((A)+(C)\) differential is a genuine \(O(1)\) that Equation 14 does not remove and no reference distribution absorbs (Proposition L4; the Vuong regime (Vuong 1989)) — though the direct measurement of Section 5.5 finds this differential small when the candidates are similar, the two large per-model \((A)+(C)\) levels nearly cancelling, so it bites in proportion to candidate dissimilarity. Second, the per-model \(\operatorname{tr}(\mathrm{RIV}_k)\) corrections are exactly the “self-contained levels” whose difference Proposition L2 shows overstates the projected complement trace whenever the tested directions mix unequally-missing information — though, remarkably, the model family studied below sits at the exact equality case (a structural fact derived there, not luck).
The pre-registered model-selection study draws four-variate normal data whose only nonzero correlations are \(\sigma_{13} = \sigma_{24} = 0.4\), imposes 40% MAR missingness, and ranks four zero-pattern covariance candidates — diagonal (MA), wrong off-diagonals (MB), the truth (MC), and saturated (MD) — by MI-AIC under congenial imputation. Across the congenial cells of the pre-registered sweep, uncorrected MI-AIC selects the true model in \(\approx 0.65\)–\(0.67\) of replicates against an oracle (complete-data AIC) rate of \(\approx 0.90\) — and 100% of its misclassifications land on the saturated, largest-RIV candidate, the directional signature of the model-specific bias. The corrected criterion Equation 14 recovers selection to \(\approx 0.81\) — substantially, but visibly short of the oracle. That gap is not sloppiness in the trace; it is the subject of the next subsection.
Source. Committed sweep, 60 cells: verification/cache/run_all-prod/phase2-w3/summary.csv (e.g. the deterministic-FIML \(N{=}500\) cell: oracle \(0.898\), uncorrected \(0.674\), corrected \(0.816\); W3C_high_RIV_frac = 1.0 in every congenial cell). Engine robustness: the Amelia/EMB congenial cells agree with deterministic FIML within \(0.011\) on the recovery gap.
Selection depends on more than a mean. Even with every candidate’s criterion correctly centered, the spread of the MI deviance differences exceeds its complete-data counterpart — at the null the limiting law is \(\sum_j \lambda_j \chi^2_1\) with every weight \(\lambda_j = 1 +
r_{\perp,j} \ge 1\) (Section 5.6), so the variance is inflated by \(\sum\lambda_j^2/q_{\mathrm{d}}\) even after the mean is fixed — and inflated spread means flipped rankings. This motivates a graded contrast, pre-registered in full before any code was written (todo/027-ladder-preregistration.md): match the first moment, the first two moments, or the entire null distribution, and measure what each step buys.
The working statistic. Anchor every candidate at the saturated model, fit to the same imputations: \(T_k = 2[\bar Q_\infty(\hat\psi^*_{\text{sat}}) - \bar Q_\infty(\hat\psi^*_k)] \ge 0\). Anchoring costs nothing — AIC ranking is invariant to the common shift — and buys three things. The anchor is the congenial imputation model, so the heavy-tailed realized \((A)+(C)\) component and the \(O_p(\sqrt n)\) noise common to both levels cancel dataset by dataset (the Proposition L3 mechanism, applied for free to every candidate); every \(T_k\) is a proper likelihood-ratio statistic with a constructible per-model null reference (the distribution of \(T_k\) when candidate \(k\) is true, available analytically from the Section 5 law without knowing the global truth); and each candidate here is a block-diagonal zero pattern, so the constrained maximizer \(\hat\psi^*_k\) is closed-form — the block restriction of the E-step moments — giving \(T_k = N[\log|\hat\Sigma_k| - \log|S^*|]\) exactly.
The three corrections. With \(\lambda_{k,1},\dots,\lambda_{k,q_{d,k}}\) the analytic null-law weights of candidate \(k\) — the nonzero eigenvalues of \([I_{\text{com}}^{-1} - G_k(G_k^\top I_{\text{com}} G_k)^{-1}G_k^\top]\, I_{\text{com}} I_{\text{obs}}^{-1} I_{\text{com}}\), with \(\sum_j \lambda_{k,j} = q_{d,k} + \operatorname{tr}(\mathrm{RIV}_{\perp,k})\) and each \(\lambda_{k,j} \ge 1\) —
Source. Scaled difference statistic \(\bar T_d = T_d/\hat c_d\) and the Satorra–Bentler (2001) combination: Satorra and Bentler (2010, :61, :79). Mean-and-variance adjustment via \(\bar T = (d/\operatorname{tr}(M))T\) with \(d' = (\operatorname{tr} M)^2/\operatorname{tr}(M^2)\): Asparouhov and Muthén (2006, :51, :67). Equipercentile equating \(e_Y(x) = G^{-1}[F(x)]\): Kolen and Brennan (2014, eq. (2.11), :1262). Imhof’s inversion integral: Imhof (1961, eq. (3.2), :237); algorithmic form: Davies (1980, :47). Monte-Carlo null reference: Chan (2022, eq. (5.1), :404). The moment identities, the uniqueness and form-equivalence of \((a_k, b_k)\), the spectrum-trace identity \(\sum_j\lambda_j = q_d + \operatorname{tr}(\mathrm{RIV}_\perp)\) with \(\lambda_j \ge 1\), and the probability-integral-transform exactness of the equated statistic are verified symbolically and numerically in verification/cas-wolfram/verify_ic_ladder.py (L-A1–L-A4, all pass).
What was predicted. The pre-registration (written before the script existed) states the split: each step should close more of the null-side calibration gap, with equating exact there up to asymptotic approximation; on the noncentral side — a misspecified candidate’s \(T_k\) against its own null map — the first moment is preserved only by the mean correction, while the two-moment and equating maps actively shrink the evidence (the affine map by \(a_k < 1\), the equating map by \(\approx 1/\lambda_{\max}\) in the far tail), because no transform built from observed data can restore destroyed Fisher information. The residual against the oracle should therefore be removable where selection errors are overfit flips, and irreducible where they are underfit retention.
Three cells, \(N = 500\), \(R = 2000\) each, 40% non-monotone MAR: the pre-registered design (\(\rho = 0.40\), truth MC), a weak-signal stress cell (\(\rho = 0.10\), near-tied candidates), and a junk-heavy cell (\(\rho = 0\), truth diagonal, so every anchored statistic is null-type). True-model selection rates (Monte-Carlo se \(\approx 0.010\)):
| Cell | oracle | uncorrected | mean | two-moment | equating |
|---|---|---|---|---|---|
| \(\rho = 0.40\) | 0.899 | 0.678 | 0.814 | 0.904 | 0.903 |
| \(\rho = 0.10\) | 0.820 | 0.582 | 0.650 | 0.648 | 0.650 |
| \(\rho = 0\) (junk) | 0.727 | 0.422 | 0.585 | 0.702 | 0.696 |
Every pre-registered criterion passed; the detailed grading is recorded in the pre-registration file. Four findings carry the story.
The null side calibrates completely. On the true model’s anchored statistic, the Kolmogorov–Smirnov distance to the paired complete-data statistic falls \(0.281 \to 0.228 \to 0.026 \to 0.019\) across the four arms — the equated statistic is statistically indistinguishable from its complete-data counterpart at this replication size — and the variance ratio falls from \(3.18\) (uncorrected and mean-corrected) to \(0.97\). The mean correction’s paired gap, \(0.076 \pm 0.100\), is itself a direct Monte-Carlo confirmation of Equation 13 at \(q_d = 4\). In the junk cell, where all three statistics are null-type, equating calibrates all of them (\(D_{KS} \le 0.035\)).
Where overfit flips drive the errors, distribution matching reaches the oracle. In the \(\rho = 0.40\) design the corrected criterion’s residual (\(0.814\) vs \(0.899\)) is variance, not mean: the two-moment and equating arms select the true model at \(0.904\) and \(0.903\) — statistically indistinguishable from complete-data AIC. The junk cell shows the same mechanism at larger scale: the full contrast closes 90% of the uncorrected-to-oracle gap.
The noncentral side shrinks as the law predicts. The underfit candidate’s mean statistic, relative to its complete-data counterpart: \(1.002\) under the mean correction (evidence preserved), \(0.586\) under two moments (predicted \(a_k = 0.570\)), \(0.362\) under equating (predicted far-tail factor \(1/\lambda_{\max} = 0.398\)). This is the toy made exact: the null map cannot see information loss, so it reads shrunken noncentrality as genuinely smaller misfit.
Where selection is genuinely hard, nothing closes the floor. In the weak-signal cell the oracle itself drops to \(0.820\) and no correction passes \(0.650\): the \(0.17\) residual is the information-loss floor, not fixable calibration error. The upper arms buy nothing there — their null-side calibration win is fully offset by underfit retention (they keep too-small models at rate \(0.270\) vs the oracle’s \(0.091\)). The light indirect-pair measurement in the same cell quantifies the coordinatewise limitation: the equated score difference for the near-tied non-nested pair has mean \(6.10\) and sd \(5.99\) against the complete-data \(\Delta\)AIC’s \(10.00\) and \(7.10\) — per-model marginal transforms cannot calibrate a difference distribution (the dependence between scores is invariant under marginal maps; here MI’s shared imputations mildly raise the score correlation, \(0.279\) vs \(0.226\)).
Source. Pre-registration with pass criteria frozen before code: todo/027-ladder-preregistration.md (§6 records the grading). Script: verification/ic-ladder.R; outputs verification/cache/ic-ladder-summary.csv, verification/cache/ic-ladder-MAIN.rds, verification/cache/ic-ladder-WEAK.rds, verification/cache/ic-ladder-JUNK.rds; log verification/cache/ic-ladder-runlog.txt. Internal asserts passed on every replicate: spectrum-trace identity to \(3\times10^{-15}\), all null-law weights \(\ge 1\), zero Imhof failures.
A structural remark closes the loop with Proposition L2. In this candidate family the naive trace difference and the exact complement trace agree to machine precision (\(3.6\times10^{-15}\)) — not luck: the candidates are block-diagonal covariance patterns, and at a block-diagonal evaluation point both \(I_{\text{com}}\) and the expected \(I_{\text{obs}}\) decompose over within-block versus cross-block coordinates, which is precisely the invariance condition under which Proposition L2 holds with equality. The overstatement is therefore structurally invisible here and material only for constraints that do not block-decompose — the \(\mu_1\) design of Section 5.2, where it is a factor above three. The earlier sweep’s scaled-shift arm (which applied its transform to the uncorrected statistic with naive-difference inputs) scores \(0.824\) in the \(\rho = 0.40\) cell against the derived two-moment arm’s \(0.904\) — the repaired construction is worth eight points of selection accuracy where calibration binds.
Honest fences. Stated once, bluntly: all-moment matching holds per pair, under that pair’s own null — nowhere else; noncentral matching is impossible in principle, and the two-moment and equating maps actively shrink evidence against misspecified candidates (their conservatism is the price of null-side calibration); per-model marginal equating cannot calibrate difference distributions in general; the affine and equating maps do not telescope, so the contrast is presented on common-reference (saturated-anchored) scores, with the indirect-pair distortion measured rather than assumed away; and none of this restores complete-data selection where the data no longer carry the information — the weak-signal floor is the honest summary. Selection-aware refinements beyond the per-model null maps — de-shrinkage of estimated noncentrality, joint (copula-aware) calibration of score vectors — are developed separately. These fences differ in standing, and it is worth being explicit about which is which: the first two are theorems with direct empirical confirmation — the noncentral shrinkage is the \(a_k\) and \(1/\lambda_{\max}\) law landed above, and difference-distribution non-calibration is the invariance of the score copula under per-model marginal transforms — whereas the last is an empirical floor, the measured \(0.17\) weak-signal residual, not a proof that no estimator built from the same observed data could do better.
Scope. As in Section 5, everything here is stated and measured for the deterministic-FIML \(\bar Q_\infty\); the anchoring cancellations are properties of sharing the imputation model, so the construction carries to proper MI unchanged, and the pre-registered sweep’s proper-MI (EMB) congenial cells track the FIML cells within \(0.011\) on the headline rates.
Consolidated provenance audit, grouped by reference. Each entry reproduces — verbatim from the literature/<citekey>.md Mathpix conversions — the specific passages the derivation invokes, annotated with the source’s own equation/theorem number and the literature/<citekey>.md line. The margin provenance blocks in the preceding sections keep a short snippet as a re-pagination-stable anchor and link here (@sec-src-<key>) for the full text.
Fair use and post-review deletion. These excerpts are reproduced for the transparency audit that JAIGP norms call for, so a reviewer can check each cited claim against its source without leaving the document. To keep each quotation interpretable we include a little surrounding context — the lead-in, the cited statement, and the following clause. These verbatim passages are temporary: they will be deleted once the blinded peer-review process concludes, leaving only the equation/theorem numbers and the literature/<citekey>.md line references. Should the venue or a reviewer object to verbatim reproduction at any earlier point, trim each entry to its number and line reference and drop the quoted text.
(Schafer 1997). Invoked by D1, D4, D6, D7, R1, L3, Term A step A1.
Complete-data likelihood factorization, eq (3.2) [:941–:947].
Viewing each term in (3.1) as a function of \(\theta\), it follows that \(l(\theta \mid Y)=l(\theta \mid Y_{obs})+\log P(Y_{mis} \mid Y_{obs}, \theta)+c\), where \(l(\theta/Y)=\log P(Y \mid \theta)\) denotes the complete-data log likelihood, \(l(\theta/Y_{obs})\) … the observed-data loglikelihood, and \(c\) an arbitrary constant. The term \(P(Y_{mis} \mid Y_{obs}, \theta)\), which we shall call the predictive distribution of the missing data given \(\theta\), plays a central role in EM because it captures the interdependence between \(Y_{mis}\) and \(\theta\).
EM decomposition of the Q-function [:949–:964].
Because \(Y_{mis}\) is unknown we cannot calculate the second term on the right-hand side of (3.2), so instead we take the average of (3.2) over the predictive distribution \(P(Y_{mis} \mid Y_{obs}, \theta^{(t)})\) … This averaging yields \(Q(\theta \mid \theta^{(t)})=l(\theta \mid Y_{obs})+H(\theta \mid \theta^{(t)})+c\), where \(Q(\theta\mid\theta^{(t)})=\int l(\theta\mid Y)\,P(Y_{mis}\mid Y_{obs},\theta^{(t)})\,dY_{mis}\) and \(H(\theta\mid\theta^{(t)})=\int \log P(Y_{mis}\mid Y_{obs},\theta)\,P(Y_{mis}\mid Y_{obs},\theta^{(t)})\,dY_{mis}\).
Missing-information principle and the information identity, eq (3.19) [:1480–:1504].
Taking the expectation of this over \(P(Y_{mis}\mid Y_{obs},\theta)\), we obtain a fundamental relationship: the complete information is equal to the observed information plus the missing information. This relationship, called the missing information principle by Orchard and Woodbury (1972) … Assuming sufficient regularity to interchange the order of differentiation and integration, we can write the complete information as \(I_c(\theta)=-\partial^2 Q(\theta\mid\theta)/\partial\theta^2\), the observed information as \(I_o(\theta)=-\partial^2 l(\theta\mid Y_{obs})/\partial\theta^2\), and the missing information as \(I_m(\theta)=-\partial^2 H(\theta\mid\theta)/\partial\theta^2\), so that \(I_c(\theta)=I_o(\theta)+I_m(\theta)\) (3.19). Note that each quantity in (3.19) is a function of \(Y_{obs}\) although this fact has been suppressed in the notation.
Observed-data MLE [:1504–:1510].
A natural large-sample estimate of the covariance matrix with incomplete data is \(I_o^{-1}(\hat\theta)\) … where \(\hat\theta\) is now the observed-data MLE, the maximizer of \(l(\theta\mid Y_{obs})\).
EM rate matrix (fraction of missing information), eq (3.21) [:1514–:1520] — cited for the D8 disambiguation; this is not RIV.
Dempster, Laird and Rubin (1977) established an important connection between these information quantities and \(D\) … the asymptotic rate matrix of EM … they showed that \(D=I_c^{-1}(\hat\theta)\, I_m(\hat\theta)\) (3.21) … for a scalar parameter (3.21) implies that \(D\) is the ratio of the missing information to the complete information. For brevity we will call this ratio the fraction of missing information.
(Magnus and Neudecker 2019). Invoked by D14, D15, and Section 3.2.
Half-vectorization [:1652] (here for \(3\times3\) symmetric \(A\)).
\(\operatorname{vech}(A)=\left(a_{11}, a_{21}, a_{31}, a_{22}, a_{32}, a_{33}\right)^{\prime}\)
Duplication matrix, eq (23) [:1658], and its inverse [:1664].
\(D_{n} \operatorname{vech}(A)=\operatorname{vec} A \quad\left(A=A^{\prime}\right)\)
\(\operatorname{vech}(A)=\left(D_{n}^{\prime} D_{n}\right)^{-1} D_{n}^{\prime} \operatorname{vec} A \quad\left(A=A^{\prime}\right)\)
Inverse of \(D_n'(A\otimes A)D_n\) [:1692].
\(\left(D_{n}^{\prime}(A \otimes A) D_{n}\right)^{-1}=\left(D_{n}^{\prime} D_{n}\right)^{-1} D_{n}^{\prime}\left(A^{-1} \otimes A^{-1}\right) D_{n}\left(D_{n}^{\prime} D_{n}\right)^{-1}\)
MVN information matrix for \((\mu,\operatorname{vech}\Omega)\) [:1804-:1811] (\(\Omega\!\leftrightarrow\!\Sigma\), \(m\!\leftrightarrow\!p\), \(n\!\leftrightarrow\!N\) in our D15).
\(\mathcal{F}=n\begin{pmatrix}\Omega^{-1} & 0 \\ 0 & \frac{1}{2} D_{m}^{\prime}\left(\Omega^{-1} \otimes \Omega^{-1}\right) D_{m}\end{pmatrix}\)
(Vaart 1998). Invoked by D5, D6, R1, R2, L1, L3, and Section 3.2.
Score identities [:2320].
“\(P_{\theta} \dot{\ell}_{\theta}=0\) (scores have mean zero), and \(P_{\theta} \ddot{\ell}_{\theta}=-I_{\theta}\) (the curvature of the likelihood is equal to minus the Fisher information).”
MLE asymptotic normality [:2325].
\(\sqrt{n}\left(\hat{\theta}_{n}-\theta\right) \rightsquigarrow N\left(0, I_{\theta}^{-1}\right)\)
MLE linear (influence-function) expansion, Theorem 5.39 [:4170, :7388-:7393]; the joint full/null-restricted form at :7392-:7393 underpins L3’s score covariance.
\(\sqrt{n}\left(\hat{\theta}_{n}-\theta\right)=\frac{1}{\sqrt{n}} \sum_{i=1}^{n} I_{\theta}^{-1} \dot{\ell}_{\theta}\left(X_{i}\right)+o_{P_{\theta}}(1)\)
Conditional expectation as \(L^2\) projection (Ch 11) [:5161, :5170].
(line 5161) “the conditional expectation is the unique measurable function \(E(X \mid Y)\) of \(Y\) that satisfies the orthogonality relation \(\mathrm{E}(X-\mathrm{E}(X \mid Y)) g(Y)=0\)”; (line 5170, Example 11.3) “\(\mathrm{EX}=E E(X \mid Y)\) … ‘the expectation of a conditional expectation is the expectation.’”
Likelihood-ratio asymptotics, Theorem 16.7 [:7457]; local alternatives and the noncentral \(\chi^2\) (§16.4) [:7553-:7570]. Invoked by Section 5.
(line 7553) “We assume that the parameter \(\vartheta\) is an inner point of the parameter set and denote the true parameter by \(\vartheta+h / \sqrt{n}\). Under the conditions of Theorem 16.7, the sequence of likelihood ratio statistics is asymptotically distributed as \(\Lambda=\|Z+I_{\vartheta}^{1 / 2} h-I_{\vartheta}^{1 / 2} H_{0}\|^{2}\)”; (line 7566) “This distribution is known as the noncentral chi-square distribution with noncentrality parameter \(\delta\).”
(Meng 1994). Invoked by D1, D13, R4, R5, R6, L4. Complete data \(Z_{c}=\{X, Y_{inc}, I\}\), observed data \(Z_{o}=\{X, Y_{obs}, I, R_{inc}\}\) [:95].
Definition 1 (congeniality), eqs (2.3.1)-(2.3.2) [:102-:110].
\(\left[\widehat{\theta}\left(Z_{o}\right), U\left(Z_{o}\right)\right] \simeq\left[E_{f}\left[\theta \mid Z_{o}\right], V_{f}\left[\theta \mid Z_{o}\right]\right]\) (2.3.1)
\(\left[\widehat{\theta}\left(Z_{c}\right), U\left(Z_{c}\right)\right] \simeq\left[E_{f}\left[\theta \mid Z_{c}\right], V_{f}\left[\theta \mid Z_{c}\right]\right]\) (2.3.2)
Definition 2 (congeniality to an imputation model), eq (2.3.3) [:120].
\(f\left(Y_{mis} \mid Z_{o}\right)=g\left(Y_{mis} \mid Z_{o}, A\right)\)
Definition 5 (self-efficiency) [:450].
“A statistical estimation procedure \(\widehat{\theta}(\cdot)\) … is said to be self-efficient … if there is no \(\lambda \in(-\infty, \infty)\) such that the mean-squared error of \(\lambda \hat{\theta}\left(W_{o}\right)+(1-\lambda) \hat{\theta}\left(W_{c}\right)\) is less than that of \(\widehat{\theta}\left(W_{c}\right)\).”
Lemma 3 and Main result (variance recovery / efficiency) [:486, :503-:511]; conditions (a) self-efficiency and (c) “second-moment proper” imputation [:504, :506].
“(i) The repeated-imputation estimator is consistent for \(\theta\), and is at least as efficient as the analyst’s incomplete-data estimator.”
(Rubin 1987). Invoked by D8, R5, R7, L4, and the D9 notation-clash callout.
Pooled estimate and the variance components [:2161, :2167, :2173, :2179].
\(\bar{Q}_{m}=\sum_{l=1}^{m} \hat{Q}_{* l} / m\) ; \(\bar{U}_{m}=\sum_{l=1}^{m} U_{* l} / m\) ; \(B_{m}=\sum_{l=1}^{m}\left(\hat{Q}_{* l}-\bar{Q}_{m}\right)^{t}\left(\hat{Q}_{* l}-\bar{Q}_{m}\right) /(m-1)\) ; \(T_{m}=\bar{U}_{m}+\left(1+m^{-1}\right) B_{m}\)
Relative increase in variance due to nonresponse [:2194] (named at :2191).
\(r_{m}=\left(1+m^{-1}\right) B_{m} / \bar{U}_{m}\)
Randomization-validity of the infinite-\(m\) inference, §4.2 — conditions (4.2.1)-(4.2.2) [:3489] and Result 4.1 [:3563]. This is the variance-recovery property L4 invokes: under proper MI the infinite-\(m\) variance estimator \(T_\infty\) is valid for the repeated-sampling variance \(T_0\) of the complete-data statistic.
\(\left(\bar{Q}_{\infty} \mid X, Y\right) \sim N\left(Q, T_{0}\right)\) (4.2.1); \(\left(T_{\infty} \mid X, Y\right) \sim\left(T_{0}, \ll T_{0}\right)\) (4.2.2), “where \(Q=Q(X, Y)\) and \(T_{0} \equiv T_{0}(X, Y)\) are fixed by the true values of \((X, Y)\) … When (4.2.1) and (4.2.2) hold, we will say that the repeated-imputation inference is randomization-valid.”
(Result 4.1) “If the Complete-Data Inference Is Randomization-Valid and the Multiple-Imputation Procedure Is Proper, Then the Infinite-\(m\) Repeated-Imputation Inference Is Randomization-Valid under the Posited Response Mechanism.”
(Rubin 1976). Invoked by D2, R3 (ignorability).
Ignorability for likelihood / Bayesian inference (abstract) [:39].
“When making direct-likelihood or Bayesian inferences about \(\theta\), it is appropriate to ignore the process that causes missing data if the missing data are missing at random and the parameter of the missing data process is ‘distinct’ from \(\theta\).”
(Little and Rubin 2019). Invoked by D1, D2, R3.
MAR definition [:438-:444].
“missingness depends on \(y_{i}\) only through the observed components \(y_{(0)i}\) … The missingness mechanism is then called missing at random (MAR).”
Ignorability — Definition 6.4 and Theorem 6.1A(a) [:3263, :3277].
(Def 6.4) “The missingness mechanism is called ignorable if inferences for \(\theta\) based on the ignorable likelihood … are the same as inferences … based on the full likelihood”; (Thm 6.1A(a), parameter distinctness) “the joint parameter space \((\theta, \psi)\) … is the product of the parameter space \(\Omega_{\theta}\) of \(\theta\) and the parameter space \(\Omega_{\psi}\) of \(\psi\).”
(Anderson 2003). Invoked by Section 3.2 (Gaussian conditional / Schur complement) and Term A’s normal Q-function.
Conditional covariance (Schur complement) [:1582].
\(\Sigma_{11 \cdot 2}=\Sigma_{11}-\Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21}\)
Conditional mean [:1598]; the \(q\)-variate normal conditional density \(f(x^{(1)} \mid x^{(2)})\) has this mean and covariance \(\Sigma_{11\cdot2}\) [:1595].
\(\mathscr{E}\left(X^{(1)} \mid x^{(2)}\right)=\mu^{(1)}+\Sigma_{12} \Sigma_{22}^{-1}\left(x^{(2)}-\mu^{(2)}\right)\)
(Petersen and Pedersen 2012). Invoked by Section 3.2 (trace and vec identities) and the Equation 3 trace reduction.
§1.1 Trace [:177-:187]. In our Mathpix conversion this section is a single unnumbered aligned display, so the Cookbook’s own inline equation numbers are not recoverable from the file; these trace identities are therefore cited by section + line.
\(\operatorname{Tr}(\mathbf{A})=\sum_{i} A_{ii}\) [
:179]; \(\operatorname{Tr}(\mathbf{A})=\operatorname{Tr}(\mathbf{A}^{T})\) [:181]; \(\operatorname{Tr}(\mathbf{A B})=\operatorname{Tr}(\mathbf{B A})\) [:182]; \(\operatorname{Tr}(\mathbf{A}+\mathbf{B})=\operatorname{Tr}(\mathbf{A})+\operatorname{Tr}(\mathbf{B})\) [:183]; \(\operatorname{Tr}(\mathbf{A B C})=\operatorname{Tr}(\mathbf{B C A})=\operatorname{Tr}(\mathbf{C A B})\) [:184]
§10.2.2 vec-operator, eqs (520)-(521) [:3345-:3346]. Equation numbers confirmed in-file: the vec rule is cited as “Eq 520” at :3692, and “Eq. 524” closes the same numbered block at :3353.
\(\operatorname{vec}(\mathbf{A X B})=\left(\mathbf{B}^{T} \otimes \mathbf{A}\right) \operatorname{vec}(\mathbf{X})\) (520)
\(\operatorname{Tr}\left(\mathbf{A}^{T} \mathbf{B}\right)=\operatorname{vec}(\mathbf{A})^{T} \operatorname{vec}(\mathbf{B})\) (521)
(Efron and Hinkley 1978). Secondary support for preferring the observed information: the observed information, not its expectation, tracks the finite-sample variance of the MLE. The MAR-specific reason — that the unconditional expected information itself departs from the naive form — is Kenward and Molenberghs (1998) (Section 7.11).
Observed vs expected Fisher information [:99].
“The quantity \(I(x)\) is aptly called the observed Fisher information by some writers, as distinguished from \(\mathscr{I}_{\theta}\), the expected Fisher information. This last name is useful even though \(E\left(\mathscr{J}_{\theta}\right) \neq \mathscr{I}_{\theta}\) in general.”
The MLE’s variance is the reciprocal of the observed information [:107–:113].
“Approximation (1•2), one over the expected Fisher information, would presumably never be applied in practice … the conditional variance \((1\cdot4)\) is a corresponding property of the observed likelihood: \((1\cdot4)\) is equal to the reciprocal of the observed Fisher information \(I(x)\).”
Observed information = negative Hessian [:224].
“the observed Fisher information \(I(x)=-\ddot{l}_{\theta}\)”
(Kenward and Molenberghs 1998). Invoked by Section 2.4, D6, and Appendix C: under MAR the unconditional expected information differs from the naive (MCAR) form, the mean–covariance orthogonality is lost, and the observed information should be used.
Under MAR the response pattern is not ancillary [:89].
“under the MAR mechanism \(\mathbf{r}\) is not an ancillary statistic for \(\boldsymbol{\theta}\) … we are not justified in restricting the sample space from that associated with the pair ( \(\mathbf{Y}, \mathbf{R}\) ).”
The naive (MCAR) expected information is biased [:60].
“The expected information in these papers is wrong because the expectation is taken under MCAR, in which the missing value mechanism is independent of the distribution of the outcome data.”
MVN: the mean–covariance cross term, vanishing iff MCAR [:374].
“In contrast to the naive information these cross-terms do not all vanish, and the orthogonality of mean and variance-covariance parameters is lost under the MAR mechanism. … the necessary and sufficient condition for the terms in (9) and (10) to be equal to zero is that \(\alpha_{0}=\alpha_{1}=0\), the condition defining, as expected, an MCAR mechanism.”
Recommendation: use the observed information [:38, :544].
(abstract) “it is shown that the classical expected information matrix is biased and the use of the observed information matrix is recommended.” (§5) “provided the observed information matrix is used, conventional likelihood based frequentist inference is applicable in the MAR setting.”
(Pawitan 2013). Invoked by Section 2.4 and the Section 4 notation note: the score-variance \(\mathcal{J}\) vs expected-Hessian \(\mathcal{I}\) distinction, information equality under correct specification, and the TIC penalty \(\operatorname{tr}(\mathcal{J}\mathcal{I}^{-1})\).
\(\mathcal{J}\) (score variance), \(\mathcal{I}\) (expected Hessian); equal under the true model [:14618, :14624, :14627].
\(\mathcal{J} \equiv E\left(\partial\log f_{\theta}/\partial\theta\right)\left(\partial\log f_{\theta}/\partial\theta'\right)\big|_{\theta_0}\) [
:14618]; \(\mathcal{I} \equiv -E\,\partial^{2}\log f_{\theta}/\partial\theta\partial\theta'\big|_{\theta_0}\) [:14624]; “These two matrices are the same if \(f_{\theta_{0}}(x)\) is the true model.” [:14627]
Sandwich variance (Theorem 13.2) and information equality [:14632, :14699].
\(\sqrt{n}(\widehat{\theta}-\theta_{0}) \xrightarrow{d} N(0, \mathcal{I}^{-1}\mathcal{J}\mathcal{I}^{-1})\) [
:14632]; “Assuming the model is correct we get \(\mathcal{J}=\mathcal{I}\), and the formula reduces to the usual inverse Fisher information \(\mathcal{I}^{-1}\) … the ‘naive’ variance formula.” [:14699]
Trace corollary \(\operatorname{tr}(\mathcal{J}\mathcal{I}^{-1})\) (the TIC penalty) [:14707]; the §13.6 AIC derivation restates \(\mathcal{J}_{k}\) at [:14977].
\(E\left\{n(\widehat{\theta}-\theta_{0})'\mathcal{I}(\widehat{\theta}-\theta_{0})\right\} \to \operatorname{trace}\left(\mathcal{I}\mathcal{I}^{-1}\mathcal{J}\mathcal{I}^{-1}\right) = \operatorname{trace}\left(\mathcal{J}\mathcal{I}^{-1}\right)\)
(Konishi and Kitagawa 1996). Invoked by Section 2.4: the TIC/GIC generalization of AIC’s penalty under misspecification — the misspecification face of the \(\mathcal{J}\)-vs-\(\mathcal{I}\) distinction, contrasted with our correctly-specified case.
Theorem 2.1 — bias under possible misspecification [:80–:90].
“Suppose that the parametric family of distributions does not necessarily contain the true distribution … the asymptotic bias of the log likelihood … \(E_{G}\left\{\eta(X_{n};\hat{G})-\eta(X_{n};G)\right\}=\frac{1}{n} b_{1}(G)+o(1/n)\)” [
:80–:84], with \(b_{1}(G)=\operatorname{tr}\left\{\int T^{(1)}(z;G)\,\partial\log f(z\mid\theta)/\partial\theta'\big|_{T(G)}\,dG(z)\right\}\) [:90].
GIC = bias-corrected log-likelihood (a trace penalty) [:100–:101].
\(\operatorname{GIC}(X_{n};\hat{G}) := -2\sum_{\alpha}\log f(X_{\alpha}\mid\hat{\theta}) + \frac{2}{n}\sum_{\alpha}\operatorname{tr}\left\{T^{(1)}(X_{\alpha};\hat{G})\,\partial\log f(X_{\alpha}\mid\theta)/\partial\theta'\big|_{\hat{\theta}}\right\}\)
(Wilks 1938). Invoked by Section 5 (the complete-data \(\chi^2_{q_{\mathrm{d}}}\) reference for \(\hat d_{\text{com}}\)).
Theorem [:99].
“…then when the hypothesis \(H\) is true that \(\theta_{i}=\theta_{0 i}, i=m+1, m+2, \cdots h\), the distribution of \(-2 \log \lambda\), where \(\lambda\) is given by (2) is, except for terms of order \(1 / \sqrt{n}\), distributed like \(\chi^{2}\) with \(h-m\) degrees of freedom.”
(Chan and Meng 2022). Invoked by Section 5.1 (the averaged-likelihood maximizer construction of \(\hat\psi^*, \hat\psi_0^*\), of which \(\bar Q_\infty\)-maximization is the \(M\to\infty\) limit).
The averaged log-likelihood and its constrained/unconstrained maximizers [:185, :191, :197].
\(\bar{L}(\psi)=\frac{1}{m} \sum_{\ell=1}^{m} L^{(\ell)}(\psi)\) [
:185]; \(\widehat{\psi}_{0}^{*}=\arg \max_{\psi \in \Psi: \theta(\psi)=\theta_{0}} \bar{L}(\psi), \quad \widehat{\psi}^{*}=\arg \max _{\psi \in \Psi} \bar{L}(\psi)\) [:191]; \(\hat{d}_{\mathrm{L}}=2\{\bar{L}(\widehat{\psi}^{*})-\bar{L}(\widehat{\psi}_{0}^{*})\}\) [:197].
(Chan 2022). Invoked by Section 5.6 (the weighted-\(\chi^2\) reference whose mean absorbs \(\operatorname{tr}(\mathrm{RIV}_\perp)\) at the null).
FMI/OMI eigenvalues [:79].
“Let the eigenvalues of \(F\) be \(f_{1} \geq \cdots \geq f_{k}\), and denote \(r_{j}:=f_{j} /(1-f_{j})\) for each \(j\). The values \(f_{j}\) and \(r_{j}\) are known as the fraction of missing information (FMI) and the odds of missing information (OMI), respectively.”
The fixed-\(m\) limit law, eq. (1.8) [:82], and its \(m\to\infty\) form, eq. (1.10) [:117].
\(\widetilde{D}_{\aleph} \Rightarrow \mathbb{D}\) [eq. (1.8),
:82]; \(\quad \mathbb{D}_{\infty}:=\frac{1}{k(1+\mu_{r})} \sum_{j=1}^{k}(1+r_{j}) G_{j}\), where \(\mu_{r}:=\frac{1}{k} \sum_{j=1}^{k} r_{j}\) [eq. (1.10),:117], with \(G_j \sim \chi^2_1\) independent.
The simulated (Monte-Carlo) null reference, eq. (5.1) [:401-:407]. Invoked by Section 6.3 (the pre-registered fallback for the equating reference).
“We propose approximating the limiting null distribution \(\mathbb{D}\) by substituting \(r_{j}=\widehat r_{j}\) into (1.8). Although it is not a named distribution, we can easily compute its quantile via Monte Carlo methods. Precisely, we first generate \(G_{j}^{(\iota)} \sim \chi_{1}^{2}\) and \(H_{j}^{(\iota)} \sim \chi_{m-1}^{2}/(m-1)\) independently …” [
:401]; \(\widehat{\mathbb{D}}^{(\iota)}:=\frac{\frac{1}{k} \sum_{j=1}^{k}\{1+(1+\frac{1}{m}) \widehat{r}_{j}\} G_{j}^{(\iota)}}{1+\frac{1}{k} \sum_{j=1}^{k}(1+\frac{1}{m}) \widehat{r}_{j} H_{j}^{(\iota)}}\) [eq. (5.1),:404]; step-by-step form in Algorithm 2.1 [:177-:179].
(Shimodaira and Maeda 2018). Invoked by Section 1.4 and Section 6.1 (the \(\mathrm{AIC}_{x;y}\) criterion whose \(\operatorname{tr}(\mathrm{RIV})\) surcharge the corrected criterion Equation 14 reproduces; the halved-penalty relation to PDIO/AICcd).
The trace decomposition (missing-information surcharge) [:59].
“Let … \(I_{z \mid y}(\boldsymbol{\theta})=I_{x}(\boldsymbol{\theta})-I_{y}(\boldsymbol{\theta})\) be the Fisher information matrix for \(p_{z \mid y}\). Since \(I_{z \mid y}(\boldsymbol{\theta})\) is nonnegative definite, we have \(\operatorname{tr}(I_{x}(\boldsymbol{\theta}) I_{y}(\boldsymbol{\theta})^{-1}) = d+\operatorname{tr}(I_{z \mid y}(\boldsymbol{\theta}) I_{y}(\boldsymbol{\theta})^{-1}) \geq d\).”
\(\mathrm{AIC}_{cd}\) as reproduced there, eq. (4) [:68]; \(\mathrm{AIC}_{x;y}\), eq. (8) [:80]; the halved surcharge, eq. (9) [:86-:89].
\(\mathrm{AIC}_{c d}=-2 Q(\hat{\boldsymbol{\theta}}_{y} ; \hat{\boldsymbol{\theta}}_{y})+2 \operatorname{tr}(I_{x}(\hat{\boldsymbol{\theta}}_{y}) I_{y}(\hat{\boldsymbol{\theta}}_{y})^{-1})\) [
:68]; \(\quad \operatorname{AIC}_{x ; y}=-2 \ell_{y}(\hat{\boldsymbol{\theta}}_{y})+d+\operatorname{tr}(I_{x}(\hat{\boldsymbol{\theta}}_{y}) I_{y}(\hat{\boldsymbol{\theta}}_{y})^{-1})\) [:80]; \(\quad \operatorname{AIC}_{x ; y}-\mathrm{AIC}=\operatorname{tr}(I_{z \mid y}(\hat{\boldsymbol{\theta}}_{y}) I_{y}(\hat{\boldsymbol{\theta}}_{y})^{-1})\) [:86]; “The additional penalty is only half the value of that in PDIO.” [:89]
(Cavanaugh and Shumway 1998). Invoked by Section 1.4 and Section 6.1 (the \(\mathrm{AIC}_{cd}\) criterion with the doubled \(2\operatorname{tr}(\mathrm{RIV})\) surcharge on the EM deviance).
\(\mathrm{AIC}_{cd}\), eq. (3.20) [:302], and the \(2d + 2\operatorname{tr}\) split, eq. (4.4) [:351-:355].
\(\mathrm{AICcd}=-2 Q(\hat{\boldsymbol{\theta}} \mid \hat{\boldsymbol{\theta}})+2 \operatorname{trace}\{\boldsymbol{I}_{o c}(\hat{\boldsymbol{\theta}} \mid \boldsymbol{Y}_{\mathrm{obs}}) \boldsymbol{I}_{o}^{-1}(\hat{\boldsymbol{\theta}} \mid \boldsymbol{Y}_{\mathrm{obs}})\}\) [
:302]; penalty \(= 2 d+2 \operatorname{trace}\{\boldsymbol{D} \boldsymbol{M}(\boldsymbol{I}-\boldsymbol{D} \boldsymbol{M})^{-1}\}\) [:354], the trace term being “the increase in variance [of \(\hat\theta\)] due to missing information” (Meng–Rubin, quoted at [:360]).
(Satorra and Bentler 2010). Invoked by Section 6.3 (the scaled difference statistic, the first-moment-matching precedent for the two-moment arm).
The scaled difference test statistic [:58-:61] and the Satorra–Bentler (2001) combination [:76-:79].
“Satorra (2000, p. 241) proposed the following scaled test statistic: \(\bar{T}_{d}:=T_{d} / \hat{c}_{d}\), where \(c_{d}:=\frac{1}{m} \operatorname{tr}\{U_{d} \Gamma\}\)” [
:61]; “Compute the scaled difference test statistic as \(\tilde{T}_{d}:=T_{d} / \tilde{c}_{d}\) with \(\tilde{c}_{d}=(r_{0} \hat{c}_{0}-r_{1} \hat{c}_{1}) / m\)” [:79].
(Asparouhov and Muthén 2006). Invoked by Section 6.3 (the mean-and-variance-adjusted difference statistic — two-moment matching via a scale and a Satterthwaite degrees-of-freedom adjustment from the moment inputs \(\operatorname{tr}(M), \operatorname{tr}(M^2)\)).
The adjustment [:48-:68].
\(\bar{T}=\frac{d}{\operatorname{tr}(M)} T\) [
:51]; “The distribution of \(\bar{T}\) is approximately a chi-square distribution with \(d^{\prime}\) degrees of freedom where \(d^{\prime}\) is the integer nearest to \(\frac{(\operatorname{tr}(M))^{2}}{\operatorname{tr}(M^{2})}\)” [:64-:67].
(Kolen and Brennan 2014). Invoked by Section 6.3 (the equipercentile equating function applied with the analytic null law as the source distribution).
The equipercentile equating function, eq. (2.11) [:1259-:1265].
“Braun and Holland (1982) indicated that the following function is an equipercentile equating function when \(X\) and \(Y\) are continuous random variables: \(e_{Y}(x)=G^{-1}[F(x)]\), where \(G^{-1}\) is the inverse of the cumulative distribution function \(G\).”
(Imhof 1961). Invoked by Section 6.3 (the closed-form evaluation of the weighted-\(\chi^2\) distribution function in the equating map).
The inversion integral, eq. (3.2) [:234-:247].
\(P\{Q>x\}=\frac{1}{2}+\frac{1}{\pi} \int_{0}^{\infty} \frac{\sin \theta(u)}{u \varphi(u)} d u\) [eq. (3.2),
:237], with \(\theta(u)=\frac{1}{2} \sum_{1}^{m}[h_{r} \tan ^{-1}(\lambda_{r} u)+\delta_{r}^{2} \lambda_{r} u(1+\lambda_{r}^{2} u^{2})^{-1}]-\frac{1}{2} x u\) and \(\rho(u)=\prod_{1}^{m}(1+\lambda_{r}^{2} u^{2})^{\frac{1}{4} h_{r}} \exp\{\cdots\}\) [:243-:246; the central case \(\delta_r = 0\), \(h_r = 1\) is the one used here].
(Davies 1980). Invoked by Section 6.3 (the standard algorithmic implementation of the weighted-\(\chi^2\) distribution function).
Description and purpose [:44-:56].
“Let \(Q=\sum_{j=1}^{r} \lambda_{j} X_{j}+\sigma X_{0}\), where \(X_{j}\) are independent random variables, \(X_{j}\) having a non-central \(\chi^{2}\) distribution … Then the purpose of this algorithm is to calculate \(\operatorname{pr}(Q<c)\). The algorithm is based on the method of Davis (1973) involving the numerical inversion of the characteristic function.” [
:47-:56]
(Meng and Rubin 1992). Invoked by Section 1.4 and Section 5.1 (the original average-the-LRT-statistics combining rule that \(\hat d_L\) refines).
The combining rule from completed-data LR statistics, eq. (2.9) [:168].
\(D_{m} \bumpeq \hat{D}_{m}^{\prime} \equiv\left(\frac{\bar{d}_{m}^{\prime}}{k}-\frac{m-1}{m+1} r_{m}\right) /\left(1+r_{m}\right)\), “where \(\bar{d}_{m}^{\prime}\) is the average of \(\{d_{* l}^{\prime} ; l=1, \ldots, m\}\)” [
:165-:171].
(Vuong 1989). Invoked by Section 5.4 and Section 5.5 (the fixed-alternative / non-nested regime, where the models’ pseudo-true values are \(O(1)\) apart).
Pseudo-true values [:81].
“The value \(\theta_{*}\) is called the pseudo-true value of \(\theta\) for the conditional model \(\boldsymbol{F}_{\theta}\) (see, e.g., Sawa (1978)). Similarly, \(\gamma_{*}\) denotes the pseudo-true value of \(\gamma\) for the conditional model \(\boldsymbol{G}_{\gamma}\).”
Scope of the framework (nested, non-nested, overlapping; LR limit laws) [:40, :55].
“The tests are directional and are derived successively for the cases where the competing models are non-nested, overlapping, or nested and whether both, one, or neither is misspecified. As a prerequisite, we fully characterize the asymptotic distribution of the likelihood ratio statistic under the most general conditions.” [
:40]
(Cox and Snell 1968). Invoked by the A4 lemma (Section 4.4, the non-monotone \(b_\Sigma\)) and by Section 4.6 (the canonical-link MLE-bias input to Lemma P): the general second-order bias of the maximum-likelihood estimate, obtained from the second-order expansion of the score equation.
The second-order score equation [:180].
“To obtain a more refined answer, we replace (10) by the second-order equation \(L^{\prime}(\beta)+(\widehat{\beta}-\beta) L^{\prime \prime}(\beta)+\frac{1}{2}(\widehat{\beta}-\beta)^{2} L^{\prime \prime \prime}(\beta)=0 .\)”
Scalar bias, order \(n^{-1}\) [:229].
“\(b \equiv E(\hat{\beta}-\beta)=\frac{1}{2 I^{2}}(K+2 J)\) which is of order \(n^{-1}\).”
Vector case, their eq. (20) [:269].
“a set of simultaneous linear equations for the biases, with solution \(b_{s} \equiv E\left(\hat{\beta}_{s}-\beta_{s}\right)=\frac{1}{2} I^{r s} I^{t u}\left(K_{r t u}+2 J_{t, r u}\right)\)”
The binomial and logit instantiations (their §§8–9) [:715, :748].
“from which we obtain the biases \(b_{s}=-\frac{1}{2} I^{r s} I^{t u} \Sigma \frac{m_{j}}{\theta_{j}\left(1-\theta_{j}\right)} \frac{\partial \theta_{j}}{\partial \beta_{r}} \frac{\partial^{2} \theta_{j}}{\partial \beta_{t} \partial \beta_{u}},\)” [
:715]
“and substitution into (44) gives \(b_{s}=-4 I^{r s} I^{t u} \sum m_{j} \theta_{j}\left(1-\theta_{j}\right)\left(1-2 \theta_{j}\right) l_{j r} l_{j t} l_{j u}\)” [
:748]
(Firth 1993). Invoked by Section 4.6: the direction of the logistic-regression MLE bias (away from \(\beta=0\)), matching the inline intercept-only \(\sinh(\eta)/n\) value; the paper also attributes the first-order bias formula to Cox & Snell.
Direction: away from zero [:211].
“maximum likelihood estimates of the canonical parameter \(\beta\) are found to be biased away from the point \(\beta=0\). Bias correction, therefore, requires some degree of ‘shrinkage’ of \(\hat{\beta}\) towards this point.”
The first-order bias attributed to Cox & Snell [:137, :140].
“The term \(-n^{-1} \kappa^{r, s} \kappa^{t, u}\left(\kappa_{s, t, u}+\kappa_{s, t u}\right) / 2=n^{-1} b_{1}^{r}(\theta)\) is the first-order bias of \(\hat{\theta}\), for example Cox & Snell (1968).”
The margin [mc] flags throughout the preceding sections link to the executable listings in this appendix. The code runs inline at render time (SymPy + NumPy); every identity is checked by reducing lhs - rhs to 0 (symbolic) or asserting a match to the R code verification/00-setup.R (numeric, to machine precision). This is the algebraic companion to the Monte Carlo study (verification/W1-bias-decomposition.R), which confirms \(\mathbb{E}[T]\approx+\tfrac12\operatorname{tr}(\mathrm{RIV})\) at many concrete configurations; here we check the chain of identities that produces it. The symbolic trace identities and the MVN covariance-block (Magnus duplication-matrix) inversion are independently re-verified by a second computer-algebra system in verification/cas-wolfram/verify_traces.py (Wolfram), so these reductions do not rest on a single CAS implementation.
Scope: the trace bookkeeping that combines the asymptotic second-moment inputs into \(+\tfrac12\operatorname{tr}(\mathrm{RIV})\) is symbolic (exact, parameterization-free); the MVN closed forms for \(I_{\text{com}},I_{\text{obs}},\operatorname{tr}(\mathrm{RIV})\) and the FIML Q-function are reimplemented in NumPy and checked to machine precision against the R code on a fixed replicate; and \(\tfrac12\operatorname{tr}(\mathrm{RIV})\) is checked within Monte-Carlo SE of the W1 cache \(\mathbb{E}[T]\). What Sympy does not prove are the second-moment inputs themselves (the asymptotic covariances L1-L4); those are likelihood asymptotics stated as lemmas in Section 3 and validated end-to-end by the Monte Carlo study.
import sympy as sp
import numpy as np
from sympy import MatMul, MatAdd, Identity, Trace, Transpose, S
from sympy.matrices.expressions.matexpr import MatrixExpr, MatrixSymbol
from numpy.linalg import inv, slogdet, eigvalsh
from pathlib import Path
import json
print("sympy", sp.__version__, "| numpy", np.__version__)sympy 1.14.0 | numpy 2.3.3MVN running example: \(X \sim \mathcal{N}_p(\mu, \Sigma)\), parameter \(\theta = (\mu,\ \operatorname{vech}\Sigma)\) of dimension \(k = p + p(p+1)/2\) (D1, D14). Sympy has no vech/duplication operator, so we build \(D_p\) explicitly, mirroring 00-setup.R (column-major lower triangle).
def vech_index_table(p):
m = np.zeros((p, p), int); c = 0
for j in range(p):
for i in range(j, p):
c += 1; m[i, j] = c; m[j, i] = c
return m
def duplication_matrix(p, sym=False):
kk = p * (p + 1) // 2
D = sp.zeros(p * p, kk) if sym else np.zeros((p * p, kk))
c = 0
for j in range(p):
for i in range(j, p):
D[j * p + i, c] = 1
if i != j:
D[i * p + j, c] = 1
c += 1
return DSympy’s Trace does not auto-apply linearity, does not evaluate Trace(Identity) = n, and does not exploit cyclic invariance Trace(AB) = Trace(BA); worse, sympy.simplify on a Trace hits a signsimp/cancel bug and raises. The helper below (i) expands matrix products inside each trace, (ii) splits Trace over sums, (iii) evaluates Trace(c*I) = c*n, and (iv) canonicalizes cyclic factor order so cyclically-equal traces collapse. Every symbolic identity is checked after trace_normal.
def symmetric_reduce(expr):
"""Treat every bare MatrixSymbol as symmetric: Transpose(X) -> X."""
return expr.replace(
lambda x: isinstance(x, Transpose) and isinstance(x.arg, MatrixSymbol),
lambda x: x.arg)
def trace_normal(expr, n):
"""Normal form for scalar Trace(MatMul) expressions (avoids sympy.simplify)."""
expr = sp.expand(symmetric_reduce(sp.expand(expr, trace=True)), trace=True)
def handle(t):
m = symmetric_reduce(sp.expand(t.arg))
if isinstance(m, MatAdd):
return sum(handle(Trace(a)) for a in m.args)
coeff = S.One
if isinstance(m, MatMul):
for a in m.args:
if not isinstance(a, MatrixExpr):
coeff *= a
mats = [a for a in m.args
if isinstance(a, MatrixExpr) and not isinstance(a, Identity)]
if not mats:
return coeff * n
rotations = [mats[i:] + mats[:i] for i in range(len(mats))]
best = min(rotations, key=lambda fs: tuple(sp.srepr(x) for x in fs))
return coeff * Trace(MatMul(*best))
if isinstance(m, Identity):
return coeff * n
return coeff * Trace(m)
return expr.replace(lambda x: isinstance(x, Trace), handle)Treat \(I_{\text{com}}, I_{\text{obs}}\) as abstract symmetric positive-definite \(n\times n\) matrices (\(n=k\)). The asymptotic second-moment inputs (Lemmas L1-L4 of Section 3, proved there, not here) feed the trace bookkeeping. First, tr(RIV) (Equation 3) in normal form:
n = sp.Symbol('n', positive=True, integer=True)
I_com = sp.MatrixSymbol('I_com', n, n)
I_obs = sp.MatrixSymbol('I_obs', n, n)
I_mis = I_com - I_obs # I_{mis|obs} := I_com - I_obs
trRIV = trace_normal(Trace(I_obs.inv() * I_mis), n)
assert sp.expand(trace_normal(Trace(I_obs.inv() * I_com), n) - n - trRIV) == 0
print("tr(RIV) =", trRIV) # = Trace(I_obs^{-1} I_com) - ntr(RIV) = -n + Trace(I_obs**(-1)*I_com)Term A \(= +\operatorname{tr}(\mathrm{RIV})\) — plug L1 into L3, then cyclic invariance (Equation 8):
A = trace_normal(Trace(I_mis * I_obs.inv()), n) # (L3) with (L1)
assert sp.expand(A - trRIV) == 0
print("A =", A, " -> A - tr(RIV) =", sp.expand(A - trRIV))A = -n + Trace(I_obs**(-1)*I_com) -> A - tr(RIV) = 0Term B \(= -\tfrac12\operatorname{tr}(\mathrm{RIV})\) — plug L2 into L4; the \(n\) (\(=k\)) terms cancel, leaving no dimension constant (Equation 7; this is the cancellation the Term B / L2 disambiguation in the body warns about):
B = trace_normal(sp.Rational(-1, 2) * Trace(I_com * (I_obs.inv() - I_com.inv())), n)
assert sp.expand(B + trRIV / 2) == 0
print("B =", B, " -> B + tr(RIV)/2 =", sp.expand(B + trRIV / 2))B = n/2 - Trace(I_obs**(-1)*I_com)/2 -> B + tr(RIV)/2 = 0Combine \(T = A + B = +\tfrac12\operatorname{tr}(\mathrm{RIV})\) (Equation 9, Equation 6):
T = trace_normal(A, n) + B
assert sp.expand(T - trRIV / 2) == 0
print("T = A+B =", sp.expand(T))
print("T - tr(RIV)/2 =", sp.expand(T - trRIV / 2), " <- exact zero")T = A+B = -n/2 + Trace(I_obs**(-1)*I_com)/2
T - tr(RIV)/2 = 0 <- exact zeroBlock structure of \(I_{\text{com}}\) (D15): block-diagonal on \((\mu,\operatorname{vech}\Sigma)\) with \(I_\mu = N\Sigma^{-1}\) and \(I_{\operatorname{vech}\Sigma} = \tfrac{N}{2}D_p^\top(\Sigma^{-1}\otimes\Sigma^{-1})D_p\).
s11, s12, s22, Nsym = sp.symbols('sigma11 sigma12 sigma22 N', positive=True)
Sig2 = sp.Matrix([[s11, s12], [s12, s22]])
Sig2_inv = Sig2.inv()
D2 = duplication_matrix(2, sym=True)
I_mu_2 = Nsym * Sig2_inv
I_sig_2 = sp.Rational(1, 2) * Nsym * D2.T * sp.Matrix(
sp.kronecker_product(Sig2_inv, Sig2_inv)) * D2
I_sig_2 = sp.simplify(I_sig_2) # safe here: ordinary Matrix, not Trace
I_mu_2, I_sig_2(Matrix([
[ N*sigma22/(sigma11*sigma22 - sigma12**2), -N*sigma12/(sigma11*sigma22 - sigma12**2)],
[-N*sigma12/(sigma11*sigma22 - sigma12**2), N*sigma11/(sigma11*sigma22 - sigma12**2)]]),
Matrix([
[ N*sigma22**2/(2*(sigma11*sigma22 - sigma12**2)**2), -N*sigma12*sigma22/(sigma11*sigma22 - sigma12**2)**2, N*sigma12**2/(2*(sigma11*sigma22 - sigma12**2)**2)],
[-N*sigma12*sigma22/(sigma11*sigma22 - sigma12**2)**2, N*(sigma11*sigma22 + sigma12**2)/(sigma11*sigma22 - sigma12**2)**2, -N*sigma11*sigma12/(sigma11*sigma22 - sigma12**2)**2],
[ N*sigma12**2/(2*(sigma11*sigma22 - sigma12**2)**2), -N*sigma11*sigma12/(sigma11*sigma22 - sigma12**2)**2, N*sigma11**2/(2*(sigma11*sigma22 - sigma12**2)**2)]]))These mirror fisher_info_com_mvn, fisher_info_obs_mvn, loglik_mvn, and barQ_fiml_analytic in verification/00-setup.R. Matching the R outputs on a fixed replicate is the gate that the symbolic construction equals the code the W1 simulations actually ran.
def fisher_com(mu, Sig, N):
p = len(mu); Si = inv(Sig); D = duplication_matrix(p)
kk = p + p * (p + 1) // 2; I = np.zeros((kk, kk))
I[:p, :p] = N * Si
I[p:, p:] = (N / 2) * D.T @ np.kron(Si, Si) @ D
return I
def fisher_obs(Sig, R):
p = Sig.shape[0]; kk = p + p * (p + 1) // 2
vi = vech_index_table(p); I = np.zeros((kk, kk))
pats = {}
for r in range(R.shape[0]):
key = tuple(R[r].astype(int)); pats[key] = pats.get(key, 0) + 1
for key, npat in pats.items():
Oi = [i for i in range(p) if key[i] == 0]; q = len(Oi)
if q == 0:
continue
Sooi = inv(Sig[np.ix_(Oi, Oi)])
for a in range(q):
for b in range(q):
I[Oi[a], Oi[b]] += npat * Sooi[a, b]
Dq = duplication_matrix(q)
Isub = 0.5 * Dq.T @ np.kron(Sooi, Sooi) @ Dq
full = np.array([p + vi[Oi[a], Oi[b]] - 1
for b in range(q) for a in range(b, q)])
I[np.ix_(full, full)] += npat * Isub
return I
def loglik(mu, Sig, X):
Si = inv(Sig); _, ld = slogdet(Sig); hl = 0.5 * np.log(2 * np.pi); p = len(mu)
d = X - mu
return float(np.sum(-p * hl - 0.5 * ld - 0.5 * np.einsum('ij,jk,ik->i', d, Si, d)))
def barQ_fiml(mu_e, Sig_e, mu_i, Sig_i, Y, R):
Np, p = Y.shape
Sei = inv(Sig_e); _, ld = slogdet(Sig_e); hl = 0.5 * np.log(2 * np.pi)
total = 0.0; pats = {}
for r in range(Np):
key = tuple(R[r].astype(int)); pats.setdefault(key, []).append(r)
for key, rows in pats.items():
Mi = [i for i in range(p) if key[i] == 1]
Oi = [i for i in range(p) if key[i] == 0]
if not Mi: # fully observed
for r in rows:
d = Y[r] - mu_e
total += -p * hl - 0.5 * ld - 0.5 * d @ Sei @ d
elif not Oi: # all missing
d = mu_i - mu_e; tr = np.trace(Sei @ Sig_i)
total += len(rows) * (-p * hl - 0.5 * ld - 0.5 * d @ Sei @ d - 0.5 * tr)
else: # partial: Schur conditional
Sooi = inv(Sig_i[np.ix_(Oi, Oi)]); Smo = Sig_i[np.ix_(Mi, Oi)]
condV = Sig_i[np.ix_(Mi, Mi)] - Smo @ Sooi @ Smo.T
tr = np.trace(Sei[np.ix_(Mi, Mi)] @ condV)
for r in rows:
cm = mu_i[Mi] + Smo @ Sooi @ (Y[r, Oi] - mu_i[Oi])
Xh = np.zeros(p)
for o in Oi:
Xh[o] = Y[r, o]
for jj, mm in enumerate(Mi):
Xh[mm] = cm[jj]
d = Xh - mu_e
total += -p * hl - 0.5 * ld - 0.5 * d @ Sei @ d - 0.5 * tr
return totalThe fixture proof-verification-fixture.json is one W1 replicate (N=500, non-monotone MAR) plus the W1-sweep-prod cell aggregates, emitted by verification/W1-sympy-fixture.R.
cands = [Path("proof-verification-fixture.json"),
Path("manuscript/proof-verification-fixture.json")]
fix = next(p for p in cands if p.exists())
J = json.load(open(fix))
def as_float(x):
try:
return float(x)
except (TypeError, ValueError):
return np.nan # jsonlite serializes missing entries of Y as "NA"
p, k, N = J["meta"]["p"], J["meta"]["k"], J["meta"]["N"]
X = np.array(J["data"]["X"], float)
Y = np.array([[as_float(v) for v in row] for row in J["data"]["Y"]])
R = np.array(J["data"]["R"], float)
mu_obs = np.array(J["theta_obs"]["mu"], float)
Sig_obs = np.array(J["theta_obs"]["Sigma"], float)
mu_com = np.array(J["theta_com"]["mu"], float)
Sig_com = np.array(J["theta_com"]["Sigma"], float)
I_com_R = np.array(J["fisher"]["I_com"], float)
I_obs_R = np.array(J["fisher"]["I_obs"], float)
sc = J["scalars"]; cache = J["cache"]
print(f"replicate: p={p} k={k} N={N}, missing entries={int(np.isnan(Y).sum())}")replicate: p=4 k=14 N=500, missing entries=400Fisher information and tr(RIV) match R (D15, Equation 3):
I_com_py = fisher_com(mu_obs, Sig_obs, N)
I_obs_py = fisher_obs(Sig_obs, R)
trRIV_py = float(np.trace(inv(I_obs_py) @ I_com_py) - k)
err_com = np.max(np.abs(I_com_py - I_com_R))
err_obs = np.max(np.abs(I_obs_py - I_obs_R))
print(f"max|I_com_py - I_com_R| = {err_com:.2e}")
print(f"max|I_obs_py - I_obs_R| = {err_obs:.2e}")
print(f"tr(RIV) py = {trRIV_py:.12f} | R = {sc['trRIV_replicate']:.12f}"
f" | diff = {trRIV_py - sc['trRIV_replicate']:.2e}")
assert err_com < 1e-8 and err_obs < 1e-8
assert abs(trRIV_py - sc["trRIV_replicate"]) < 1e-8max|I_com_py - I_com_R| = 2.73e-12
max|I_obs_py - I_obs_R| = 1.36e-12
tr(RIV) py = 7.600274370088 | R = 7.600274370088 | diff = 7.11e-15Information identity and \(I_{\text{mis}\mid\text{obs}}\succeq 0\) (D7, Equation 1):
I_mis_py = I_com_py - I_obs_py
print("max|I_com - (I_obs + I_mis)| =",
f"{np.max(np.abs(I_com_py - (I_obs_py + I_mis_py))):.2e}")
print("min eigenvalue of I_mis|obs =", f"{eigvalsh(I_mis_py).min():.4e} (>= 0)")
assert eigvalsh(I_mis_py).min() > -1e-8max|I_com - (I_obs + I_mis)| = 4.55e-13
min eigenvalue of I_mis|obs = -2.5625e-13 (>= 0)FIML Q-function and the A / B / T decomposition match R (D9, D10-D12):
barL_py = barQ_fiml(mu_obs, Sig_obs, mu_obs, Sig_obs, Y, R)
ell_com_at_com_py = loglik(mu_com, Sig_com, X)
ell_com_at_obs_py = loglik(mu_obs, Sig_obs, X)
A_py = barL_py - ell_com_at_obs_py
B_py = ell_com_at_obs_py - ell_com_at_com_py
T_py = barL_py - ell_com_at_com_py
for name, py, r in [("barL_obs", barL_py, sc["barL_obs"]),
("A_r", A_py, sc["A_r"]),
("B_r", B_py, sc["B_r"]),
("T_r", T_py, sc["T_r"])]:
print(f"{name:9s} py = {py:14.8f} | R = {r:14.8f} | diff = {py - r:.2e}")
assert abs(py - r) < 1e-6barL_obs py = -2421.57536151 | R = -2421.57536151 | diff = 0.00e+00
A_r py = 30.60215855 | R = 30.60215855 | diff = 0.00e+00
B_r py = -5.01518599 | R = -5.01518599 | diff = 0.00e+00
T_r py = 25.58697257 | R = 25.58697257 | diff = 0.00e+00\(\tfrac12\operatorname{tr}(\mathrm{RIV})\) within Monte-Carlo SE of \(\mathbb{E}[T]\) (Equation 6 against the W1 production cache):
half_pred_rep = 0.5 * trRIV_py # from this replicate's theta_obs
half_pred_mc = cache["half_trRIV_pred"] # 0.5 * MC-mean of analytic tr(RIV)
T_hat, T_mcse = cache["T_hat"], cache["T_mcse"]
lo, hi = T_hat - 1.96 * T_mcse, T_hat + 1.96 * T_mcse
print(f"cache E[T] = {T_hat:.4f} +/- {T_mcse:.4f} (MCSE)")
print(f"95% MCSE band = [{lo:.4f}, {hi:.4f}]")
print(f"1/2 tr(RIV) [MC mean] = {half_pred_mc:.4f}")
print(f"1/2 tr(RIV) [replicate] = {half_pred_rep:.4f}")
print(f"cache E[A] = {cache['A_hat']:.3f} +/- {cache['A_mcse']:.3f}"
f" vs tr(RIV)_MC = {cache['trRIV_mean']:.3f}")
print(f"cache E[B] = {cache['B_hat']:.3f} +/- {cache['B_mcse']:.3f}"
f" vs -1/2 tr(RIV)_MC = {-half_pred_mc:.3f}")
assert lo <= half_pred_mc <= hi
print("\nGATE PASSED: 1/2 tr(RIV) lies within the MCSE band of E[T].")cache E[T] = 3.6701 +/- 0.4525 (MCSE)
95% MCSE band = [2.7832, 4.5570]
1/2 tr(RIV) [MC mean] = 3.8367
1/2 tr(RIV) [replicate] = 3.8001
cache E[A] = 7.856 +/- 0.474 vs tr(RIV)_MC = 7.673
cache E[B] = -4.185 +/- 0.058 vs -1/2 tr(RIV)_MC = -3.837
GATE PASSED: 1/2 tr(RIV) lies within the MCSE band of E[T].The \(+\operatorname{tr}(\mathrm{RIV})\) of Term A (Section 4.4) enters the FIML Q-function through the per-observation conditional-variance penalty \(-\tfrac12\operatorname{tr}(\Sigma_{\text{eval}}^{-1}\bar V_i)\), where for a partial pattern \(V = \Sigma_{MM} - \Sigma_{MO}\Sigma_{OO}^{-1}\Sigma_{OM}\) is the Schur complement (Section 3.2). Shown symbolically for a 3-variable pattern, \(X_1\) missing:
syms = sp.symbols('s11 s12 s13 s22 s23 s33', positive=True)
s11, s12, s13, s22, s23, s33 = syms
Sig3 = sp.Matrix([[s11, s12, s13], [s12, s22, s23], [s13, s23, s33]])
Mi, Oi = [0], [1, 2]
Smm = Sig3[Mi, Mi]; Smo = Sig3[Mi, Oi]; Soo = Sig3[Oi, Oi]
V = sp.simplify(Smm - Smo * Soo.inv() * Smo.T) # Schur complement, scalar here
print("conditional Var(X_1 | X_2, X_3) =")
sp.pprint(V[0, 0])conditional Var(X_1 | X_2, X_3) =
⎛ 2⎞
s₁₁⋅⎝s₂₂⋅s₃₃ - s₂₃ ⎠ - s₁₂⋅(s₁₂⋅s₃₃ - s₁₃⋅s₂₃) + s₁₃⋅(s₁₂⋅s₂₃ - s₁₃⋅s₂₂)
────────────────────────────────────────────────────────────────────────
2
s₂₂⋅s₃₃ - s₂₃ The pattern-sum-to-\(\operatorname{tr}(\mathrm{RIV})\) reduction is a conditional-expectation calculation Sympy does not carry out cleanly; it is verified numerically (A_py matches R to machine precision in Listing 7, and the cache confirms \(\mathbb{E}[A]\approx\operatorname{tr}(\mathrm{RIV})\)).
Appendix C claims the observed information carries a mean \(\times\) (co)variance cross block under MAR that the naive block-diagonal form omits, and that it is the selection-induced effect of Kenward and Molenberghs (1998). Both halves are checked inline here. First the observed information in closed form (the realized negative Hessian, with per-pattern residual \(r_P\) and cross-product \(M_P\)), validated against a finite-difference Hessian of \(\ell_{\text{obs}}\) on the fixture replicate:
def loglik_obs(mu, Sig, Y, R):
p = len(mu); hl = 0.5 * np.log(2 * np.pi); total = 0.0; pats = {}
for r in range(Y.shape[0]):
key = tuple(R[r].astype(int)); pats.setdefault(key, []).append(r)
for key, rows in pats.items():
Oi = [i for i in range(p) if key[i] == 0]; q = len(Oi)
if q == 0:
continue
S = Sig[np.ix_(Oi, Oi)]; Si = inv(S); _, ld = slogdet(S)
for r in rows:
d = Y[np.ix_([r], Oi)][0] - mu[Oi]
total += -q * hl - 0.5 * ld - 0.5 * (d @ Si @ d)
return total
def observed_info_obs(mu, Sig, Y, R):
# Realized negative Hessian of l_obs -- the OBSERVED information (mirrors
# observed_info_obs_mvn in 00-setup.R): mean block n_P S^-1; mean x vech
# (r_P' S^-1 (x) S^-1) D_q (propto r_P); vech x vech with M_P.
p = len(mu); kk = p + p * (p + 1) // 2; vi = vech_index_table(p)
J = np.zeros((kk, kk)); pats = {}
for r in range(Y.shape[0]):
key = tuple(R[r].astype(int)); pats.setdefault(key, []).append(r)
for key, rows in pats.items():
Oi = [i for i in range(p) if key[i] == 0]; q = len(Oi)
if q == 0:
continue
S = Sig[np.ix_(Oi, Oi)]; Si = inv(S); nP = len(rows)
d = Y[np.ix_(rows, Oi)] - mu[Oi] # nP x q residuals
rP = d.sum(axis=0); MP = d.T @ d # q ; q x q
Dq = duplication_matrix(q)
full = np.array([p + vi[Oi[a], Oi[b]] - 1
for b in range(q) for a in range(b, q)])
J[np.ix_(Oi, Oi)] += nP * Si # mean block
cross = np.kron((rP @ Si)[None, :], Si) @ Dq # q x |vech|
J[np.ix_(Oi, full)] += cross
J[np.ix_(full, Oi)] += cross.T
SiMSi = Si @ MP @ Si
vv = Dq.T @ (-0.5 * nP * np.kron(Si, Si)
+ 0.5 * (np.kron(SiMSi, Si) + np.kron(Si, SiMSi))) @ Dq
J[np.ix_(full, full)] += vv
return Jdef pack(mu, Sig):
return np.concatenate([mu, [Sig[i, j] for j in range(p) for i in range(j, p)]])
def unpack(th):
mu = th[:p]; Sig = np.zeros((p, p)); c = p
for j in range(p):
for i in range(j, p):
Sig[i, j] = Sig[j, i] = th[c]; c += 1
return mu, Sig
th0 = pack(mu_obs, Sig_obs); hh = 1e-4
def fobs(th):
mu_, Sig_ = unpack(th); return loglik_obs(mu_, Sig_, Y, R)
Hfd = np.zeros((k, k))
for a in range(k):
for b in range(a, k):
ea = np.zeros(k); ea[a] = hh; eb = np.zeros(k); eb[b] = hh
Hfd[a, b] = Hfd[b, a] = (fobs(th0 + ea + eb) - fobs(th0 + ea - eb)
- fobs(th0 - ea + eb) + fobs(th0 - ea - eb)) / (4 * hh * hh)
O_obs_py = observed_info_obs(mu_obs, Sig_obs, Y, R)
rel_fd = np.max(np.abs(O_obs_py - (-Hfd))) / np.max(np.abs(O_obs_py))
cross_obs = np.max(np.abs(O_obs_py[:p, p:])) # observed: mean x vech block
cross_naive = np.max(np.abs(I_obs_py[:p, p:])) # naive fisher_obs: block-diagonal
scale = np.max(np.abs(np.diag(O_obs_py)))
trRIV_obs = float(np.trace(inv(O_obs_py) @ I_com_py) - k)
print(f"observed info vs -Hessian(l_obs): max rel diff = {rel_fd:.2e}")
print(f"cross block observed = {cross_obs:8.3f} ({cross_obs/scale:.3f} of diag)"
f" | naive = {cross_naive:.2e}")
print(f"tr(RIV) observed = {trRIV_obs:.4f} | naive = {trRIV_py:.4f}")
assert rel_fd < 1e-3 # observed info IS the realized Hessian
assert cross_naive < 1e-9 # naive form is block-diagonal
assert cross_obs > 0.01 * scale # observed form has a real cross blockobserved info vs -Hessian(l_obs): max rel diff = 1.48e-06
cross block observed = 158.300 (0.063 of diag) | naive = 0.00e+00
tr(RIV) observed = 8.2460 | naive = 7.6003The cross block is non-zero because, under MAR, selection on the observed values shifts each pattern’s observed mean off the marginal, so \(\mathbb{E}[r_P]\neq0\). The selected-Gaussian closed form (Appendix C) — non-zero under MAR, exactly zero under MCAR — against a seeded simulation:
import math
rng = np.random.default_rng(20260531)
mu0 = np.zeros(3)
Sig0 = np.array([[1.0, 0.3, 0.5], [0.3, 1.0, 0.2], [0.5, 0.2, 1.0]])
a_sel, b_sel, Nbig = -0.5, 0.4, 2_000_000
Xs = rng.multivariate_normal(mu0, Sig0, size=Nbig)
def Phi(z): return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0)))
def phi(z): return math.exp(-0.5 * z * z) / math.sqrt(2.0 * math.pi)
for label, b in [("MAR ", b_sel), ("MCAR", 0.0)]:
V = rng.standard_normal(Nbig) # R=1 (coord 0 missing) iff V - b*X[:,2] < a
sel = (V - b * Xs[:, 2]) < a_sel
zc = a_sel / math.sqrt(1 + b * b); lam = phi(zc) / Phi(zc)
pred = (b * Sig0[:, 2] / math.sqrt(1 + b * b)) * lam # E[X_k - mu_k | selected]
emp = Xs[sel].mean(axis=0) - mu0
print(f"{label}: shift empirical = {np.round(emp, 4)} | closed form = {np.round(pred, 4)}")
if b == 0.0:
assert np.max(np.abs(emp)) < 5e-3 # MCAR: no selection shift
else:
assert np.max(np.abs(emp - pred)) < 5e-3 # MAR: matches inverse-Mills formMAR : shift empirical = [0.2075 0.0824 0.415 ] | closed form = [0.2071 0.0828 0.4141]
MCAR: shift empirical = [ 0.0007 0.0004 -0.001 ] | closed form = [0. 0. 0.]The trace bookkeeping above used abstract MatrixSymbol objects with no MVN commitment, so the leading-order \(T = +\tfrac12\operatorname{tr}(\mathrm{RIV})\) holds for any regular likelihood given L1-L4 and an estimated scale — the assumption behind Term A’s \(\operatorname{tr}(I_{\text{mis}\mid\text{obs}}V_{\text{obs}})\) input, which needs the entropy-gradient cancellation of Section 4.4 and so requires \(\alpha\neq0\). For a known-scale, location-only fit that input halves and \(T=0\) (scope remark after the theorem). Under MAR, add the design-imbalance \((A)+(C)\) of Section 4.4, which is \(0\) under MCAR. What remains MVN-free to verify is that the score-equation expansion delivers the covariance input L2 (Equation 5): the observed score is the conditional expectation of the complete-data score under MAR, so \(\mathbb{E}[U_{\text{com}}U_{\text{obs}}^\top]=I_{\text{obs}}\).
# (L2) from the efficient-estimator orthogonality, symbolically.
Vobs = I_obs.inv() # Var(theta_obs) (L1)
Vcom = I_com.inv() # Var(theta_com)
Cov_oc = (I_obs.inv() * I_obs * I_com.inv()).doit() # Cov(obs,com) = I_com^{-1}
print("Cov(theta_obs, theta_com) =", Cov_oc, " (= I_com^{-1} = Var(theta_com))")
Vdiff = symmetric_reduce(sp.MatAdd(Vobs, -Cov_oc, -Cov_oc.T, Vcom).doit())
target = (I_obs.inv() - I_com.inv()).doit()
print("Var(diff) - (I_obs^{-1} - I_com^{-1}) =", sp.expand(symmetric_reduce(Vdiff - target)))Cov(theta_obs, theta_com) = I_com**(-1) (= I_com^{-1} = Var(theta_com))
Var(diff) - (I_obs^{-1} - I_com^{-1}) = 0Feeding the symbolic L2 straight back into Term B and recombining — no separately-asserted covariance, a pure consequence of the score relations (L2 Equation 5):
B_from_scores = trace_normal(sp.Rational(-1, 2) * Trace(I_com * Vdiff), n)
A_from_L1 = trace_normal(Trace(I_mis * Vobs), n)
T_B = trace_normal(A_from_L1, n) + B_from_scores
assert sp.expand(T_B - trRIV / 2) == 0
print("Tier B: T = A + B =", sp.expand(T_B), " == tr(RIV)/2 [exact]")Tier B: T = A + B = -n/2 + Trace(I_obs**(-1)*I_com)/2 == tr(RIV)/2 [exact]| Identity | Status | Listing |
|---|---|---|
| \(\operatorname{tr}(\mathrm{RIV}) = \operatorname{tr}(I_{\text{obs}}^{-1}I_{\text{com}}) - k\) | symbolic | Listing 1 |
| Term A \(= +\operatorname{tr}(\mathrm{RIV})\) (given L1, L3) | symbolic | Listing 2 |
| Term B \(= -\tfrac12\operatorname{tr}(\mathrm{RIV})\) (given L2, L4) | symbolic | Listing 3 |
| \(T = A + B = +\tfrac12\operatorname{tr}(\mathrm{RIV})\) | symbolic | Listing 4 |
| MVN \(I_{\text{com}}, I_{\text{obs}}, \operatorname{tr}(\mathrm{RIV})\) | numeric (machine precision) | Listing 5 |
| Information identity, \(I_{\text{mis}\mid\text{obs}}\succeq 0\) | numeric | Listing 6 |
| FIML Q-function, \(A/B/T\) | numeric (machine precision) | Listing 7 |
| \(\tfrac12\operatorname{tr}(\mathrm{RIV})\) vs \(\mathbb{E}[T]\) | within MCSE | Listing 8 |
| Observed info \(=\) realized \(-\)Hessian; cross block (naive omits); larger RIV | numeric (FD, MAR fixture) | Listing 9 |
| MAR cross term \(\mathbb{E}[r_P]\) vs inverse-Mills closed form; MCAR \(=0\) | numeric (seeded sim) | Listing 10 |
| Tier B (abstract, MVN-free): L2 from scores | symbolic | Listing 11 |
The RIV’s denominator is the observed-data information \(I_{\text{obs}}\) (Section 2.4), which under MAR is not the block-diagonal naive form: the unconditional expectation carries a mean–covariance cross term (Kenward and Molenberghs 1998, sec. 3.2), :374, and is estimated by the observed information \(-\partial^2\ell_{\text{obs}}/\partial\theta\partial\theta^\top\). This appendix gives that observed information in closed form for the MVN under missingness — so the RIV is computable analytically, with no imputations and no numerical differentiation — and shows where the cross term comes from.
Observed information. The observed information is the realized negative Hessian \(\mathcal{J}_{\text{obs}}(\theta;Y)=-\partial^2\ell_{\text{obs}}/\partial\theta\,\partial\theta^\top\), the quantity whose expectation is \(I_{\text{obs}}\) (D6). It is the second differential \(-d^2\ell\) before the expectation that yields the Fisher information (Magnus and Neudecker 2019, Ch 15 §3, §6), :1800; the component derivatives (of \(\log|S|\) and of the quadratic form) are (Petersen and Pedersen 2012, sec. 2.4.4), :510, the Gaussian second derivative carrying the residual structure \(S^{-1}ee^\top S^{-1}-S^{-1}\) [:2013].
Per missingness pattern \(P\) (observed set \(O\), \(n_P\) observations, \(S=\Sigma_{OO}\), \(q=|O|\), \(D_q\) the duplication matrix), with residual sum \(r_P=\sum_{i\in P}(y_i^O-\mu_O)\) and cross-product \(M_P=\sum_{i\in P}(y_i^O-\mu_O)(y_i^O-\mu_O)^\top\), the pattern’s contribution to \(\mathcal{J}_{\text{obs}}\) in \((\mu,\operatorname{vech}\Sigma)\) coordinates is, blockwise, \[
\begin{aligned}
\text{mean}:\quad & n_P\,S^{-1} \quad(\text{data-free; equals the expected-info block}),\\
\text{mean}\times\operatorname{vech}:\quad & (r_P^\top S^{-1}\otimes S^{-1})\,D_q \quad(\propto r_P),\\
\operatorname{vech}\times\operatorname{vech}:\quad & D_q^\top\!\Big[-\tfrac{n_P}{2}(S^{-1}\!\otimes S^{-1})+\tfrac12\big((S^{-1}M_PS^{-1}\otimes S^{-1})+(S^{-1}\otimes S^{-1}M_PS^{-1})\big)\Big]D_q,
\end{aligned}
\] summed over patterns and embedded in the full index set exactly as in fisher_info_obs_mvn. Taking expectations, \(\mathbb{E}[M_P]=n_PS\) collapses the \(\operatorname{vech}\) block to \(\tfrac{n_P}{2}D_q^\top(S^{-1}\!\otimes S^{-1})D_q\), and the cross block is governed by \(\mathbb{E}[r_P]\). Setting \(\mathbb{E}[r_P]=0\) zeroes it and recovers the block-diagonal \(I_{\text{obs}}\) of D15 (fisher_info_obs_mvn) — but \(\mathbb{E}[r_P]=0\) holds only under MCAR. [mc] the observed information equals the realized \(-\)Hessian of \(\ell_{\text{obs}}\) (finite differences, Listing 9 here; numerical Hessian in verification/00-test-primitives.R), it carries a cross block the naive form omits, and its RIV is the larger one.
The cross term under MAR. Under MAR the missingness depends on the observed values, so each pattern’s observed mean is selected away from the marginal and \(\mathbb{E}[r_P]\neq 0\). For the running example’s mechanism — coordinate \(j\) made missing with probability \(\Phi(a+b\,y_{\text{sel}})\), equivalently selection on a Gaussian \(S_{\text{sel}}=V-b\,y_{\text{sel}}\) with \(V\sim N(0,1)\) independent — the selected-Gaussian moment gives, per selected coordinate, \[ \mathbb{E}[\,y_k-\mu_k \mid \text{selected}\,] \;=\; \frac{b\,\Sigma_{k,\text{sel}}}{\sqrt{1+b^2}}\;
\lambda\!\Big(\tfrac{a}{\sqrt{1+b^2}}\Big), \qquad \lambda(z)=\phi(z)/\Phi(z) \] (the inverse Mills ratio). The shift is \(\propto b\,(\Sigma\text{-coupling})\) — non-zero under MAR, exactly zero under MCAR (\(b=0\)) — so \(\mathbb{E}[r_P]\neq0\) and the cross block \((\mathbb{E}[r_P]^\top S^{-1}\!\otimes S^{-1})D_q\) survives in the unconditional \(I_{\text{obs}}\). This is the MVN, probit-mechanism instance of Kenward and Molenberghs (1998), §3.2, whose logistic-mechanism bivariate version has the same structure \(\mathbb{E}[y_2-\mu_2\mid r]=\beta\alpha_1\), \(\beta=\sigma_{12}\sigma_{11}^{-1}\), :374. [mc] the closed-form shift matches a seeded simulation under MAR and vanishes under MCAR (Listing 10 here; the project mechanism against \(4\times10^6\) draws, MAR shift \(\approx[0.124,0,0.414,0.207]\), in todo/008); the cross-block contribution to \(\operatorname{tr}(\mathrm{RIV})\) does not decay with \(n\) (todo/008).
The RIV and the IC correction. The (single) RIV is \(\operatorname{tr}(\mathcal{J}_{\text{obs}}^{-1}I_{\text{mis}\mid\text{obs}})\) — the observed information \(\mathcal{J}_{\text{obs}}\) above in place of the naive block-diagonal \(I_{\text{obs}}\), the quantity tr_riv_observed_general computes. The information-criterion correction (manuscript §5) subtracts this per candidate model. For model selection specifically, using the naive block-diagonal information instead is, empirically, close to a wash: the naive-vs-unconditional gap is roughly proportional across candidate models, so it nearly cancels in the \(\arg\min\) (w3-sweep.R); the residual gap to the oracle is finite-sample MI-AIC noise. The correction is not optional, however, for the bias level — the deviance-bias estimate itself and likelihood-ratio test calibration — where the cross term must be carried.