[Feature Request]:

[professor-cathy-chen]

Describe the feature you want to propose or implement @SvenKlaassen @PhilippBach

I would like to contribute a new model class for sample selection models estimated with a control function (Heckman-type) correction, where the target parameter is the sample selection bias coefficient itself.

Model. An outcome equation and a participation (selection) equation,

Y_i = g_0(X_i) + eps_i           (Y_i observed only if D_i = 1)
D_i = 1{ Z_i'beta_0 + v_i > 0 },  Z_i = (X_i, U_i)

with jointly normal (eps_i, v_i). X_i are the covariates of the outcome equation and U_i are additional variables that shift participation but are excluded from the outcome equation. Under joint normality, E[eps_i | X_i, Z_i, D_i = 1] = theta_0 * h_i with the inverse Mills ratio h_i = phi(Z_i'beta_0) / Phi(Z_i'beta_0), so the observed outcome equation becomes

Y_i = g_0(X_i) + theta_0 * h_i + u_i,   E[u_i | X_i] = 0

and the parameter of interest is theta_0 = **sigma_{eps,v}.**

Why this is useful. theta_0 determines whether outcomes observed in a self-selected subsample can be used to say anything about the units that do not report. If theta_0 != 0, imputing outcomes for non-reporters from a model fitted on reporters is systematically biased. The class provides both a sqrt(n)-consistent estimate of the bias and a test of its absence. Typical applications are voluntary corporate disclosure (our motivation is carbon emissions disclosure), survey non-response, and wage equations with participation decisions — settings where selection is on unobservables and the researcher has covariates that shift participation but not the outcome.

Relation to the existing DoubleMLSSM. These are different estimands with different identifying assumptions, not competing implementations of the same thing. DoubleMLSSM targets an ATE under outcome attrition using inverse probability weighting, under MAR or under nonignorable nonresponse with an instrument. The proposed class targets the selection bias coefficient in a structural outcome equation via a control function, and identifies coefficients of the outcome equation under selection on unobservables. In practice they answer different questions and would typically be used by different literatures.

Propose a possible solution or implementation

A new class following the standard OOP structure — only the nuisance functions and the score components are specified, everything else inherited from the abstract base class.

Linear Neyman-orthogonal score, evaluated on the selected subpopulation {D = 1}:

psi(W; theta, eta) = (Y - g(X) - theta*h) * (h - m(X))
psi_a = -h * (h - m)
psi_b = (Y - g) * (h - m)

Nuisance functions:

learner	type	target
ml_pi	classifier	pi(Z) = P(D = 1 | Z)
ml_g	regressor	g(X) = E[Y | X, D = 1]
ml_m	regressor	m(X) = E[h | X, D = 1]

Three implementation points that I think are the interesting ones:

(a) The inverse Mills ratio is link free. Since pi_i = Phi(Z_i'beta_0), the index is recoverable from the propensity score through the probit link, giving h_i = phi(Phi^{-1}(pi_i)) / pi_i. Under joint normality this is exactly the inverse Mills ratio, so any probabilistic classifier can be used for the selection equation and a probit is just the special case. This is what makes a control function correction usable inside the ML-learner interface at all. Note that h diverges as pi -> 0, so strict overlap is a genuine identification requirement here rather than only numerical hygiene; the class uses the standard PSProcessorConfig / init_ps_processor machinery for trimming.

(b) Nested cross-fitting. h_hat is a generated regressor that is simultaneously the target of m. Each training fold is therefore split in half: the first half estimates the selection equation, the second half (restricted to D = 1) estimates g and m. This is the same pattern as the nonignorable score of DoubleMLSSM. Exposed as nested_cross_fitting=True (default), which can be switched off.

(c) Zero padding gives the effective sample size for free. The score elements are set to zero for D_i = 0. With n = sum_i D_i observed outcomes, the base class variance estimator mean(psi^2) / J^2 / N then reduces algebraically to sigma^2_theta / n with -J_0 = E_n[(h - m)^2 | D = 1]. No override of the variance estimation is needed. There is a unit test asserting this identity against a manual implementation, in both languages.

Status of the implementation. Working and tested against DoubleML 0.11.4, with an R6 twin against DoubleML 1.0.2:

• 60 passing pytest tests. The core tests compare the coefficient, standard error, both score elements, all four cross-fitted nuisance vectors and the score-test statistic against a manual line-by-line implementation of the algorithm, across n_folds in {3, 5} and with/without nested cross-fitting. Plus the zero-padding variance identity, a theta_0 = 0 centering check, and a full exception battery.
• black --check clean, ruff clean.
• Monte Carlo evidence that the orthogonalization is doing the work: bias falls as -0.208 -> -0.108 -> -0.056 for N = 1000, 2000, 4000 (roughly O(1/N), consistent with the second-order product-of-errors term), while the naive two-step plug-in sits at +0.33, +0.31, +0.30 and does not improve. At theta_0 = 0 the estimator is centred (bias -0.004) and the score test has empirical size 0.08 against a nominal 5%, with power 0.84 at theta_0 = -0.25 and 1.00 at -0.50.
• An executed example-gallery notebook in the style of py_double_ml_ssm.ipynb, built entirely on a simulated DGP (see the last field).

Open questions I'd like your view on before writing PR code:

1. Naming. I have been using DoubleMLSSCF (Sample Selection, Control Function). I'm not attached to it, and SSCF vs the existing SSM may be easy to misread in review. DoubleMLHSM (Heckman selection model) or DoubleMLSSB (sample selection bias) are alternatives. Your call.
2. Data backend. Currently plain DoubleMLData with d_cols holding the selection indicator, x_cols the outcome covariates and z_cols the exclusion restrictions. This needs no new data class, but it overloads "treatment" to mean "selection indicator". The alternative is extending DoubleMLSSMData to allow s_col with no treatment variable, which is semantically honest but touches an existing class. Which do you prefer?
3. Scope. Should the two extra methods above stay on the class or move to the example gallery?
4. Sensitivity analysis. Currently _sensitivity_implemented = False. I'm genuinely unsure whether a Chernozhukov-style sensitivity analysis is well defined for this parameter and would welcome a view.
5. Tuning. _nuisance_tuning is implemented; _nuisance_tuning_optuna raises NotImplementedError. Is that acceptable given the nested structure, as with DoubleMLSSM?
6. R twin. Open the doubleml-for-r PR at the same time, or after the Python class is settled?

Did you consider alternatives to the proposed solution. If yes, please describe

Yes, four.

Reusing DoubleMLSSM with a new score string. Rejected. The nuisance set is different (this model needs the projection m(X) = E[h | X, D = 1], which has no counterpart in an IPW score), the estimand is a structural coefficient rather than an ATE, and the identifying assumption is joint normality of the error terms rather than a conditional independence or instrument condition. Bolting a third score onto DoubleMLSSM would make that class harder to reason about for the model it already serves.

A callable score on DoubleMLPLR. Technically possible — h can be treated as an endogenous-looking regressor and the score is linear. Rejected because the user would have to construct the cross-fitted inverse Mills ratio themselves, get the nested cross-fitting right by hand, and remember to zero-pad off {D = 1}. Those are precisely the three things that are easy to get wrong, and they are the reason a dedicated class is worth having.

Estimating theta_0 by naive two-step Heckman with ML first stages. This is the obvious baseline and it is what the class exists to replace. It is inconsistent here: without partialling out E[h | X], the regularization bias of the high-dimensional first stages contaminates the estimate. In the simulation it is attenuated to less than half the true value and, unlike the orthogonal estimator, does not improve with sample size (bias +0.33, +0.31, +0.30 at N = 1000, 2000, 4000). It also has a tighter sampling distribution than the orthogonal estimator — tightly wrong — which is a good cautionary illustration for the docs.

A new data class vs. reusing DoubleMLData. Considered both; this is open question 2 above. I implemented the lighter option so that the contribution touches as little existing code as possible, but I don't think it is obviously the right one.

On scope: the paper also contains an adaptive kernel group lasso for post-selection estimation of the outcome equation. I deliberately propose to leave it out of this contribution. It is a variable selection device applied after theta_hat is available, not part of the orthogonal score, and bundling it would roughly double the review surface for no benefit to the model class. I'd ship it with the example notebook and propose it separately only if there's appetite.

Comments, context or references

Data availability. The method is developed in a working paper (Chen, Lioui and Scaillet, Green Silence: Double machine learning carbon emissions under sample selection bias). The empirical data are licensed from a commercial vendor and cannot be redistributed, so the contribution ships a simulated DGP and the example notebook is built entirely on it — in the same spirit as the existing py_double_ml_ssm example, which uses the Bia, Huber and Lafférs design. The DGP generator (characteristics with Toeplitz covariance, Gaussian kernel features on characteristic-sorted portfolio grids, group-sparse coefficients, an active set in the selection equation disjoint from the outcome active set, and an intercept calibrated to a target participation rate) would be contributed alongside the class.

Compatibility. Developed against DoubleML 0.11.4 (Python) and 1.0.2 (R). Happy to rebase onto main.

Proposed file layout, if the design is agreed (sample selection models under irm/, next to ssm.py):

doubleml/irm/sscf.py                       # model class
doubleml/irm/datasets/dgp_green_silence_data.py
doubleml/irm/tests/test_sscf.py
doubleml/irm/tests/test_sscf_exceptions.py
doubleml/irm/tests/_utils_sscf_manual.py   # manual reference implementation

plus exports and API doc stubs; the notebook goes to doubleml-docs as a separate PR.

References

• Heckman, J. J. (1979), Sample selection bias as a specification error, Econometrica 47(1), 153-161.
• Chernozhukov, Chetverikov, Demirer, Duflo, Hansen, Newey and Robins (2018), Double/debiased machine learning for treatment and structural parameters, The Econometrics Journal 21(1), C1-C68.
• Bia, Huber and Lafferty (2023), Double machine learning for sample selection models, JBES — the basis of the existing DoubleMLSSM, included here for contrast.
• Yuan and Lin (2006), Model selection and estimation in regression with grouped variables, JRSS-B 68(1), 49-67 — background for the post-selection step, not part of this proposal.

I'm happy to open a draft PR straight away if the design looks reasonable, or to adjust it first based on your answers to the questions in the previous field.

[SvenKlaassen]

Thanks you!
This looks like a very well-thought-out contribution.
Regarding you open points:

1. Class Name: DoubleMLSSCF is fine for me, even if it is close to DoubleMLSSM.
2. Data Backend: I generally would prefer to extend DoubleMLSSMData (using s_cols for selection) than overloading d_cols in DoubleMLData. Keeping treatment ($D$) and selection ($S$) semantically distinct makes user code much clearer.
3. Scope: Could you specify what exactly you are refering to?
4. Sensitivity Analysis: Setting _sensitivity_implemented = False fine. I would have to check whether this is possible, but that could also be added at a later stage.
5. Hyperparameter Tuning: Raising NotImplementedError for _nuisance_tuning_optuna makes complete sense.
6. R Twin: I think opening doubleml-for-py PR first to settle the API and core design before opening the doubleml-for-r PR is best approach.

Further, with the ML-based index estimation wouldn't it be clearer to actually write

$$D_i = \mathbf{1} \{ f_0(Z_i) + v_i > 0 \}$$

where $f_0(Z_i)$ is an unknown non-linear function and $h_i = \frac{\phi(\Phi^{-1}(\pi(Z_i)))}{\pi(Z_i)}$ is the generalized Inverse Mills Ratio constructed from the propensity score $\pi(Z_i) = P(D_i = 1 \mid Z_i) = \Phi(f_0(Z_i))$.

You should add your paper (Chen, Lioui and Scaillet, Green Silence: Double machine learning carbon emissions under sample selection bias) to the References for the PR :)

Looking forward to seeing the draft PR!

[professor-cathy-chen]

@SvenKlaassen Thanks for the detailed feedback — very helpful. I've pushed the revisions in commits 6c975d3 and d75dd9a. Let me respond point by point.

On your six points

1. Class name (DoubleMLSSCF). Agreed, keeping it.

2. Data backend (DoubleMLSSMData with s_col). I agree with the principle — keeping treatment and selection semantically distinct is clearly better than overloading d_cols.

One complication I ran into when implementing it. In this model there is no treatment variable at all — the data are (Y, S, X, U), where S is the selection indicator and U are the exclusion restrictions. Moving the selection indicator into s_col leaves d_cols with nothing to put in it, and on 0.11.4 d_cols is a required positional argument:

DoubleMLSSMData(df, y_col='y', d_cols=None, x_cols=..., z_cols=..., s_col='s')
# TypeError: The treatment variable(s) d_cols must be of str or list type.

Three options as I see them:

(a) Relax d_cols in DoubleMLSSMData to accept None (with n_treat = 0). Cleanest semantically and a small change, but it touches a class DoubleMLSSM relies on, so it's your call.

(b) A thin subclass DoubleMLSSCFData(DoubleMLSSMData) that permits d_cols=None and inherits everything else. Fully contained in this contribution, at the cost of one extra data class.

(c) Require a dummy treatment column — I'd rather not, since it reintroduces the overloading you're rightly trying to avoid.

I'd lean towards (a), with (b) as the low-risk fallback. Happy to implement either. I'm holding off the refactor until you've weighed in, since the choice changes the constructor, _check_data, every test fixture and the notebook.

Also a small thing: on 0.11.4 the argument is s_col (singular, str) rather than s_cols. Is s_cols something planned for a future release that I should target instead?

3. Scope. Sorry that was too vague. Two separate things:

(i) The kernel group lasso. The paper also contains an adaptive kernel group lasso for variable selection in the outcome equation after theta_hat is available. It is not part of the orthogonal score and the model class does not depend on it. It is currently in the branch as doubleml/sscf/kgl.py. My inclination is to drop it from this PR and ship it with the example notebook instead, to keep the review surface to the model class. Happy to keep it if you'd find it useful in the package.

(ii) Two extra public methods on the class, beyond the base API:

• selection_bias_test(level=0.05) — orthogonal score test of H₀: θ₀ = 0, evaluating the score at θ = 0 with the cross-fitted nuisances; asymptotically χ²₁.
• bias_quantification() — the implied bias term θ̂ĥᵢ, the bias-corrected outcomes, and summary statistics split by selected / non-selected.

There is precedent in robust_confset() on DoubleMLIIVM, but if you'd prefer a lean class that sticks to the base API, both move to the notebook without loss.

4. Sensitivity analysis. Agreed, _sensitivity_implemented = False for now.

5. Hyperparameter tuning. Good, NotImplementedError for optuna stays.

6. R twin. Agreed — Python first, R after the API is settled. The R6 class is ready and waiting.

On the f₀(Z) notation

You're right, and I've adopted it in the latest push. The derivation only needs two things:

• v_i | Z_i ~ N(0, 1), which gives the probit link π(Z) = Φ(f₀(Z));
• joint normality of (ε_i, v_i), which gives E[ε | v] = θ₀v.

Linearity of the index is a parametric convention from Heckman (1979) and is not an identifying assumption. For any f₀,

h_i = φ(f₀(Z_i)) / Φ(f₀(Z_i)) = φ(Φ⁻¹(π(Z_i))) / π(Z_i)

and the code has always implemented the general version. Writing Z'β in the docs was understating the scope of the contribution.

I've also added a nonlinear_index=True option to the DGP so the test suite can demonstrate the class works where a probit would fail (a linear fit to the true index explains only R² = 0.36 under this option). Two new tests cover it (the suite is now 63 passing).

Calibration matters — a caveat worth documenting

Adopting f₀(Z) naturally invites users to plug in flexible learners. One thing I found while testing this is that the generalized IMR construction is link free but not calibration free. Because Φ⁻¹ steepens in the tails, θ̂ is sensitive to the calibration of π̂, not only to its accuracy. This is a stronger requirement than IPW-type scores impose.

Small experiment on the shipped DGP (θ₀ = −0.5, N = 4000, 6 replications, linear index so the specification is not the issue):

ml_pi	bias in θ̂
ℓ₁-logistic (correctly specified)	−0.003
HistGradientBoosting, raw	+0.275
HistGradientBoosting + Platt/sigmoid calibration	+0.032
HistGradientBoosting + isotonic calibration	+0.069

So an uncalibrated tree ensemble is badly biased even when the index is linear, and CalibratedClassifierCV with Platt scaling removes most of it. I've documented this in the class docstring (Notes section) and I'm flagging it here because it's the natural follow-up to your f₀ suggestion.

These numbers are only 6 replications — I'll run them properly before putting them in the docs. But the direction and magnitude are clear.

Should the calibration support in PSProcessorConfig be the recommended mechanism for this, or is CalibratedClassifierCV the right user-facing recommendation?

Reference

Added per your suggestion:

Chen, C. Y.-H., Lioui, A. and Scaillet, O. (2025), Green Silence: Double Machine Learning Carbon Emissions under Sample Selection Bias. Available at SSRN: https://ssrn.com/abstract=5368583

What's in this push

• Notation: Z'β → f₀(Z) throughout docstrings and class description.
• New: nonlinear_index option in make_green_silence_data, with two new tests (63 total, all passing).
• New: Notes section in the class docstring on choosing ml_pi (calibration).
• Fixed: deprecated n_alphas argument removed from all LassoCV calls.
• Fixed: test imports now use from doubleml.sscf import ... (removed the doubleml_sscf shim).
• Added: DoubleMLSSCF exported from doubleml/__init__.py.
• SSRN reference added.

Holding off on the data-backend refactor until you've weighed in on (a) vs (b) above.

[professor-cathy-chen]

@SvenKlaassen Thanks for the clear guidance. All three points are now addressed in the latest push (57a7cc3, CI green):

1. Data backend → DoubleMLSSMData with s_col: selection indicator moved from d_cols to s_col; d_cols made optional in DoubleMLSSMData (defaults to None, n_treat = 0). Minimal change to ssm_data.py (3 deletions). Existing SSM tests unaffected.

2. Kernel group lasso dropped: kgl.py and test_kgl.py removed. Will ship with the example notebook in a separate PR on doubleml-docs after this is merged.

3. Output-structure tests added for selection_bias_test() and bias_quantification(): key/type/shape checks and pre-fit() error handling.

Also in this push:

• f_0(Z) notation adopted throughout (per your earlier suggestion)
• nonlinear_index option added to the DGP, with tests
• ml_pi calibration caveat documented in the class docstring
• SSRN reference added
• Deprecated n_alphas removed
• notebooks/ and scripts/ removed from the package

10 changed files, 1,599 additions, 3 deletions. Ready for review whenever you are.