This vignette documents the method implemented by oal()
at the level of formulas, explains the one nontrivial piece of numerical
machinery (the exact penalty-scale correction for glmnet),
records the provenance of every default (which choices are Shortreed
& Ertefaie’s, which are Schnitzer’s, which are the package’s),
documents the differences from the author’s own legacy code, situates
oalasso among related software, and closes with an honest
section on nonlinear extensions. Throughout, “S&E” is Shortreed
& Ertefaie (2017) and “the rejoinder” is Jones, Ertefaie &
Shortreed (2023), whose companion code is the reference implementation
this package reproduces.
For subjects \(i = 1, \dots, n\): covariates \(X_i\) (a \(p\)-column numeric matrix after full dummy expansion of factors), binary treatment \(A_i \in \{0, 1\}\), and outcome \(Y_i\). The propensity score model is logistic, \(e(X_i) = \Pr(A_i = 1 \mid X_i) = \operatorname{expit}(\alpha_0 + X_i^\top \alpha)\). Covariates are conceptually partitioned into confounders (predict \(A\) and \(Y\)), pure outcome predictors (predict \(Y\) only), instruments (predict \(A\) only), and noise. All fitting is done on the standardized matrix \(X^s\) (below); coefficients are back-transformed to the original scale for output.
An unpenalized regression of the outcome on treatment and all PS-model covariates, fit on the full sample (both arms pooled — the S&E/rejoinder convention):
\[ Y_i = \eta A_i + \sum_{j=1}^{p} \beta_j X^s_{ij} + \varepsilon_i \quad\text{(via } \texttt{lm} \text{ for } \texttt{gaussian()}\text{, } \texttt{glm} \text{ for } \texttt{binomial()}\text{)}. \]
Only the covariate coefficients \(\hat\beta_1, \dots, \hat\beta_p\) are used;
the treatment and intercept coefficients are discarded, and the
treatment is never penalized anywhere. A binomial() outcome
family is allowed with a one-time warning: S&E’s validating
simulations used continuous outcomes; the GLM extension of the theory is
Baldé (2025, doi:10.1016/j.spl.2025.110379).
Because OLS coefficients are almost surely nonzero, an exactly-zero
or NA \(\hat\beta_j\) from
the internal model signals a data pathology (rank deficiency, aliased
columns) that would silently corrupt the penalty scaling below —
oal() therefore stops with an error naming the offending
columns rather than continuing.
\[ \hat\alpha(\lambda_n, \gamma) \;=\; \arg\min_\alpha \; -\ell_n(\alpha; A, X^s) \;+\; \lambda_n \sum_{j=1}^{p} \hat w_j\, |\alpha_j|, \qquad \hat w_j = |\hat\beta_j|^{-\gamma}, \]
where \(\ell_n\) is the total (not per-observation) binomial log-likelihood. The mechanism: \(\hat\beta_j\) is \(\sqrt n\)-consistent, so for outcome-related covariates \(\hat w_j = O_p(1)\), while for outcome-unrelated covariates (instruments, noise) \(|\hat\beta_j| = O_p(n^{-1/2})\) and the penalty weight explodes at rate \(n^{\gamma/2}\) — the outcome model decides which coefficients get an exploding penalty in the treatment model. S&E’s Theorem 1 gives the oracle-type result under the rate conditions
\[ \lambda_n / \sqrt{n} \to 0, \qquad \lambda_n\, n^{\gamma/2 - 1} \to \infty, \qquad \gamma > 1 : \]
outcome-unrelated coefficients are set to zero with probability tending to one, and the retained coefficients are \(\sqrt n\)-asymptotically normal with the oracle variance. Note the second condition is heavier than the standard adaptive-lasso rate (\(\lambda_n n^{(\gamma-1)/2} \to \infty\)): an instrument has a genuinely nonzero treatment-model coefficient, so the treatment likelihood pulls it in, and only a heavier penalty overrides that pull. \(\gamma > 1\) is exactly what makes the window between the two conditions nonempty. These guarantees are conditional on the linear outcome model being correctly specified — see the final section.
lambda supplies exponents \(\delta\), with \(\lambda_n = n^\delta\) and \(n\) always the number of rows of the model
matrix. The default is S&E’s grid:
\[ \delta \in \{-10, -5, -2, -1, -0.75, -0.5, -0.25, 0.25, 0.49\}. \]
Each \(\delta\) is paired with a penalty exponent (never crossed, so there is no silent grid explosion):
\[ \gamma_\delta = 2\,(g_{cf} - \delta + 1), \qquad g_{cf} = 2 \text{ (the \texttt{gamma.factor} argument)}, \]
which is the value solving \(\lambda_n \,
n^{\gamma/2 - 1} = n^{g_{cf}}\) — every grid point sits on the
same divergence rate for the second oracle condition while \(\lambda_n\) itself ranges from negligible
to nearly \(\sqrt n\). This is the
S&E / rejoinder convention. Supplying a scalar gamma
(e.g. gamma = 2.5) instead fixes \(\gamma\) and crosses it with the whole
\(\delta\) grid — Schnitzer et al.’s
(2025) configuration; their raw lambda grid converts to exponents via
\(\delta = \log(\lambda_n)/\log(n)\). A
guard warns if lambda values look like raw penalties rather
than exponents (max(lambda) > 3).
For each grid point, the fitted scores are clipped to
clip (default \([0.01,
0.99]\)) and converted to IPW weights at the requested
estimand:
\[ \text{ATE:}\quad w_i = \frac{A_i}{e_i} + \frac{1 - A_i}{1 - e_i}, \qquad\qquad \text{ATT:}\quad w_i = A_i + (1 - A_i)\,\frac{e_i}{1 - e_i}. \]
The weighted absolute mean difference is then, over all \(p\) standardized columns (including excluded ones),
\[ \mathrm{wAMD} \;=\; \sum_{j=1}^{p} |\hat\beta_j| \cdot \left| \frac{\sum_i w_i X^s_{ij} A_i}{\sum_i w_i A_i} - \frac{\sum_i w_i X^s_{ij} (1 - A_i)}{\sum_i w_i (1 - A_i)} \right|, \]
and the selected grid point is the argmin — the first minimum within
a \(10^{-9}\) relative tolerance, in
the documented grid order (lambda2 ascending with 0 first,
then \(\delta\) ascending), so ties
prefer plain OAL and then the smallest penalty. This formula is
implemented once, natively, in the exported oal_wamd(); it
is never delegated to cobalt (formula fidelity), while the
display balance table ($balance,
bal.tab()) uses cobalt’s native conventions, exactly as in
psAve.
Known issues of the criterion, stated plainly: wAMD is a balance
heuristic — there is no theorem that it lands \(\lambda_n\) inside the oracle window; it
reuses the same \(|\hat\beta_j|\) as
the penalty and so inherits the same outcome-model-misspecification
vulnerability; and it is unstable under near-positivity violations,
where a few huge weights dominate the sums — the failure mode the
rejoinder identifies in correlated designs, and the motivation for GOAL.
print()/summary() therefore warn when more
than 5% of the scores sit at a clipping bound.
With the default refit = FALSE, the score is \(\operatorname{expit}\) of the
penalized coefficients (S&E/rejoinder). With
refit = TRUE (the Schnitzer variant, adapted from their
longitudinal setting and provenance-labeled as such), an unpenalized
logistic regression is refit on the selected covariates at every grid
point, and both the wAMD and the final score come from the refit;
per-grid-point convergence is recorded in
path$refit.converged. All scores are clipped to
clip before the wAMD and in the returned
$ps.
method = "goal" (Baldé, Yang & Lefebvre 2023) adds a
ridge term to the objective:
\[ \hat\alpha(\lambda_n, \gamma, \lambda_2) \;=\; \arg\min_\alpha \; -\ell_n(\alpha) \;+\; \lambda_n \sum_j \hat w_j |\alpha_j| \;+\; \lambda_2 \sum_j \alpha_j^2 . \]
The \(\ell_2\) term supplies an elastic-net grouping effect and numerical stability, so correlated confounders are selected together rather than arbitrarily, and it behaves like propensity score truncation under near-positivity — the rejoinder’s own remedy for OAL’s failure mode. Baldé (2025) proves the oracle property under the OAL rate conditions plus \(\lambda_2/\sqrt n \to 0\).
The solver uses the Zou–Hastie data-augmentation trick, exactly as in the rejoinder code and as described by Baldé (2025): for \(\lambda_2 > 0\), augment
\[ X^s_{\text{aug}} = \begin{pmatrix} X^s \\ \sqrt{\lambda_2}\, I_p \end{pmatrix}, \qquad A_{\text{aug}} = (A^\top, 0_p^\top)^\top, \]
fit the lasso on the augmented data, then rescale
all coefficients — including the intercept —
by \((1 + \lambda_2)\) (Baldé’s
official code:
expit(cbind(1, X) %*% (1 + lambda2) * coef)); scores are
computed from the original \(n\) rows
only. On augmented grid points the raw penalty constant uses the
augmented sample size, \(\lambda_n = (n + q)^\delta\) with \(q = p\) augmentation rows (the author’s
adaptive.lasso(lambda = n.q^(il)),
n.q = n + q). Grid points with \(\lambda_2 = 0\) run the plain OAL code
path — no augmentation, \(\lambda_n =
n^\delta\) — because even zero-valued augmented rows would
contribute intercept terms to the binomial likelihood and break exact
nesting. (Baldé’s script literally augments with \(q\) zero rows and keeps the \((n+q)^\delta\) base even at \(\lambda_2 = 0\); oalasso
deliberately prefers exact nesting there.) lambda2 is
selected jointly with \((\delta,
\gamma)\) by the same wAMD: the flat first-minimum over the \(\lambda_2\)-outer \(\times\) \(\delta\)-inner grid, equivalent to Baldé’s
nested rule (per-\(\lambda_2\) minima
over \((\delta, \gamma)\), then the
first minimum over \(\lambda_2\)).
Provenance. The default grid
lambda2 = c(0, 10^c(-2, -1.5, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 1))
is the author’s published grid (Baldé 2025 supplement,
“taken from Zou and Hastie (2005)”), verified against the official
lqa-based GOAL code on 2026-07-02; the printed provenance line says so.
Supplying your own lambda2 switches the label to
“user-specified lambda2 grid”.
Clipping. Baldé’s reference code runs
unclipped IPW weights (\(1/e\)
and \(1/(1-e)\) with no truncation).
The oalasso default clip = c(0.01, 0.99) is a
Shortreed–Ertefaie-lineage safety choice (the suite constant shared with
psAve); to reproduce the reference behavior, effectively
disable clipping with clip = c(1e-12, 1 - 1e-12).
\(X^s = \texttt{scale}(X, \texttt{TRUE},
\texttt{TRUE})\) is computed once, with training
statistics, over all columns including dummies (the
rejoinder convention). The outcome model and every glmnet
fit use the same \(X^s\) with
standardize = FALSE. This matters because adaptive penalty
weights are scale-dependent: letting glmnet standardize internally, or
scaling different subsets separately, changes \(|\hat\beta_j|^{-\gamma}\) and hence the
selections. Output coefficients are back-transformed (\(\alpha_j = \alpha^s_j / s_j\); intercept
adjusted by \(-\sum_j \alpha^s_j\, c_j /
s_j\)), and both scales are returned (coefficients,
coefficients.std); the centers and scales are stored in
info.
glmnet alone is
not enough)This is the core correctness requirement of the package, and the
reason a naive glmnet(penalty.factor = ...) translation of
OAL does not reproduce the published method.
Two scale mismatches stand between glmnet and the OAL objective:
nvars:
\(\tilde v_j = v_j / \overline{v}\)
where \(\overline v = \tfrac1p \sum_j
v_j\). The documented consequence (“the penalty factors are
internally rescaled to sum to nvars, and the lambda sequence will
reflect this change”) is that the lambda you pass is silently
reinterpreted. With adaptive weights this is not benign: one near-zero
\(\hat\beta_j\) inflates \(\overline v\) enormously and thereby
relaxes the penalty on every other covariate.Both are undone at once by evaluating each grid point at
\[ s \;=\; \overline{v}_{\text{eff}} \cdot \frac{\lambda_n}{n}, \qquad \overline{v}_{\text{eff}} = \frac{\sum_{j:\,v_j < \infty} v_j \;+\; \#\{j : v_j = \infty\}}{p}, \]
via
coef(fit, s = s, exact = TRUE, x = Xs, y = A, penalty.factor = pen)
— exact = TRUE re-solves at exactly this \(s\) rather than interpolating along the
path, and requires re-supplying the data and penalty factors. The
effective mean \(\overline{v}_{\text{eff}}\) reflects
glmnet’s actual behavior with infinite factors: an infinite \(v_j\) (\(\hat b_j
= 0\) under the user hook below) is converted to
exclude — a hard drop — and then internally reset to 1
before the sum-to-nvars rescaling, so each excluded
column contributes 1 (not nothing) to the rescaling mean on the full
matrix. (The plain \(\operatorname{mean}\{v_j
: v_j < \infty\}\) would be correct only on a reduced matrix
from which excluded columns are physically removed; oalasso
keeps the full matrix and uses \(\overline{v}_{\text{eff}}\), which reduces
exactly to \(\operatorname{mean}(v)\)
when all factors are finite.) Under GOAL augmentation both constants
switch to the augmented sample size: \(\lambda_n = (n + q)^\delta\) and the \(1/n\) becomes the augmented per-observation
weight \(1/(n + q)\), \(q = p\) augmentation rows.
The correction is always on — it is not a knob — and
its acceptance criterion is not the formula but the
Karush–Kuhn–Tucker conditions of the target objective: the test
suite verifies, at the exact solution of every grid-point type, that
\(|\partial \ell / \partial \alpha_j| \le
\lambda_n \hat w_j\) for zero coefficients with equality for
nonzero ones. This recipe and its KKT verification are ported from the
rejoinder’s companion code (jmiahjones/OAL-rejoinder-2022,
tests/test-glmnet.R); a drift-guard test pins the behavior
against the installed glmnet version.
A minimal illustration of the machinery on simulated data
(x5, x6 are instruments):
library(oalasso)
set.seed(20260702)
n <- 300
X <- matrix(rnorm(n * 6), n, 6, dimnames = list(NULL, paste0("x", 1:6)))
A <- rbinom(n, 1, plogis(X %*% c(1, 1, 0, 0, 1, 1)))
Y <- as.numeric(2 * A + X %*% c(0.6, 0.6, 0.6, 0.6, 0, 0) + rnorm(n))
sim <- data.frame(A, Y, X)
fit <- oal(A ~ x1 + x2 + x3 + x4 + x5 + x6, data = sim, outcome = ~ Y)
## the path records, per grid point: lambda2, delta, lambda.n = n^delta,
## the paired gamma = 2*(2 - delta + 1), the exact glmnet s actually used,
## the wAMD, and the selected-set size
fit$path[, c("delta", "lambda.n", "gamma", "s", "wamd", "n.selected")]
#> delta lambda.n gamma s wamd n.selected
#> 1 -10.00 1.693509e-25 26.00 9.012153e+19 0.07692437 4
#> 2 -5.00 4.115226e-13 16.00 8.995054e+13 0.07691585 4
#> 3 -2.00 1.111111e-05 10.00 2.386841e+10 0.07691109 4
#> 4 -1.00 3.333333e-03 8.00 1.553156e+09 0.08201828 4
#> 5 -0.75 1.387264e-02 7.50 7.853034e+08 0.10329190 4
#> 6 -0.50 5.773503e-02 7.00 3.972424e+08 0.16404414 3
#> 7 -0.25 2.402811e-01 6.50 2.010352e+08 0.32596912 3
#> 8 0.25 4.161791e+00 5.50 5.155981e+07 0.70276946 0
#> 9 0.49 1.636023e+01 5.02 2.684984e+07 0.70276946 0
all.equal(fit$path$gamma, 2 * (2 - fit$path$delta + 1))
#> [1] TRUE| Setting | Default | Source |
|---|---|---|
| Penalty objective, adaptive weights \(\hat w_j = \lvert\hat\beta_j\rvert^{-\gamma}\) | — | S&E (2017) |
| Outcome model \(Y \sim A + X\), full sample, OLS/GLM | family = gaussian() |
S&E (2017) / rejoinder code |
| \(\lambda\) grid (exponents \(\delta\)) | c(-10, -5, -2, -1, -0.75, -0.5, -0.25, 0.25, 0.49) |
S&E (2017) |
| \(\gamma\) pairing \(2(g_{cf} - \delta + 1)\), \(g_{cf} = 2\) | gamma = NULL, gamma.factor = 2 |
S&E (2017) / rejoinder code |
| Fixed scalar \(\gamma\) crossed with the \(\delta\) grid | opt-in (gamma = 2.5) |
Schnitzer et al. (2025) |
| Tuning criterion: wAMD with \(\lvert\hat\beta_j\rvert\) weights | the only criterion | S&E (2017) |
| Estimand for the wAMD weights | "ATE" |
S&E (2017) (note: psave() defaults to ATT) |
| PS from the penalized fit | refit = FALSE |
S&E (2017) / rejoinder code |
| Post-selection unpenalized refit | opt-in (refit = TRUE) |
Schnitzer et al. (2025), adapted from a longitudinal setting |
| GOAL \(\ell_2\) term, augmentation + \((1+\lambda_2)\) rescale (all coefficients incl. intercept) | method = "goal" |
Baldé et al. (2023); Baldé (2025) supplement |
GOAL lambda2 grid |
c(0, 10^c(-2, -1.5, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 1)) |
Author’s published grid (Baldé 2025 supplement, “taken from Zou and Hastie (2005)”; verified against the official code 2026-07-02) |
| GOAL penalty base \((n+q)^\delta\) on augmented grid points | fixed | Baldé (2025) supplement (n.q = n + q) |
Standardize once, glmnet(standardize = FALSE) |
fixed protocol | Rejoinder code convention |
| Exact penalty-scale correction \(s = \overline v_{\text{active}} \lambda_n / n\) | always on | Rejoinder code (KKT-verified) |
| Clipping constants | clip = c(0.01, 0.99) |
Package choice, S&E lineage (suite constant shared with
psAve; same constants in the legacy code). Baldé’s
reference runs unclipped — use
clip = c(1e-12, 1 - 1e-12) to reproduce it |
| Ties: first minimum, \(\lambda_2\) then \(\delta\) ascending | fixed | Package choice (consistent with Schnitzer et al.’s smallest-lambda tie note and psAve’s tolerant first-minimum) |
NA anywhere in used variables is an error |
fixed | Package/suite policy (never silent dropping) |
lqa (Ulbricht 2010) — the package
S&E’s original implementation was built on — is archived and removed
from CRAN; oalasso does not depend on or reference it as a
dependency. The glmnet route with the exact-scale
correction replaces it.jmiahjones/OAL-rejoinder-2022 (GitHub)
— the rejoinder’s companion code — is the reference implementation: its
oal_funs.R recipe (standardization, pairing, \(s\)-correction, GOAL augmentation) and its
KKT test are what oalasso reproduces and tests
against.refit = TRUE variant and the fixed-\(\gamma\) mode are adapted. Their SE-scaled
balance criterion is not implemented in v1 (its validation is
longitudinal-only, and its bootstrap would break the package’s
determinism).QianGao-SXMU/GOAL (Gao et
al.) is a different method — a generalized propensity score for
continuous exposures — that happens to share the acronym.
method = "goal" here is Baldé’s.ctmle (CRAN) implements collaborative
TMLE (van der Laan & Gruber 2010, doi:10.2202/1557-4679.1181; scalable versions Ju et
al. 2019, doi:10.1177/0962280217729845) — the target-driven,
learner-agnostic parent framework of outcome-aware nuisance
selection.hal9001 (CRAN) implements the highly
adaptive lasso, on which the outcome-highly-adaptive lasso (oHAL; Ju,
Benkeser & van der Laan 2020, doi:10.1111/biom.13121) is built — see the next
section.tom-beer/Outcome-Adaptive-LASSO (GitHub).lqa-based referencesDuring development we ran oalasso head-to-head against
the official GOAL/OAL reference code from the supplement of
Baldé (2025) on identical data, with the archived lqa
engine installed in an isolated library (the full report ships with the
package sources under verification/). The outcome-model
coefficients agree to machine precision, but the selected \((\lambda, \gamma)\) can differ. The root
cause is reproducible and instructive: lqa.default() is
called with its default standardize = TRUE, which rescales
the (already SD-standardized) columns to unit Euclidean norm
internally. On SD-standardized data that norm is exactly \(\sqrt{n-1}\), so for plain
OAL the reference code effectively solves the published
objective with the L1 constant inflated to \(\lambda_n\sqrt{n-1}\). Bridging
glmnet with that factor reproduces the lqa OAL
grid points and the reference’s exact OAL selection (the residual gap
traces to lqa’s loose conv.eps = 1e-3;
tightening it collapses the difference to \(\sim 10^{-6}\)). For GOAL with
\(\lambda_2 > 0\) the story is
messier and a single scale factor does not suffice:
lqa standardizes the augmented design (realized
column norms \(\sqrt{n - 1 + \lambda_2(1 -
1/(n+q))}\), with nonzero augmented-row column means interacting
with the intercept), and its modified PIRLS treats the augmentation rows
as weight-1/response-0 ridge terms inside each IRLS step —
structurally different from the binomial pseudo-row objective a
glmnet augmentation solves. GOAL grid points therefore
agree only approximately with the reference solver even after bridging;
selections on our test data agreed at \(\lambda_2 = 0\) and differed in PS by \(\sim 4\times 10^{-2}\) at the selected
point overall (full details in the equivalence report).
oalasso deliberately solves the objective as printed
in the papers — that is what the KKT acceptance tests pin down —
rather than reproducing the historical solver’s scale. Because \(\lambda_n = n^{\delta}\) and \(\sqrt{n-1}\approx n^{1/2}\), users who want
selections on the historical lqa scale can simply shift the
exponent grid:
oal(A ~ ., data = dat, outcome = ~ y,
lambda = c(-10, -5, -2, -1, -0.75, -0.5, -0.25, 0.25, 0.49) + 0.5)In practice the wAMD criterion re-selects the tuning adaptively, so
downstream propensity scores are typically close; but simulation
comparisons against numbers published with lqa-based code
should keep this scale difference in mind.
Everything above — the penalty weights, the wAMD weights, and both
published oracle theorems — rests on the linear outcome model.
Under misspecification, \(\hat\beta\)
converges to best-linear-projection coefficients, and the set it zeroes
need not be the outcome-unrelated set: a confounder whose outcome effect
is purely nonlinear can have a zero projection coefficient and be
wrongly dropped, while a null covariate correlated with a strong
predictor can be wrongly kept. Selection consistency for the confounder
set is therefore not guaranteed under misspecification, and to our
knowledge no dedicated published study isolates OAL’s selection behavior
in that regime (a literature reading, not a theorem); the rejoinder
explicitly used a correctly specified linear model and pointed to
flexible cross-fitted methods as the way forward. Hand-coding basis
terms (splines, interactions) into formula and
outcome is the within-scope remedy, since both models then
share the expanded columns.
For genuinely nonparametric outcome-adaptive selection, use the
published tools: the outcome-highly-adaptive lasso (Ju,
Benkeser & van der Laan 2020, doi:10.1111/biom.13121), which selects over HAL’s
spline-interaction basis via hal9001 and pairs with TMLE,
and the C-TMLE framework (ctmle; van der
Laan & Gruber 2010), which selects nuisance models by their value
for the target parameter with any learner. Those methods carry their own
theory; oalasso does not attempt to reimplement them.
outcome.coef hookoal()’s single extension seam is the
outcome.coef argument: a named numeric vector \(b_j\) over the standardized model-matrix
columns that replaces the internal outcome model in both the
penalty (\(\hat w_j =
|b_j|^{-\gamma}\)) and the wAMD weights. The values are
interpreted on the standardized-\(X^s\)** scale** — adaptive weights are
scale-dependent, so coefficients from an unstandardized fit will not
mean what you think. Zeros are allowed here (unlike in the internal
model): \(|0|^{-\gamma} = \infty\)
becomes a hard exclude in glmnet. Using the hook switches
the printed provenance to an explicit experimental label — screening use
only.
The reference construction (reproducing the internal weights by hand)
and the ridge variant for \(p\) close
to \(n\) (the natural analogue of the
author’s legacy shrink = TRUE branch):
## Reference construction: what oal() computes internally
f <- A ~ x1 + x2 + x3 + x4 + x5 + x6
X <- model.matrix(f, sim)[, -1, drop = FALSE]
Xs <- scale(X, TRUE, TRUE) # the package's fixed protocol
b <- coef(lm(sim$Y ~ sim$A + Xs))[-(1:2)] # drop intercept and treatment
names(b) <- colnames(Xs)
fit.hook <- oal(f, data = sim, outcome.coef = b) # same selections as outcome = ~ Y
# (provenance label differs)
## Ridge outcome model for p close to n
set.seed(1) # cv.glmnet folds are random
cv <- glmnet::cv.glmnet(x = cbind(A = sim$A, Xs), y = sim$Y, alpha = 0,
penalty.factor = c(0, rep(1, ncol(Xs))))
b <- as.numeric(coef(cv, s = "lambda.min"))[-(1:2)]
names(b) <- colnames(Xs)
fit.ridge <- oal(f, data = sim, outcome.coef = b)A tempting idea is to make the weights nonlinear-aware by replacing \(|\hat\beta_j|\) with a tree-ensemble variable importance. The package deliberately does not export this: no published, peer-reviewed method plugs tree importance into an OAL-style adaptive penalty, and the theory does not transfer. The following framing, reproduced verbatim from the package’s theory review, is the ceiling of what may be claimed for such a construction:
- Claim it as a screening / dimension-reduction procedure targeting outcome-relevant covariates (X_C ∪ X_P), not an oracle selector.
- State the target guarantee as a sure screening property (Fan & Lv 2008) — no false negatives w.p.→1 — with a second-stage refit for false positives.
- Do not claim the OAL/GOAL oracle property (√n asymptotic normality with oracle variance): tree importance is not a √n-consistent signed coefficient, so the Zou/Shortreed–Ertefaie proof does not transfer.
- Do not claim λ_n rate conditions (λ_n/√n→0, λ_n n^{γ/2−1}→∞) confer selection consistency here — those depend on |β̃ⱼ|=O_p(n^{−1/2}) under the null, which importance measures do not satisfy.
- Explicitly warn that impurity/permutation importance is biased with mixed-scale (Strobl 2007) and correlated (Strobl 2008) covariates, which can let an instrument masquerade as relevant.
- Default to conditional permutation importance and unbiased/honest trees to mitigate those biases; document this as a requirement, not an option.
- Offer model-X knockoffs (Candès et al. 2018) as the finite-sample FDR-controlled confounder-selection mode, flagging its dependence on an accurate covariate model.
- Keep the instrument-exclusion rationale (Brookhart 2006; Myers 2011; Pearl 2011): beneficial under possible unmeasured confounding and IPW estimation; note nonparametric outcome modeling does not remove Z-bias or positivity concerns.
- Recommend pairing selection with a doubly robust estimator (AIPW/TMLE) + cross-fitting, so misspecification/selection error degrades gracefully rather than breaking a fragile oracle claim.
- Summarize honestly: “empirically competitive, screening-consistent, FDR-controllable confounder selection for nonlinear settings” — not “oracle-consistent causal variable selection.”
Under that framing, the recipe below feeds conditional
permutation importance from a conditional-inference forest
(party::cforest with cforest_unbiased() — the
Strobl-recommended configuration; see points 5–6 above) into the
outcome.coef hook. It is intentionally left unevaluated:
run it only with the framing above in mind, and pair the result with a
doubly robust downstream estimator (point 9).
## EXPERIMENTAL — screening use only, per the ten points above. Requires `party`.
library(party)
f <- A ~ x1 + x2 + x3 + x4 + x5 + x6
X <- model.matrix(f, sim)[, -1, drop = FALSE]
Xs <- scale(X, TRUE, TRUE) # same standardization as oal()
df <- data.frame(Y = sim$Y, Xs, check.names = FALSE)
set.seed(20260702) # forests are stochastic
cf <- cforest(reformulate(colnames(Xs), response = "Y"), data = df,
controls = cforest_unbiased(ntree = 1000,
mtry = ceiling(sqrt(ncol(Xs)))))
vi <- varimp(cf, conditional = TRUE) # conditional permutation importance
vi <- pmax(vi[colnames(Xs)], 0) # negative importance -> 0 = hard exclusion
fit.imp <- oal(f, data = sim, outcome.coef = vi)
fit.imp # provenance line flags the user-supplied weights; screening use onlyNotes on the recipe. The names of the vector must exactly cover the
standardized model-matrix columns (factor covariates must be expanded
via model.matrix first, as above). Importance is
nonnegative, which is fine: only \(|b_j|\) enters the penalty and the wAMD.
Zeroed importances become hard exclusions. And the choice of importance
measure is load-bearing, not cosmetic — plain impurity or marginal
permutation importance can hand an instrument a small penalty precisely
when it is correlated with a confounder, defeating the method’s purpose
(Strobl et al. 2007, 2008).
Baldé, I. (2025). The oracle property of the generalized outcome-adaptive lasso. Statistics & Probability Letters, 221, 110379. doi:10.1016/j.spl.2025.110379
Baldé, I., Yang, Y. A., & Lefebvre, G. (2023). Reader reaction to “Outcome-adaptive lasso: Variable selection for causal inference” by Shortreed and Ertefaie (2017). Biometrics, 79(1), 514–520. doi:10.1111/biom.13683
Brookhart, M. A., Schneeweiss, S., Rothman, K. J., Glynn, R. J., Avorn, J., & Stürmer, T. (2006). Variable selection for propensity score models. American Journal of Epidemiology, 163(12), 1149–1156. doi:10.1093/aje/kwj149
Candès, E., Fan, Y., Janson, L., & Lv, J. (2018). Panning for gold: ‘model-X’ knockoffs for high dimensional controlled variable selection. Journal of the Royal Statistical Society: Series B, 80(3), 551–577. doi:10.1111/rssb.12265
Fan, J., & Lv, J. (2008). Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B, 70(5), 849–911. doi:10.1111/j.1467-9868.2008.00674.x
Jones, J., Ertefaie, A., & Shortreed, S. M. (2023). Rejoinder to “Reader reaction to ‘Outcome-adaptive lasso: Variable selection for causal inference’ by Shortreed and Ertefaie (2017)”. Biometrics, 79(1), 521–525. doi:10.1111/biom.13681
Ju, C., Benkeser, D., & van der Laan, M. J. (2020). Robust inference on the average treatment effect using the outcome highly adaptive lasso. Biometrics, 76(1), 109–118. doi:10.1111/biom.13121
Ju, C., Gruber, S., Lendle, S. D., Chambaz, A., Franklin, J. M., Wyss, R., Schneeweiss, S., & van der Laan, M. J. (2019). Scalable collaborative targeted learning for high-dimensional data. Statistical Methods in Medical Research, 28(2), 532–554. doi:10.1177/0962280217729845
Myers, J. A., Rassen, J. A., Gagne, J. J., Huybrechts, K. F., Schneeweiss, S., Rothman, K. J., Joffe, M. M., & Glynn, R. J. (2011). Effects of adjusting for instrumental variables on bias and precision of effect estimates. American Journal of Epidemiology, 174(11), 1213–1222. doi:10.1093/aje/kwr364
Pearl, J. (2011). Invited commentary: Understanding bias amplification. American Journal of Epidemiology, 174(11), 1223–1227. doi:10.1093/aje/kwr352
Schnitzer, M. E., Talbot, D., Liu, Y., Berger, C., Wang, G., O’Loughlin, J., Sylvestre, M.-P., & Ertefaie, A. (2025). Adaptive sparsening and smoothing of the treatment model for longitudinal causal inference using outcome-adaptive LASSO and marginal fused LASSO. Statistics in Medicine. doi:10.1002/sim.70316
Shortreed, S. M., & Ertefaie, A. (2017). Outcome-adaptive lasso: Variable selection for causal inference. Biometrics, 73(4), 1111–1122. doi:10.1111/biom.12679
Strobl, C., Boulesteix, A.-L., Zeileis, A., & Hothorn, T. (2007). Bias in random forest variable importance measures: Illustrations, sources and a solution. BMC Bioinformatics, 8, 25. doi:10.1186/1471-2105-8-25
Strobl, C., Boulesteix, A.-L., Kneib, T., Augustin, T., & Zeileis, A. (2008). Conditional variable importance for random forests. BMC Bioinformatics, 9, 307. doi:10.1186/1471-2105-9-307
van der Laan, M. J., & Gruber, S. (2010). Collaborative double robust targeted maximum likelihood estimation. The International Journal of Biostatistics, 6(1), Article 17. doi:10.2202/1557-4679.1181
Zou, H. (2006). The adaptive lasso and its oracle properties. Journal of the American Statistical Association, 101(476), 1418–1429. doi:10.1198/016214506000000735