Response Surfaces and the Bias They Create

Python · PyMC · pymc_bart  ·  LaLonde NSW/CPS, 614 observations

The Disciplined Version Has the Same Disease

The first example left this group in an awkward position. A coherent Bayesian analysis cannot use the propensity score — the likelihood factorises and the assignment model is ancillary — so robustness has to come from the outcome model. But the obvious way to make an outcome model robust failed: a flexible learner attenuated a true effect of 2.0 down to 0.48, because the prognostic signal — everything that predicts earnings regardless of treatment — dwarfed the treatment effect and the trees spent their splits elsewhere.

BART (Chipman, George & McCulloch 2010) is the disciplined version of that idea, and Hill (2011) brought it to causal inference. Each tree is regularized to be a weak learner contributing a small share of the fit, so the ensemble is flexible without any single tree dominating. The site already has a BART example, so the machinery is not new here. What is new is that the disciplined version turns out to have the same disease, and that the cure is a change of parameterisation rather than of tuning.

One caveat governs every number on this page. pymc_bart’s sampler is not reproducible at a fixed seed — the same configuration returns different answers on re-runs — so each figure is one draw from a distribution, and the claims below are stated as ranges wherever a range is what the evidence supports.

The standard recipe puts treatment in the design matrix alongside the eight covariates, and takes counterfactuals afterwards by evaluating the same fitted trees with that column flipped. No propensity score is estimated anywhere. It returns USD 283 against a randomized benchmark of USD 1,794 — 84% low, and the worst adjusted estimate in the observational group, beaten only by the naive comparison that gets the sign wrong.

Is That Just My Tuning?

A result that bad, on a method with a good reputation, is more likely to be the analyst’s fault than the method’s. So the same fit is repeated with four times the trees, twice the draws, and three different seeds — and the check found something I was not looking for.

treesdrawsseedATT95% interval
505000716[−1,137, 2,753]
505001154[−1,588, 1,647]
505002290[−1,700, 2,087]
2005001632[−1,072, 2,559]
20010001274[−1,776, 2,095]
randomized benchmark1,794—

Seed alone moves the estimate by USD 562; trees and draws together by USD 477.

The estimate is not merely attenuated but unstable, and the instability exceeds the tuning sensitivity. Changing only the random seed moves it by USD 562; quadrupling the trees and doubling the draws moves it by USD 477. Any single run is close to uninformative about its own answer. What is robust is the direction and scale of the miss: across every configuration the estimate lands between USD 154 and USD 716, never within a thousand dollars of the benchmark.

The mechanism is the one the previous example exposed. The prognostic signal here moves 1978 earnings by thousands while the treatment effect is under two thousand dollars. A regularized ensemble spends its splits where the variance is, and the treatment indicator is left to soak up what remains. The shrinkage that makes the model good at predicting is what destroys it as an estimator of a causal contrast. Hahn, Murray & Carvalho (2020) named this regularization-induced confounding.

The per-unit effects show the same thing from another angle. A response surface returns a posterior for every unit’s τi\tau_i, which nothing else in the observational group can do — but only if that variation is identified rather than imagined. Here it is systematic, with an R2R^2 of 0.67 on the covariates, and tiny beside its own uncertainty, at a signal-to-noise ratio of 0.12. Ranking individuals by estimated benefit would be reading sampling error.

Bayesian Causal Forests

The repair is structural. Instead of one surface in which treatment is a covariate like any other, write two:

E[Y∣X,W]  =  μ(X, e^(X))⏟prognostic  +  τ(X)⏟treatment effect⋅WE[Y \mid X, W] \;=\; \underbrace{\mu\big(X,\ \hat e(X)\big)}_{\text{prognostic}} \;+\; \underbrace{\tau(X)}_{\text{treatment effect}} \cdot W

Two changes do the work. The treatment effect gets its own surface, with its own priors, so shrinking the prognostic fit no longer shrinks the effect. And the estimated propensity score e^(X)\hat e(X) is admitted to the prognostic part, where it can absorb the selection that would otherwise leak into τ\tau. There is a bonus: τ(X)\tau(X) is the treatment effect, so the ATT reads straight off its posterior — no counterfactual evaluation needed at all.

BCF returns USD 1,328. The obvious objection is that it simply had more capacity, and the earlier sweep gave that objection teeth by showing single-surface estimates that sometimes climbed with tree count. So both models were run at a matched budget of 75 trees, two seeds each.

model, 75 trees eachseed 11seed 12mean
single surface254724489
BCF 50 + 251,3231,1401,232
randomized benchmark1,794

At equal capacity the gap is USD 743. Splitting the surface buys what buying more trees does not — and the apparent “more trees fixes it” signal in the earlier sweep was itself seed noise, which is precisely why the matched comparison uses two seeds rather than one.

The per-unit effects improve in kind as well: signal-to-noise goes from 0.12 to 1.12. Under the single surface the between-unit variation was a ninth of the uncertainty within any one unit; now they are comparable. Comparable is not resolved. A ratio near 1 says the model expresses heterogeneity that is not obviously imaginary, which is a weaker and more defensible claim than a ranking.

Every Estimator Against the Benchmark

estimatorATTvs benchmark
naive difference in means−635−135%
BART, single surface283−84%
Mahalanobis matching719−60%
TMLE1,199−33%
AIPW1,225−32%
entropy balancing1,273−29%
IPW inverse-propensity weighting1,316−27%
BART, causal forests1,328−26%
regression adjustment1,548−14%
NN-propensity matching1,792−0%
randomized benchmark1,794—

The two BART rows are the same algorithm on the same data with the same priors, differing only in whether the treatment effect has a surface of its own.

Where the propensity score went. The first example proved a coherent Bayesian cannot use it as a weight. Here it returns as a predictor, admitted to the prognostic surface — the one door the likelihood leaves open — and that alone is worth a factor of four on this dataset. The design information shown to be unusable in one role turns out to be decisive in another.

And the caution this group keeps earning. Even repaired, BCF sits 26% below the benchmark, in the same band as every doubly-robust method here — while plain nearest-neighbour propensity matching, the least sophisticated estimator in the collection, lands within USD 2 of the truth. One dataset is not a ranking of methods. But on this problem sophistication has bought better properties rather than better answers, which is the third time this group has reached that conclusion by a different route.

Where this sits

This closes the Bayesian selection-on-observables group. The estimators it is measured against come from matching and weighting and modern balancing; the untestable assumption underneath all of them is priced in Confounding as a Parameter.

Notebook

Python A single BART response surface on LaLonde with counterfactuals by tree evaluation, a seed-and-tuning sweep separating instability from attenuation, the per-unit effects and their signal-to-noise, then Bayesian Causal Forests as two BART components with the propensity score in the prognostic part — and a matched 75-tree comparison to rule out capacity as the explanation (Python, PyMC, pymc_bart) Open →

Downloads

lalonde_obs.csv 185 NSW treated against 429 CPS controls — the sample every model here is fitted to lalonde_exp.csv The randomized NSW sample supplying the USD 1,794 benchmark

References