What a Double-Robustness Simulation Can Establish
A simulation of longitudinal targeted maximum likelihood estimation (TMLE) should do more than compare a single estimate with a known answer. Its central purpose is to test what happens when the nuisance models are correct in different combinations, while also recording whether estimates and standard errors remain numerically usable. These are separate questions: robustness is an asymptotic property under assumptions, whereas finite-sample behavior depends on sample size, positivity, model fitting, and the targeting procedure.
Consider a fixed treatment regime \(\bar a=(a_0,\ldots,a_{K-1})\), with covariate history measured before each treatment and a final outcome \(Y\). Let \(\psi_0=E[Y^{\bar a}]\) denote the mean outcome under the regime, identified under consistency, sequential exchangeability, and positivity. In a simulation, the data-generating mechanism (DGM) specifies the joint distribution of baseline covariates, time-varying covariates, treatment, and outcome. The regime-specific mean can then be calculated from that DGM, independently of the estimator being evaluated.
For a fixed regime, longitudinal TMLE is commonly described as doubly robust with respect to the sequential outcome regressions and treatment (and, when relevant, censoring) mechanisms: under the usual identification, positivity, regularity, and targeting conditions, consistency can hold if the complete sequential outcome-regression sequence is consistently estimated, or if the treatment/censoring mechanism sequence is consistently estimated. This is not a guarantee that every arbitrary combination of local model errors works, nor does it guarantee good finite-sample precision. The simulation should therefore define exactly which nuisance components are correct in each scenario. Product-rate conditions and cross-fitting address additional requirements for asymptotic linearity; they do not make positivity or finite-sample instability disappear.
Specify the DGM and Its Target Before Fitting Models
The following two-time-point DGM gives a fully calculable example. Let \(X\) be baseline, \(A_0\) and \(A_1\) binary treatments, \(L_1\) and \(L_2\) binary covariates measured after the respective preceding treatments, and \(Y\) a binary final outcome. Generate \(X\sim\mathrm{Bernoulli}(0.5)\). Treatment assignment is confounded: for example, use \(P(A_0=1\mid X)=\operatorname{expit}(-0.8+1.2X)\) and \(P(A_1=1\mid X,L_1,A_0)=\operatorname{expit}(-1+0.8X+1.1L_1+0.5A_0)\). Generate the covariates and outcome according to
All displayed probabilities lie between zero and one for binary inputs. For the regime \(\bar a=(1,1)\), the target can be derived without simulating treatment assignment. Under the intervention, \(E[X]=0.5\), so
Substituting these expectations into the outcome equation gives
This derivation is valid because the conditional outcome mean is linear in \(X,L_1,L_2,A_0,A_1\). It does not assume that the observed treatment groups resemble the intervention groups; the intervention sets both treatments to one, and the covariate expectations are propagated forward under that setting.
Worked Example: Building the Four Nuisance Scenarios
Suppose the simulation uses the DGM above, targets \(\psi_0=0.6015\), and compares four prespecified model scenarios. In the “correct” models, the fitted conditional means or probabilities contain the variables and functional forms required by the DGM. In deliberately misspecified models, omit specified predictors rather than relying on an accidental fitting failure.
Both sequences correct. Fit the treatment probabilities with the stated logistic linear predictors and fit the covariate and outcome regressions with their stated conditional-mean forms. This is the reference scenario for checking implementation and sampling variability.
Outcome sequence misspecified, treatment sequence correct. Keep both treatment models correctly specified. For example, omit \(L_2\) from the final outcome regression and omit \(L_1\) from the regression for \(L_2\). These omissions deliberately make the sequential outcome regressions inadequate. This scenario probes the treatment-mechanism robustness branch, subject to the estimator’s stated assumptions.
Outcome sequence correct, treatment sequence misspecified. Keep the sequential outcome regressions correct, but fit the second treatment probability without \(L_1\), which is a predictor in its DGM. This probes the outcome-regression branch.
Both sequences misspecified. Use both sets of omissions. Double robustness does not predict protection when both required sequences fail; this is a negative-control scenario for the robustness claim, not a scenario where inconsistency is guaranteed to be large in every finite sample.
Before interpreting results, verify that the stated omissions really make the model wrong under the chosen DGM, that the correct models are estimable, and that the same regime, sample size, fitting rules, and interval construction are used across scenarios. In particular, do not call an estimator robust merely because one finite simulation happens to yield a small bias.
Performance Measures and Their Derivations
Let \(\hat\psi_r\) be the estimate from replicate \(r\), and let \(\widehat{\mathrm{SE}}_r\) be its estimated standard error. The simulation summaries below use all \(R\) replicates unless a fit fails; how failures are handled must be declared in advance.
The variance formula is the ordinary sample variance: its deviations are measured from the replicate mean, not from the known truth. Root mean squared error can also be reported; with the denominator \(R\), its squared value decomposes exactly into squared empirical bias plus the replicate variance with denominator \(R\):
This follows by writing \(\hat\psi_r-\psi_0=(\hat\psi_r-\bar{\hat\psi})+(\bar{\hat\psi}-\psi_0)\), squaring, and summing: the cross-term is zero because \(\sum_r(\hat\psi_r-\bar{\hat\psi})=0\). Coverage is the fraction of intervals that include the truth:
Coverage should be interpreted alongside bias and interval width. An interval can cover poorly because it is centered incorrectly, because its standard error is too small, or both. Conversely, very wide intervals may cover well but be uninformative.
Worked Example: Calculating Bias and Variance
Suppose five illustrative replicate estimates are \(0.58,0.60,0.62,0.59,0.61\), with truth \(0.6015\). These values are a hand calculation for the summaries, not claimed output from a simulation run.
Their sum is \(3.00\), so \(\bar{\hat\psi}=3.00/5=0.60\). The empirical bias is \(0.60-0.6015=-0.0015\).
The deviations from the replicate mean are \(-0.02,0,0.02,-0.01,0.01\); their squared values sum to \(0.001\). Therefore \(\widehat{\operatorname{Var}}=0.001/(5-1)=0.00025\), and the empirical standard deviation is \(\sqrt{0.00025}\approx0.0158\).
If the five estimated standard errors are \(0.018,0.017,0.020,0.016,0.019\), the Wald interval for the estimate \(0.62\), for example, is \(0.62\pm1.96(0.020)=[0.5808,0.6592]\), which includes \(0.6015\). Applying the same check to all five intervals gives five covers out of five. That observed coverage of \(1.00\) is highly imprecise with only five replicates; it is not evidence that the long-run coverage is exactly one.
Compare Performance Without Confounding the Comparison
A useful simulation varies sample size as well as nuisance-model scenario. Increasing \(n\) helps distinguish persistent bias from ordinary sampling noise, while repeated independent replicates reveal variability. Hold the DGM and target fixed across scenarios. If comparing multiple regimes, retain the same simulated data within a replicate where appropriate, so differences are not needlessly obscured by unrelated random draws.
For each scenario and sample size, report at least the number of attempted replicates, the number of successful fits, bias, empirical standard deviation, mean estimated standard error, coverage, and interval width. Monte Carlo uncertainty matters: the standard error of the estimated bias is approximately the replicate standard deviation divided by \(\sqrt{R}\). If the estimated coverage is \(\hat c\), its binomial Monte Carlo standard error is approximately \(\sqrt{\hat c(1-\hat c)/R}\). Thus, a coverage difference of a few percentage points may be simulation noise when \(R\) is modest.
Worked Example: Separating Accuracy From Stability
Imagine that, at a particular sample size, the correct-outcome/wrong-treatment scenario has bias near zero among successful fits but a 12% failure rate caused by extreme treatment probabilities and unstable targeting. The correct-treatment/wrong-outcome scenario has a small failure rate and modest bias. Reporting only the first scenario’s successful-fit bias would conceal a substantial practical problem.
Report the failure rate as \(0.12\) and state precisely how estimates from failed replicates were handled. Do not silently remove failures and then describe the remaining estimates as performance over all generated data sets. Conditional summaries among successful fits may be useful, but label them as conditional. Also report diagnostics such as extreme estimated probabilities, non-finite estimates, convergence warnings, and unusually large estimated standard errors. These expose finite-sample instability even when the successful estimates appear unbiased.
The main alternatives answer related but different questions. Longitudinal inverse-probability weighting emphasizes the treatment/censoring mechanism and can become highly variable with extreme weights. Untargeted sequential regression emphasizes outcome regressions and can retain bias when those regressions are misspecified. Longitudinal TMLE combines sequential regression with targeting, but does not eliminate positivity limitations or guarantee better finite-sample performance in every scenario. A fair comparison applies all methods to the same simulated data sets, uses the same target and failure accounting, and reports uncertainty rather than ranking methods from one summary alone.
Interpretation, Diagnostics, and Limits
The expected pattern is a useful diagnostic, not a pass/fail theorem for a finite run: both-correct estimates should approach the truth; either one-sequence-correct scenario may approach it under the applicable double-robust conditions; and both-wrong scenarios may show persistent bias. At small sample sizes, a nominally correct model can still fit poorly, and the empirical pattern can depart substantially from this expectation.
- Check implementation separately from robustness. In the both-correct case, compare estimated nuisance functions and targeted estimates with known DGM quantities where possible. An indexing or regime-coding error can affect every scenario.
- Inspect positivity and targeting. Summarize estimated treatment probabilities along the regime and record boundary predictions, separation, extreme clever covariates, and failed or unstable updates. A correct model does not repair weak positivity.
- Preserve the full replicate record. Store failures and warnings as outcomes of the procedure. Do not replace failed estimates with convenient values or omit them without disclosure.
- Separate Monte Carlo error from estimator error. Increase \(R\) when performance summaries themselves are noisy; increase \(n\) to study the estimator’s sampling behavior. These changes solve different problems.
Finally, a simulation supports conclusions only for its specified DGM, nuisance learners, sample sizes, and implementation. Correctly specified parametric models in a controlled DGM do not establish robustness to every form of misspecification encountered with machine learning. Likewise, poor performance in a severe positivity scenario does not by itself refute double robustness, whose usual guarantees require positivity and other regularity conditions.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Under the displayed DGM, why can the target for regime \((1,1)\) be computed without using the treatment-assignment probabilities?
- If both the sequential outcome regressions and treatment models are misspecified, does double robustness predict consistency?
- In a simulation with \(R=400\) replicates and estimated coverage \(0.90\), what is the approximate Monte Carlo standard error of coverage?
- Why is it misleading to report bias only among successful fits when some replicates fail?
Solutions. (1) The target is defined under an intervention setting \(A_0=A_1=1\); propagate the covariate distributions under that intervention and then average the outcome regression. Treatment-assignment probabilities describe observed assignment, not the intervention distribution. (2) No. The standard double-robust claim protects when one required nuisance sequence is consistently estimated, under its assumptions; it does not generally protect when both fail. (3) \(\sqrt{0.90(0.10)/400}=\sqrt{0.000225}=0.015\), or about 1.5 percentage points. (4) Failed replicates are part of the procedure’s finite-sample behavior. Excluding them can hide instability and changes the summary to performance conditional on successful fitting.