Bernoulli longitudinal imputation in LongMI

Proposed release: 0.2.0a1
Scope: package work only
Status: implemented — this is the design record for the released 0.2.0a1 feature set. The authoritative descriptions of what the software does are the algorithm page and the project status; where this plan and the implementation differ (e.g. the constructor argument is n_quad), the implementation governs.

Purpose. Add principled multiple imputation for one incomplete longitudinal binary response while preserving LongMI's existing contract:

validate data → fit imputation model → generate completed datasets
→ fit the complete-data analysis repeatedly → pool → report validity conditions

This feature is not a study-specific analysis. LongMI must not contain references to cannabis, anxiolytics, ORR, medication classes, or a manuscript-specific estimand.

The plan follows established multiple-imputation theory and recent guidance for longitudinal data, while retaining LongMI's explicit distinction between assumptions, mechanically verified properties, and empirical validation (Rubin, 1976; Rubin, 1987; Wijesuriya et al., 2025).

flowchart LR
    A["Validated binary longitudinal data"] --> B["Fit Bernoulli random-intercept model"]
    B --> C["Draw parameters and participant effects"]
    C --> D["Draw missing 0/1 outcomes"]
    D --> E["Analyze each completed dataset"]
    E --> F["Pool estimates and covariance"]
    F --> G["Report assumptions, diagnostics, and sensitivity"]

1. Release objective

LongMI currently supports continuous and count outcomes. Binary support adds a new outcome contract, a new stochastic backend, new numerical failure modes, and a new compatibility question when a subject-specific imputer is followed by a marginal logistic GEE. Longitudinal MI methods must preserve the repeated-measures structure and the relationships required by the intended substantive analysis (White, Royston, and Wood, 2011; Wijesuriya et al., 2025).

The release should therefore be a new minor-feature alpha:

longmi 0.2.0a1

Supported in 0.2.0a1

  • one incomplete longitudinal binary response;
  • observed outcomes restricted to 0 and 1;
  • fully observed numeric predictors;
  • an explicit participant-by-wave grid;
  • independent participants or declared independent clusters;
  • categorical wave effects and declared predictor-by-wave interactions;
  • a participant random intercept;
  • MAR imputation;
  • logit-scale delta sensitivity for specified MNAR departures;
  • group- and wave-targeted delta rules;
  • arbitrary complete-data analysis through the existing adapter interface;
  • multivariate Rubin pooling (Rubin, 1987; Barnard and Rubin, 1999);
  • a validated Bernoulli-MI-to-logistic-GEE workflow in prespecified simulation settings.

Explicitly outside 0.2.0a1

  • binomial totals such as successes out of n trials;
  • incomplete predictors;
  • multiple jointly incomplete responses;
  • random slopes;
  • structural eligibility or competing events;
  • automatic causal interpretation;
  • automatic resolution of complete or quasi-separation;
  • outcome-scale delta adjustment for binary responses;
  • a claim that every Bernoulli imputer and every logistic GEE are congenial.

2. Public data contract

The user-facing constructor becomes:

from longmi import LongitudinalData

long_data = LongitudinalData(
    frame,
    id_col="participant",
    time_col="wave",
    outcome_col="outcome",
    predictor_cols=("exposure", "age", "baseline_outcome"),
    outcome_type="binary",
    times=(0, 3, 6, 9, 12),
)

For every observed participant-wave response,

\[ Y_{ij}^{\mathrm{obs}} \in \{0,1\}. \]

For every completed dataset,

\[ Y_{ij}^{(m)} \in \{0,1\}. \]

Observed outcomes must remain unchanged:

\[ Y_{ij}^{(m)} = Y_{ij}^{\mathrm{obs}} \qquad \text{whenever } R_{ij}=1. \]

Mechanically enforced invariants

The constructor and completion checks must reject:

  • values outside 0 and 1;
  • nonfinite outcomes;
  • malformed strings or silent categorical recoding;
  • duplicate participant-wave rows;
  • absent scheduled rows when a complete grid is required;
  • empty datasets;
  • missing, nonnumeric, or nonfinite imputation predictors;
  • completed datasets that alter observed outcomes;
  • completed datasets that retain missing values.

Boolean outcomes may be normalized to 0 and 1. Character or categorical outcomes must not be silently converted.

3. Validity conditions: assumptions A1–A8

LongMI's canonical assumptions are defined in docs/theory/assumptions.md. The Bernoulli backend should map each assumption to a concrete package obligation.

Assumption Meaning for the Bernoulli backend Package status Required response when doubtful
A1 — Independent clusters Participants or declared clusters are independent; repeated outcomes within a participant may be dependent. Partly verified. LongMI records the cluster identifier and preserves complete trajectories. Independence itself is substantive. Resample the independent cluster, not rows. If participants are nested in sites and sites induce dependence, the site may be the correct sampling unit.
A2 — MAR Conditional on observed outcomes and predictors, missingness does not depend further on the missing binary response. Declared, not verified. MAR cannot be distinguished from MNAR using the observed data alone. Use logit-scale delta sensitivity; report MAR as an assumption rather than a finding.
A3 — Distinct parameters Outcome-model and missingness-model parameters are variation independent, so the missingness model is ignorable for outcome-model inference under MAR. Declared. The Bernoulli imputer does not fit a shared-parameter dropout model. A shared latent process for response and dropout requires a joint/shared-parameter model outside this release.
A4 — Positivity / empirical support Relevant participant histories must have some observed follow-up support. Entire design strata cannot be structurally unobserved. Partly diagnosed. The backend can summarize events, nonevents, and observed responses by wave and covariate pattern. Treat structural non-observation as a design or eligibility problem, not an imputation problem.
A5 — Adequate imputation model The random-intercept Bernoulli model must represent the conditional means, wave effects, interactions, and dependence needed for the estimand. Not provable from one dataset. Evaluated through diagnostics and simulation. Use predictive checks, richer interactions, alternative dependence structures, and deliberate misspecification simulations.
A6 — Parameter and outcome uncertainty Every imputation must vary both model parameters and missing binary outcomes. Backend verified. Parameter draws and Bernoulli outcome draws occur for every completed dataset. Refuse deterministic probability filling or rounded fitted values.
A7 — Valid complete-data analysis The estimator fitted to each completed dataset must be valid, with an appropriate covariance estimator. Partly verified by adapters. For GEE, use the robust sandwich covariance. Reject nonconverged fits; allow alternate complete-data analyses when needed.
A8 — Compatibility / congeniality The imputation model must preserve all analysis variables, wave terms, exposure-by-wave interactions, and within-participant dependence relevant to the substantive estimand. Conditionally supported. A conditional logistic mixed imputer is not automatically congenial with a marginal logistic GEE. Use categorical wave and declared interactions; evaluate the full MI-to-GEE workflow by simulation; consider bootstrap-then-impute when standard Rubin variance is doubtful (Bartlett and Hughes, 2020).

MAR is never marked as verified. LongMI may verify that observed outcomes were preserved and that both sources of imputation uncertainty were propagated. It cannot verify MAR, the adequacy of the imputation model, or congeniality from observed data alone.

The distinction between MAR ignorability and the truth of MAR follows the classical missing-data framework (Rubin, 1976; Little and Rubin, 2019). The additional qualifications for repeated-imputation inference and estimating-equation analyses are grounded in the congeniality and imputation-estimator literature (Meng, 1994; Wang and Robins, 1998; Robins and Wang, 2000; Yang and Kim, 2016).

What the software should distinguish

Type of statement Examples
Verified by construction binary support; observed outcomes preserved; all eligible missing outcomes filled; parameter draws generated; Bernoulli outcomes drawn rather than rounded
Declared by the analyst or backend MAR; independent sampling unit; adequacy of the fixed-effect structure; compatibility with the analysis model
Evaluated empirically optimizer stability; predictive behavior; bias and coverage in simulation; cross-language numerical or statistical agreement

4. Statistical model

For participant \(i\) and wave \(j\),

\[ Y_{ij}\mid b_i \sim \mathrm{Bernoulli}(p_{ij}), \]

with

\[ \mathrm{logit}(p_{ij}) = x_{ij}^{\top}\alpha + b_i, \qquad b_i \sim N(0,\tau^2). \]

This is a standard subject-specific logistic mixed-model parameterization. Its regression coefficients are conditional on the participant effect and therefore need not equal marginal, population-averaged logistic coefficients (Zeger, Liang, and Albert, 1988; Heagerty and Zeger, 2000).

The fixed-effect design should support:

  • categorical wave indicators;
  • fully observed predictors;
  • declared predictor-by-wave interactions;
  • additional auxiliary variables associated with the response or missingness process.

A categorical-wave imputation model is the safe default because it contains a linear-time analysis as a restricted special case and better protects A8.

Stable Bernoulli log likelihood

Let

\[ \eta_{ij} = x_{ij}^{\top}\alpha + b_i. \]

The row contribution is

\[ \ell_{ij} = Y_{ij}\eta_{ij} - \log\!\left(1 + e^{\eta_{ij}}\right), \]

implemented stably as:

y * eta - np.logaddexp(0.0, eta)

Observed-data likelihood

The participant contribution is

\[ L_i(\alpha,\tau) = \int \left[ \prod_{j:R_{ij}=1} p_{ij}(b)^{Y_{ij}} \{1-p_{ij}(b)\}^{1-Y_{ij}} \right] \phi(b;0,\tau^2)\,db. \]

LongMI should approximate this one-dimensional integral using the shared Gauss–Hermite quadrature machinery. Adaptive or carefully centered Gauss–Hermite quadrature is a standard numerical strategy for random-effect integrals of this form (Liu and Pierce, 1994).

5. Fit and imputation algorithm

The imputation loop must propagate uncertainty in both the fitted parameters and the missing responses; fixing the parameters and only drawing outcomes generally understates uncertainty (Rubin, 1987; Wang and Robins, 1998).

Fit

  1. Build the fixed-effect design while preserving declared participant and wave order.
  2. Verify that at least one event and one nonevent are observed.
  3. Verify observed-data design rank.
  4. Obtain stable starting values.
  5. Maximize the quadrature-approximated observed-data likelihood.
  6. check optimizer success and the per-observation score norm;
  7. validate curvature and reject a materially singular or indefinite information matrix;
  8. construct the documented large-sample parameter covariance;
  9. return an immutable BernoulliFit containing the model specification, fitted quantities, numerical diagnostics, data fingerprint, and ValidityDeclaration.

Impute

For imputation \(m\):

  1. Draw

$$ (\alpha^{(m)},\log\tau^{(m)}) $$

from the documented large-sample normal approximation.

  1. For a participant with observed outcomes, draw

$$ b_i^{(m)} \sim p!\left( b_i\mid Y_i^{\mathrm{obs}},X_i, \alpha^{(m)},\tau^{(m)} \right) $$

using the adaptive numerical grid.

  1. For a participant with no observed outcomes, draw

$$ b_i^{(m)} \sim N!\left(0,{\tau^{(m)}}^2\right). $$

  1. Form

$$ p_{ij}^{(m)} = \mathrm{expit}!\left( x_{ij}^{\top}\alpha^{(m)} + b_i^{(m)} \right). $$

  1. Draw each missing outcome:

$$ Y_{ij}^{\mathrm{mis},(m)} \sim \mathrm{Bernoulli}!\left(p_{ij}^{(m)}\right). $$

Never round probabilities. The backend must draw Bernoulli outcomes. Replacing a missing outcome with a fitted probability, threshold, or rounded probability is not proper multiple imputation.

6. Separation, identifiability, and numerical safeguards

Binary mixed models add failure modes that deserve explicit handling. In particular, complete and quasi-complete separation can make ordinary logistic maximum-likelihood estimates nonexistent or effectively infinite (Albert and Anderson, 1984).

Hard failures

The backend should refuse to impute when:

  • no outcomes are observed;
  • all observed outcomes are 0;
  • all observed outcomes are 1;
  • the observed fixed-effect design is rank deficient;
  • the optimizer fails materially;
  • the normalized score remains large;
  • the information matrix is singular or materially indefinite;
  • repeated parameter-draw attempts produce nonfinite quantities;
  • the conditional random-effect grid retains material boundary mass.

Diagnostics

BernoulliFit.diagnostics should expose:

  • observed events and nonevents overall and by wave;
  • fixed-effect design rank;
  • optimizer success, message, iteration count, and objective;
  • raw and per-observation gradient norms;
  • curvature eigenvalues and condition number;
  • maximum absolute fixed-effect estimate;
  • quadrature node count and sensitivity;
  • random-effect grid boundary mass and expansion count;
  • rejected parameter draws;
  • accepted \(\tau\) range;
  • observed and imputed prevalence by wave.

A large coefficient or weak curvature should produce an explicit separation warning or failure, never silent near-deterministic imputations.

7. Delta-adjusted MNAR sensitivity

For binary outcomes, delta acts only on the logit scale. This is a controlled pattern-mixture sensitivity parameter: it encodes a specified departure from MAR rather than an estimated correction (Leacy et al., 2017; Cro et al., 2020).

For binary outcomes, delta acts only on the logit scale:

\[ \mathrm{logit}p_{ij}^{\delta} = \mathrm{logit}p_{ij}^{\mathrm{MAR}} + \delta_{ij}. \]

Therefore,

\[ e^{\delta_{ij}} \]

multiplies the conditional odds of the missing response. It does not multiply the probability.

from longmi import DeltaAdjustment

scenario = DeltaAdjustment(
    delta=np.log(0.80),
    scale="linear_predictor",
    where={"exposure": 1},
    times=(6, 9, 12),
    label="group 1 follow-up odds multiplied by 0.80",
)

The package should:

  • apply delta only to missing responses;
  • reject scale="outcome" for binary outcomes;
  • validate referenced columns, values, and waves;
  • reject overlapping rules by default;
  • record the complete scenario and realized per-row delta vector;
  • make clear that delta values are sensitivity assumptions, not estimated corrections.

Why targeted deltas matter. A common delta may change the overall prevalence trajectory without materially changing a group-by-time contrast. Sensitivity of an interaction often requires differential assumptions across groups or waves.

8. Conditional imputer versus marginal logistic GEE

The imputer is subject-specific, whereas the intended GEE analysis is marginal. This distinction is fundamental for longitudinal binary outcomes (Liang and Zeger, 1986; Zeger, Liang, and Albert, 1988; Heagerty and Zeger, 2000).

The imputer is subject-specific:

\[ \mathrm{logit} P(Y_{ij}=1\mid X_{ij},b_i) = X_{ij}^{\top}\alpha+b_i. \]

A logistic GEE is marginal:

\[ \mathrm{logit} P(Y_{ij}=1\mid X_{ij}) = X_{ij}^{\top}\beta. \]

After integrating over the random intercept,

\[ P(Y_{ij}=1\mid X_{ij}) = \int \mathrm{expit}\!\left(X_{ij}^{\top}\alpha+b\right) \phi(b;0,\tau^2)\,db, \]

which is not generally equal to \(\mathrm{expit}(X_{ij}^{\top}\beta)\). Thus,

\[ \alpha \neq \beta \qquad \text{in general}. \]

Compatibility is evaluated, not assumed. The Bernoulli random-intercept imputer followed by a marginal logistic GEE is not automatically congenial. LongMI should describe this workflow as conditionally compatible and evaluate its bias, standard-error calibration, and coverage against the marginal GEE target in simulation.

The relevant package-level qualification is documented in GEE after multiple imputation. MI followed by GEE has been studied directly for incomplete longitudinal binary outcomes, but its validity depends on the imputation and analysis specifications rather than on the workflow label alone (Beunckens, Sotto, and Molenberghs, 2008; Lipsitz, Fitzmaurice, and Weiss, 2020).

9. Validity declaration

A fitted Bernoulli backend should expose a declaration equivalent to:

ValidityDeclaration(
    missingness_assumption="MAR",
    mar_empirically_testable=False,
    sampling_unit="participant",
    parameter_uncertainty_propagated=True,
    outcome_uncertainty_propagated=True,
    observed_outcomes_preserved=True,
    analysis_terms_in_imputation_model=True,
    longitudinal_dependence_modeled=True,
    analysis_nested_in_imputation_model=False,
    congeniality_status=(
        "conditionally supported; marginal GEE workflow "
        "evaluated by simulation"
    ),
    supported_outcome_types=("binary",),
    notes=(
        "Parameter uncertainty uses a large-sample normal approximation; "
        "the participant random effect is drawn from a controlled numerical "
        "approximation to its conditional distribution."
    ),
)

For a delta scenario:

missingness_assumption = "MNAR(delta)"
mnar_sensitivity_performed = True

10. Validation program

Validation should be layered rather than summarized by one test suite. The simulation design and reporting should predefine estimands, data-generating mechanisms, performance measures, Monte Carlo error, and failure classifications (Morris, White, and Crowther, 2019).

Unit and invariant tests

  • binary support and rejection of invalid values;
  • observed-value preservation;
  • completion of every eligible missing response;
  • participant and wave ordering;
  • participants with no observed outcomes;
  • deterministic reproducibility under a fixed seed;
  • separation, rank, optimizer, curvature, and grid failures;
  • rejection of outcome-scale delta;
  • targeted-rule selection, overlap detection, and metadata.

Numerical tests

  • row log likelihood against a trusted scalar calculation;
  • Gauss–Hermite likelihood against high-precision integration on a tiny fixture;
  • gradient agreement with finite differences;
  • convergence as quadrature nodes increase;
  • conditional random-effect density normalization;
  • adaptive-grid boundary-mass control.

Repeated-sampling simulations

The target of the end-to-end workflow is the marginal GEE solution \(\beta^{\star}\), not the conditional mixed-model coefficient \(\alpha\).

Required validity scenarios:

  • MCAR;
  • monotone MAR;
  • intermittent MAR;
  • high missingness;
  • low prevalence;
  • moderate and strong within-participant dependence.

Required expected-failure scenarios:

  • omitted predictor-by-wave interaction;
  • omitted auxiliary predictor related to outcome and response;
  • random-slope data fitted with a random-intercept imputer;
  • MNAR data analyzed under MAR;
  • rare-event or practical-separation stress cases.

Required sensitivity studies:

  • common delta curve;
  • group-specific delta curve;
  • group-by-wave delta surface;
  • monotonicity of the relevant pooled contrast as delta changes.

Report:

\[ \mathrm{Bias},\quad \mathrm{MCSE}(\mathrm{Bias}),\quad \mathrm{RMSE},\quad SD_{\mathrm{emp}},\quad \overline{SE}, \]
\[ \frac{\overline{SE}}{SD_{\mathrm{emp}}},\quad \mathrm{Coverage}_{95\%},\quad \mathrm{MCSE}(\mathrm{Coverage}),\quad \overline{\mathrm{FMI}},\quad \text{failure rate}. \]

Maintain:

  • a small CI smoke suite;
  • a release-scale archived simulation suite;
  • the classification vocabulary validated in evaluated scenarios, expected statistical failure, numerical failure, and not yet evaluated.

Cross-language validation

Use a shared binary fixture and distinguish:

Comparison Standard
Data validation, ordering, masks, design matrices exact
Row and marginal likelihoods at fixed parameters numerical
Quadrature and conditional random-effect density numerical
Rubin pooling exact
Independent stochastic imputations distributional
Bias and coverage statistical

An initial comparison with mice(method="logreg") is a cross-method statistical comparison, not backend parity. It complements, rather than replaces, the longitudinal-model validation described in the recent tutorial literature (Wijesuriya et al., 2025).

11. Public API sketch

from longmi import LongitudinalData, pool_rubin
from longmi.analysis import StatsmodelsGEE
from longmi.impute import BernoulliImputer

long_data = LongitudinalData(
    frame,
    id_col="participant",
    time_col="wave",
    outcome_col="outcome",
    predictor_cols=("exposure", "age", "baseline_outcome"),
    outcome_type="binary",
    times=(0, 3, 6, 9, 12),
)

imputer = BernoulliImputer(
    time_interactions=("exposure",),
    n_quad=25,
)

fit = imputer.fit(long_data)
completed = fit.impute(m=50, random_state=20260724)

analysis = StatsmodelsGEE(
    formula="outcome ~ exposure * wave + age + baseline_outcome",
    groups="participant",
    time="wave",
    family="binomial",
    cov_struct="exchangeable",
)

estimates = completed.analyze(analysis)
pooled = pool_rubin(estimates, validity=completed.declaration)

print(pooled.summary())
print(pooled.validity_report())

12. Package-only pull-request sequence

Stage 0 — shared-core hardening

Finish outstanding ordering, numeric-validation, covariance, GEE-ordering, and stochastic-draw safeguards.

Done when: the continuous and count backends retain their validated behavior and all shared contracts are stable.

Stage 1 — binary data and scenario contracts

Implement:

  • outcome_type="binary";
  • strict binary completion checks;
  • targeted DeltaAdjustment rules;
  • DeltaScenario;
  • deterministic tests for selection, overlap, order, and metadata.

Done when: all binary and scenario invariants are tested without introducing a stochastic backend.

Stage 2 — shared random-intercept infrastructure

Extract reusable quadrature, optimization, curvature, participant aggregation, adaptive-grid, and parameter-draw helpers from the negative-binomial backend.

Done when: existing negative-binomial reference results remain unchanged within prespecified numerical tolerance.

Stage 3 — Bernoulli backend

Implement BernoulliImputer, BernoulliFit, diagnostics, failure handling, and posterior-predictive binary draws.

Done when: unit, numerical, and invariant tests pass.

Stage 4 — end-to-end statistical validation

Add the marginal-GEE target calculation, validity scenarios, expected-failure scenarios, delta studies, and CI smoke suite.

Done when: prespecified smoke gates pass and deliberate assumption violations fail in the expected direction.

Stage 5 — R comparison and public example

Add a public binary example, shared fixtures, R logistic-imputation comparison, and machine-readable outputs.

Done when: the documentation and cross-language classifications are reproducible from a clean environment.

Stage 6 — 0.2.0a1 release validation

Run the release-scale simulation suite, strict documentation build, supported-Python matrix, wheel/sdist audit, clean-install smoke, release bundle, and Trusted Publishing workflow.

13. Release gates

Gate Requirement
Binary contract invalid observed or imputed values are impossible to pass silently
Observed-value preservation mechanically guaranteed
Parameter uncertainty propagated and documented as approximate
Outcome uncertainty propagated through Bernoulli draws
Numerical stability convergence, curvature, grid, and draw diagnostics exposed
Delta sensitivity logit-scale, targetable, auditable, and missing-outcome-only
A1–A8 mapping explicit in documentation and ValidityDeclaration
Conditional/marginal distinction documented and directly evaluated
MI-to-GEE simulation bias, SE calibration, coverage, FMI, and failure rates archived
Expected failures visible under model misspecification and MNAR
R comparison complete and correctly labeled
Public example open, reproducible, and study-neutral
CI and packaging green across supported environments
Nonclaims explicit in release notes and documentation

14. Documentation additions

The feature should add or update:

docs/algorithms/bernoulli_glmm.md
docs/explanation/binary_outcomes.md
docs/explanation/conditional_vs_marginal_logistic.md
docs/how_to/impute_binary_outcomes.md
docs/how_to/run_binary_delta_sensitivity.md
docs/how_to/choose_an_imputer.md
docs/theory/assumptions.md
docs/theory/gee_after_imputation.md
docs/project-status.md

Each page should link back to the canonical A1–A8 assumptions rather than restating them inconsistently.

15. Package nonclaims

LongMI 0.2.0a1 should not claim that it:

  • verifies MAR;
  • identifies the correct MNAR delta;
  • recovers information that was never observed;
  • guarantees compatibility between arbitrary imputation and analysis models;
  • equates subject-specific and marginal logistic coefficients;
  • solves structural eligibility or competing-event problems;
  • supports arbitrary separation or rare-event settings;
  • establishes validity outside the evaluated simulation scenarios.

The defensible release statement is:

LongMI implements a random-intercept Bernoulli imputation backend with explicit validity declarations, numerical diagnostics, logit-scale MNAR sensitivity, and repeated-sampling evaluation of selected marginal logistic GEE workflows under prespecified data-generating scenarios.

LongMI documentation

Methodological references cited in this plan

  1. Albert A, Anderson JA. On the existence of maximum likelihood estimates in logistic regression models. Biometrika. 1984;71(1):1–10. doi:10.1093/biomet/71.1.1.
  2. Barnard J, Rubin DB. Small-sample degrees of freedom with multiple imputation. Biometrika. 1999;86(4):948–955. doi:10.1093/biomet/86.4.948.
  3. Bartlett JW, Hughes RA. Bootstrap inference for multiple imputation under uncongeniality and misspecification. Stat Methods Med Res. 2020;29(12):3533–3546. doi:10.1177/0962280220932189.
  4. Beunckens C, Sotto C, Molenberghs G. A simulation study comparing weighted estimating equations with multiple imputation based estimating equations for longitudinal binary data. Comput Stat Data Anal. 2008;52(3):1533–1548. doi:10.1016/j.csda.2007.04.020.
  5. Cro S, Morris TP, Kenward MG, Carpenter JR. Sensitivity analysis for clinical trials with missing continuous outcome data using controlled multiple imputation: a practical guide. Stat Med. 2020;39(21):2815–2842. doi:10.1002/sim.8569.
  6. Heagerty PJ, Zeger SL. Marginalized multilevel models and likelihood inference. Stat Sci. 2000;15(1):1–19. doi:10.1214/ss/1009212671.
  7. Leacy FP, Floyd S, Yates TA, White IR. Analyses of sensitivity to the missing-at-random assumption using multiple imputation with delta adjustment. Am J Epidemiol. 2017;185(4):304–315. doi:10.1093/aje/kww107.
  8. Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22. doi:10.1093/biomet/73.1.13.
  9. Lipsitz SR, Fitzmaurice GM, Weiss RD. Using multiple imputation with GEE with non-monotone missing longitudinal binary outcomes. Psychometrika. 2020;85(4):890–904. doi:10.1007/s11336-020-09729-y.
  10. Little RJA, Rubin DB. Statistical Analysis with Missing Data. 3rd ed. Wiley; 2019. doi:10.1002/9781119482260.
  11. Liu Q, Pierce DA. A note on Gauss–Hermite quadrature. Biometrika. 1994;81(3):624–629. doi:10.1093/biomet/81.3.624.
  12. Meng XL. Multiple-imputation inferences with uncongenial sources of input. Stat Sci. 1994;9(4):538–573. doi:10.1214/ss/1177010269.
  13. Morris TP, White IR, Crowther MJ. Using simulation studies to evaluate statistical methods. Stat Med. 2019;38(11):2074–2102. doi:10.1002/sim.8086.
  14. Robins JM, Wang N. Inference for imputation estimators. Biometrika. 2000;87(1):113–124. doi:10.1093/biomet/87.1.113.
  15. Rubin DB. Inference and missing data. Biometrika. 1976;63(3):581–592. doi:10.1093/biomet/63.3.581.
  16. Rubin DB. Multiple Imputation for Nonresponse in Surveys. Wiley; 1987. doi:10.1002/9780470316696.
  17. Wang N, Robins JM. Large-sample theory for parametric multiple imputation procedures. Biometrika. 1998;85(4):935–948. doi:10.1093/biomet/85.4.935.
  18. White IR, Royston P, Wood AM. Multiple imputation using chained equations: issues and guidance for practice. Stat Med. 2011;30(4):377–399. doi:10.1002/sim.4067.
  19. Wijesuriya R, Moreno-Betancur M, Carlin JB, White IR, Quartagno M, Lee KJ. Multiple imputation for longitudinal data: a tutorial. Stat Med. 2025;44(3–4):e10274. doi:10.1002/sim.10274.
  20. Yang S, Kim JK. A note on multiple imputation for method of moments estimation. Biometrika. 2016;103(1):244–251. doi:10.1093/biomet/asv073.
  21. Zeger SL, Liang KY, Albert PS. Models for longitudinal data: a generalized estimating equation approach. Biometrics. 1988;44(4):1049–1060. doi:10.2307/2531734.