Skip to content

Handle missing outcomes via the delta missingness mechanism (#32) - #34

Open
ck37 wants to merge 2 commits into
masterfrom
claude/nifty-albattani-qae38s
Open

Handle missing outcomes via the delta missingness mechanism (#32)#34
ck37 wants to merge 2 commits into
masterfrom
claude/nifty-albattani-qae38s

Conversation

@ck37

@ck37 ck37 commented Sep 14, 2026

Copy link
Copy Markdown
Owner

Fixes #32.

Confirming the bug

Running the simulation from #32 verbatim on master, every CV-TMLE fold dies with:

Error in if (max(Y_star) > 1 | min(Y_star) < 0) { : missing value where TRUE/FALSE needed

and the run ends with "No VIMs could be calculated."

Two details matter for the fix:

  • The block Alan flagged isn't the one that runs. The simulation is all-numeric, so it goes through vim-numerics.R, which had its own copy of the threshold at line 298. estimate_tmle2.R had a third copy at lines 41-52.
  • Removing that block is necessary but not sufficient. The < 10 rule only ever subset the training fold. Yv/deltav were never touched, so apply_tmle_to_validation() always received NA outcomes — it had a bare # TODO: handle delta in some way? where the missingness handling belonged. That is why the error persists above 10 missing, and why it also fails at 5 missing (verified: 5 and 11 both fail identically on master).

The fix

Missing observations are now retained everywhere and handled by the missingness mechanism, which is what the TODO in vim-factors.R intended:

  • estimate_tmle2() no longer drops observations below the threshold of 10, and returns g_delta_model so the fit can be applied to a validation fold. Its CV fold count for the missingness model is reduced (and stratification skipped) when the delta indicator is sparse, which otherwise made SuperLearner error out and silently fall back to glm.
  • apply_tmle_to_validation() predicts P(Delta = 1 | A = 1, W) on the validation fold and builds the clever covariate as A * delta / (g1W * g.Delta). Observations missing Y get zero weight instead of poisoning the bounds check and every later calculation.
  • estimate_pooled_results() fits the fluctuation on the observed outcomes and uses that clever covariate in the influence curve, so missing outcomes no longer turn the standard errors into NAs.
  • vim_factors() / vim_numerics() drop their copies of the threshold. Their binary-outcome and minYs checks now ignore NAs — previously length(unique(Yt)) == 2 returned 3 with any missing Y, silently routing a binary outcome down the continuous branch.
  • create_cv_folds() keeps stratifying a binary outcome that has NAs, and spreads the missing outcomes evenly across folds.

Separate defect fixed in the same code path

estimate_pooled_results() used stats::offset(logit_Q_hat) in its fluctuation formula. Namespace-qualifying offset() stops R's terms machinery from treating it as an offset, so logit_Q_hat was fit as an ordinary covariate and epsilon came back with two elements, which then recycled across observations in epsilon * data$H1W. Introduced in 386f6a8; the recycling warning only surfaces when the row count is odd, so it went unnoticed.

On a fixed synthetic fixture (tests/testthat/test-pooled-fluctuation.R), the difference is unambiguous:

stats::offset() (before) offset() (after)
length(epsilon) 2 1
max abs influence-curve mean 0.345 0.069
fold estimate vs mean(Q_hat) off by 0.20–0.22 within 0.065

End-to-end on null data (Y independent of every X), the improvement is real but distributional, not per-seed — over 10 seeds:

n mean abs estimate sd 95% CI coverage
before 27 0.162 0.202 21/27 (78%)
after 30 0.102 0.119 27/30 (90%)

Two honest caveats on that table: individual seed batches vary a lot (an earlier 4-seed batch happened to give 7/12 vs 12/12, which overstates the separation), and 90% is still short of the nominal 95% — this PR improves the fluctuation but does not make the intervals exact.

Results on complete data change with this PR because of that fix, so it is worth a look before merging.

Verification

  • The simulation from handling situations where the number of observations missing either A or Y is >0 and < 10. #32: 0 errors, full results table.
  • Also verified with factor variables, a continuous (gaussian) outcome, and mixed factor/numeric data with missing values in both Y and X.
  • Estimates on data with no missing values at all are bit-identical to master with the missing-data changes alone; the offset fix then moves them.
  • R CMD check passes (only pre-existing NOTEs) and the existing test suite passes.
  • New tests/testthat/test-missing-outcome.R: 6 tests covering the issue's simulation, the 9-vs-11 missing boundary that the old threshold straddled, estimate_tmle2() with and without missingness, the zero-weighting in apply_tmle_to_validation(), and stratified fold creation with NA outcomes.
  • New tests/testthat/test-pooled-fluctuation.R: 3 tests / 8 assertions on the offset defect, deterministic and sub-second. Verified to fail on master (5 failures) and pass here.

Note on #33

#33 removes only the vim-factors.R block and adds several scratch test scripts at the repository root. It does not touch the numeric path or the validation fold, so it does not resolve the issue. Closing in favor of this PR.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Xo3yCjZifUfAfLdHL1c1Hu

varimpact() could not estimate any VIMs when Y (or A) contained missing
values: every CV-TMLE fold died in apply_tmle_to_validation() with
"missing value where TRUE/FALSE needed", and the run ended with "No VIMs
could be calculated". The arbitrary "drop observations if fewer than 10
are missing" rule that issue #32 points at only ever applied to the
training fold, so the validation fold always received NA outcomes -
which is why the failure happened with 5 missing outcomes as well as
with 11.

Missing observations are now retained everywhere and handled by the
missingness mechanism, as the TODO in vim-factors.R intended:

- estimate_tmle2() no longer drops observations below the threshold of
  10, and returns g_delta_model so the fit can be applied to a
  validation fold. Its CV fold count for the missingness model is
  reduced (and stratification skipped) when the delta indicator is
  sparse, which otherwise made SuperLearner fail and silently fall back
  to glm.
- apply_tmle_to_validation() predicts P(Delta = 1 | A = 1, W) on the
  validation fold and builds the clever covariate as
  A * delta / (g1W * g.Delta). Observations missing Y get zero weight
  instead of poisoning the bounds check and every later calculation.
- estimate_pooled_results() fits the fluctuation on the observed
  outcomes and uses that clever covariate in the influence curve, so
  missing outcomes no longer turn the standard errors into NAs.
- vim_factors()/vim_numerics() drop their copies of the threshold, and
  their binary-outcome and minYs checks now ignore NAs rather than
  misclassifying a binary outcome with missingness as continuous.
- create_cv_folds() keeps stratifying a binary outcome that has NAs, and
  spreads the missing outcomes evenly across folds.

Also fixes a separate defect in the same code path: the pooled
fluctuation used stats::offset(logit_Q_hat) in its glm formula.
Namespace-qualifying offset() stops the terms machinery from treating it
as an offset, so logit_Q_hat was fit as an ordinary covariate and
epsilon came back with two elements, which then recycled across
observations in "epsilon * data$H1W". Under a pure null the estimates
had sd 0.21 and ranged over [-0.35, 0.59]; with offset() called
unqualified they have sd 0.08 and range over [-0.17, 0.12]. Results on
complete data therefore change with this commit.

Verified against the simulation in the issue (0 errors, full results
table), and on factor, continuous-outcome and mixed data with missing
outcomes. Estimates on data with no missing values at all are
bit-identical to before, apart from the offset fix. R CMD check and the
existing test suite pass.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Xo3yCjZifUfAfLdHL1c1Hu
The stats::offset() defect this branch fixes had no test. These cover it
directly, using a fixed synthetic pair of validation folds rather than a
varimpact() run, so they are deterministic and take under a second:

- the fluctuation estimates exactly one coefficient, named HAW. Two
  elements means offset() was not honored and epsilon recycles across
  observations downstream.
- per-fold influence curves are finite and centered near zero, within
  the 1/sqrt(n) Monte Carlo error expected of a correctly targeted
  influence curve.
- targeting does not move the estimate far from the initial Q, which is
  well specified in this fixture.

Verified to fail on master (5 failures across the 3 tests: epsilon has
length 2, influence-curve means reach 0.35, and the fold estimates land
0.20 away from mean(Q_hat)) and to pass on this branch.

The fold row counts are odd on purpose. An even total made the recycled
epsilon a silent multiple of the row count, so only odd counts produced
R's "longer object length" warning - which is why this went unnoticed.

Note on what is not tested here: an end-to-end check on null data
separates the two code paths only in distribution, not per seed. Over 10
seeds, master gives 78% CI coverage with sd 0.202 and this branch gives
90% with sd 0.119 - a real improvement, but not a per-seed guarantee, so
asserting it would be flaky. Worth recording that 90% is still short of
the nominal 95%: this branch improves the fluctuation but does not make
the intervals exact.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Xo3yCjZifUfAfLdHL1c1Hu
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

handling situations where the number of observations missing either A or Y is >0 and < 10.

2 participants