--- title: "Method details and provenance" author: "Daijiro Kabata" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Method details and provenance} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` 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. ## Notation 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. ## The OAL estimator ### Step 1: the outcome (weight) model 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. ### Step 2: the adaptively penalized treatment model $$ \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. ### Step 3: the tuning grid `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`). ### Step 4: the wAMD selection criterion 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. ### Step 5: the propensity score 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`. ## GOAL: the generalized outcome-adaptive lasso `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)`. ## Standardization protocol (fixed, not an argument) $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`. ## The exact penalty-scale correction (why `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: 1. **Per-observation vs total loss.** glmnet minimizes $-\tfrac{1}{n}\ell_n(\alpha) + \lambda \sum_j \tilde{v}_j |\alpha_j|$ — the *average* log-likelihood — while OAL's $\lambda_n = n^\delta$ grid is defined against the *total* log-likelihood. Ignoring this changes the effective penalty by a factor of $n$, and the $n^\delta$ grid loses its theoretical meaning. 2. **Penalty-factor rescaling.** glmnet internally rescales the supplied penalty factors to sum to `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): ```{r, message = FALSE} 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")] all.equal(fit$path$gamma, 2 * (2 - fit$path$delta + 1)) ``` ## Provenance of every default | 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) | ## Differences from the author's legacy code `oalasso` supersedes the OAL implementations (`fAL`, `fOAL`) in the author's earlier simulation code, which were built on the archived `lqa` package. Four defects of that code are deliberately not reproduced, so outputs *will* differ where those code paths were exercised: 1. **The `fAL` gamma formula.** `fAL` computed `gamma = 2*(gcf - delta) + 1`, which disagrees with its own comment and gives $\gamma$ smaller by 1 than the S&E convention — a transcription bug. Only the correct pairing $\gamma = 2(g_{cf} - \delta + 1)$ (as in `fOAL` and the rejoinder) is implemented; the buggy formula appears nowhere in the package, and a regression test asserts its absence from the grid. 2. **The `varlist` wAMD bug.** One `fOAL` branch passed the *treatment* variable name as the covariate list to its wAMD function, so the criterion was computed over the treatment column alone and the lambda selection in that branch was meaningless. `oal_wamd()` always sums over all $p$ covariate columns. 3. **Global-`n` scoping.** In the legacy code, the `n` in `lambda = n^delta` was never defined inside the functions — it resolved from the calling (global) environment, silently using whatever `n` the surrounding script happened to hold. In `oalasso`, $n$ is always `nrow(model.matrix)`, by construction. 4. **Per-subset `scale()`.** The legacy code standardized fitting and evaluation subsets separately, by different means and SDs — which changes the adaptive weights and the scores. The fixed protocol above (scale once, training statistics) makes this unreproducible by design. ## Relation to other software * **`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. * **Schnitzer-Biostats-Lab/Longitudinal-outcome-adaptive-LASSO** (GitHub; Schnitzer et al. 2025, doi:10.1002/sim.70316) — the longitudinal OAL pipeline from which the `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). * **GOAL disambiguation.** Baldé's GOAL has no public repository; the official code lives in the Wiley supplement of Baldé et al. (2023). The GitHub repository `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. * A Python OAL implementation exists at `tom-beer/Outcome-Adaptive-LASSO` (GitHub). ## A verified scale difference against the historical `lqa`-based references During 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: ```{r, eval = FALSE} 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. ## Beyond linearity: honest notes on nonlinear extensions ### The linear limitation 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. ### The `outcome.coef` hook `oal()`'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): ```{r, eval = FALSE} ## 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) ``` ### Experimental: importance-based weights (a recipe, not a feature) 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: > 1. Claim it as a **screening / dimension-reduction** procedure targeting outcome-relevant covariates (X_C ∪ X_P), not an oracle selector. > 2. 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. > 3. 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. > 4. 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. > 5. 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. > 6. **Default to conditional permutation importance and unbiased/honest trees** to mitigate those biases; document this as a requirement, not an option. > 7. 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. > 8. 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. > 9. 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. > 10. 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). ```{r, eval = FALSE} ## 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 only ``` Notes 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). ## References 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