Fix Issue #8: Transform continuous outcome estimates back to original scale

[ck37]

🎯 Overview

Addresses Issue #8: variable importance estimates for a continuous outcome were reported on the internal [0, 1] scale rather than the scale of Y itself.

🔧 Problem

varimpact() maps a continuous outcome into [0, 1] with Y_star = (Y - Qbounds[1]) / diff(Qbounds) before running the CV-TMLE, and never mapped the results back. If the outcome ranged 10–50, estimates came back between 0 and 1.

A second, independent problem sat on top of it. estimate_tmle2() applied plogis() to Qstar twice when mapping the targeted predictions back to the outcome scale — a duplicated line; tmle::tmle() has it once. After the first line Qstar is already on the outcome's scale, so the second plogis() distorts it. This is pre-existing on master, not caused by the rescaling.

How bad that was depends on the outcome's range

This corrects an earlier version of this description, which said gaussian "produced no results whatsoever". That is true for some outcomes and not others, because plogis() only saturates once its argument is far from zero. Measured against master:

outcome	master
ranging 31 to 68	results_all is NULL	fully saturated: every prediction collapses to the upper bound, the training estimate is identical in every bin, which.max() and which.min() pick the same bin, every fold is discarded as "min and max level are the same"
ranging −2 to 3	returns results	doesn't collapse — comes back distorted, roughly the original-scale estimate divided by the outcome's range

The second case is the harder one to notice: nothing errors, you just get plausible-looking numbers on the wrong scale.

For a binary outcome the bounds are c(0, 1) and the extra transform was plogis() of a probability, which shifted theta into (0.5, 0.731) but left bin selection intact.

💡 Solution

R/estimate_tmle2.R — removed the duplicated plogis() line, matching tmle::tmle().

R/estimate_pooled_results.R — gained a Qbounds argument (default c(0, 1)) and applies the inverse of the Y -> Y_star map to the fluctuated Q_star and to Y_star, immediately after the fluctuation. One transformation site: the per-fold thetas, the risk difference, the risk ratio and the influence curves are all built from those two, so they inherit the correct scale.

theta_original = theta_star * diff(Qbounds) + Qbounds[1]
IC_original    = diff(Qbounds) * IC_star

The location shift cancels out of the influence curve because both of its terms are differences: HAW * (Y_star - Q_star) and Q_star - mean(Q_star).

R/vim-numerics.R, R/vim-factors.R — pass down the Qbounds that varimpact() computed.

No map_to_ystar flag

An earlier revision gated the transform behind a family/uniqueness heuristic duplicated in three places. It isn't needed: Qbounds is c(0, 1) for a binary outcome, which makes the transform the identity. The flag was also read in vim-numerics.R at the per-bin call but only assigned ~65 lines later, and it isn't a formal of vim_numerics() — so that was an unbound-variable error, swallowed by the surrounding tryCatch(..., error = ...), silently dropping every bin's pooled result.

🔄 Brought onto current master

Opened against 938867d, which master has long since left behind — #40, #41, #43, #46, #31, #51 and #52 have all landed, and all of them bear on the files this PR touches.

The conflicts had a cause worth recording. An earlier "merge master" commit on this branch, d496f70, was not actually a merge commit — it has a single parent. The merge had been staged with --no-commit, but the verification that followed checked out master to compare results, which cleared MERGE_HEAD; the commit then captured master's content without its history. So the merge base stayed pinned at 938867d and every file master had added since looked like an add/add conflict.

Four files conflicted for that reason alone — R/tmle_estimate_g.R and three test files — and none is touched by this PR. Confirmed by diffing this PR's own commit against its base (938867d..a1bdd49): six files changed, none of those four. Master's version was taken for all four, and each now matches master byte for byte. The current merge commit has two parents, so this won't recur.

R/estimate_tmle2.R auto-merged and carries every side correctly: this PR's de-duplicated plogis(), #46's min_cell_size guard, and #52's removal of the pDelta1/g.Deltaform locals.

🧪 Testing

tests/testthat/test-continuous-outcome-scale.R checks the rescaling algebraically against synthetic fold results (thetas shift and scale by Qbounds, influence curves scale only, epsilon untouched), runs a gaussian model end to end, and keeps a binomial run as a regression check.

Verified on the merged tree, against current master

check	result
Full suite	252 passing, 0 failures
pkgdown::build_site()	exits 0, no problems
Binary results vs current master	identical() on the same seeded data
Wide-range gaussian (Y 31–68)	9.45 / 7.06 / 3.34, correctly ordered, original scale — master returns nothing
Project coverage	83.00% → 84.51%

⚠️ About the warning count

The suite reports 16 warnings against master's 7. All 9 new ones are NaNs produced, from the three gaussian runs in test-varimpact.R, whose outcome is Y + rnorm() on a binary Y and therefore straddles zero.

That is the limitation recorded below: the risk ratio takes log(EY1 / EY0), undefined when a mean is negative. They appear because gaussian runs produce results at all — on master those runs yield no risk ratio to take a log of. Not a regression, and the risk difference columns are unaffected.

⚠️ Known limitation (not addressed here)

The risk ratio EY1 / EY0 is only interpretable for a strictly positive outcome. Now that continuous estimates are on the original scale, an outcome whose range straddles zero can give a negative or undefined RR, since compile_results() takes log(psi_rr). Previously both means were in [0, 1], so the log was finite but meaningless. The risk difference columns are unaffected, and compile_results() filters on thetaV rather than thetaV_rr, so RR NaNs do not drop variables. Recorded under "Known limitations" in NEWS.md.

📝 Files changed

Against current master: NEWS.md, R/estimate_pooled_results.R, R/estimate_tmle2.R, R/vim-factors.R, R/vim-numerics.R, tests/testthat/test-continuous-outcome-scale.R.

📚 Related

• Fixes Transform continuous outcome estimates #8

[openhands-ai]

Looks like there are a few issues preventing this PR from being merged!

• GitHub Actions are failing:
  • R-CMD-check

If you'd like me to help, just leave a comment, like

@OpenHands please fix the failing actions on PR #30

Feel free to include any additional details that might help me get this PR into a better state.

_{^{You can manage your notification settings}}

[ck37]

Pushed a1bdd49: merged current master in and reworked this. The "Known Issues" section — vim$results_all coming back NULL — had a specific cause, and it wasn't the rescaling.

estimate_tmle2() applies plogis() to Qstar twice.

if (map_to_ystar) {
  Qstar <- plogis(Qstar)*diff(stage1$ab)+stage1$ab[1]
  Qstar <- plogis(Qstar)*diff(stage1$ab)+stage1$ab[1]   # <- duplicated

tmle::tmle() has this line once (checked against the installed source). After the first line Qstar is already on the outcome's scale, so the second plogis() saturates it. For a continuous outcome with a wide range every value collapses to stage1$ab[2], which makes the training theta identical in every bin — which.max() and which.min() pick the same bin, every fold is discarded as "min and max level are the same", and varimpact() returns nothing for family = "gaussian". So the rescaling could never be exercised end to end. Confirmed pre-existing on master, not introduced by this branch.

With that line removed, a gaussian run on N = 300 gives:

      Type Estimate          CI95      P-value  Est. RR
X1 Ordered 6.359346 (5.39 - 7.33) 0.000000e+00 1.137994
X2 Ordered 5.095404 (4.12 - 6.07) 0.000000e+00 1.109116
X3 Ordered 3.285424 (2.33 - 4.24) 9.118151e-12 1.066816

against Y = 50 + 15*(0.3*X1 + 0.2*X2 - 0.1*X3) + N(0, 3) — original scale, right ordering, and close to the 7.2 / 4.8 / 2.4 you'd expect from a median split.

The rescaling itself is now one transform, not three. estimate_pooled_results() takes a Qbounds argument and inverts varimpact()'s Y -> Y_star map on the fluctuated Q_star (and on Y_star) right after the fluctuation. The per-fold thetas, the risk difference, the risk ratio and the influence curves are all built from those two, so they inherit the correct scale:

theta_original = theta_star * diff(Qbounds) + Qbounds[1]
IC_original    = diff(Qbounds) * IC_star

The location shift cancels out of the influence curve because both of its terms are differences (HAW * (Y_star - Q_star) and Q_star - mean(Q_star)).

Dropped the map_to_ystar flag. Two reasons:

• Qbounds is c(0, 1) for a binary outcome, so the transform is already the identity. No family check is needed, and the heuristic didn't need to be duplicated in three places.
• In vim-numerics.R the flag was read at the per-bin estimate_pooled_results() call (~line 727) but only assigned ~65 lines later. It isn't a formal of vim_numerics(), so that was an unbound-variable error — swallowed by the surrounding tryCatch(..., error = function(error) ...), which silently dropped every bin's pooled result. That was a second, independent reason the tests couldn't pass.

Binary is unchanged. Verified a binomial run on this branch reproduces origin/master's estimates exactly (0.15814193, 0.14445406, 0.07116556).

Tests rewritten. The old file asserted that fewer than 90% of estimates fell in [0, 1], which passes vacuously when results_all is NULL. The new one:

• checks the rescaling algebraically against synthetic fold results — thetas shift and scale by Qbounds, influence curves scale only, epsilon is untouched;
• runs a gaussian model end to end and checks the estimates land near the known contrasts on the original scale;
• keeps a binomial regression check.

Full suite passes.

One limitation worth recording, not fixed here: the risk ratio EY1 / EY0 is only interpretable for a strictly positive outcome. Now that continuous estimates are on the original scale, an outcome whose range straddles zero can give a negative or undefined RR (compile_results() takes log(psi_rr)). Previously both means were in [0, 1] so the log was finite but meaningless. The risk difference columns are unaffected, and compile_results() filters on thetaV, not thetaV_rr, so RR NaNs don't drop variables. Noted under "Known limitations" in NEWS.md.

This is still marked draft — happy to take it out of draft if the above looks right to you.

---

Generated by Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 84.51%. Comparing base (f931e42) to head (4715011).
⚠️ Report is 14 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #30      +/-   ##
==========================================
+ Coverage   83.00%   84.51%   +1.50%     
==========================================
  Files          30       30              
  Lines        2366     2331      -35     
==========================================
+ Hits         1964     1970       +6     
+ Misses        402      361      -41     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.