Selecting Covariates for PATED Using Historical Data

Statistic-specific variance reduction, artificial grouping, and sparse influence-score projection

Abstract

Historical data can inform both the choice of covariates for a future randomized trial and the variance reduction expected from their adjustment. We develop this planning approach for prognostic variables assisted treatment effect detection (PATED), as implemented in the R package multipleOutcomes. PATED subtracts a covariance-weighted linear combination of baseline-variable contrasts from a primary marginal statistic. Its population variance reduction defines a statistic-specific \(R^2\), which need not equal a conventional regression goodness-of-fit measure. We construct the corresponding historical-data target using a common independent artificial treatment label and then eliminate label generation analytically. This yields a projection of the primary statistic’s influence score onto candidate covariate scores, allowing different variable types without specifying a joint parametric response model.

Adaptive penalization selects nonzero additional projection contributions. We establish a fixed-dimensional, pointwise oracle property under a root-\(n\) expansion of the estimated score covariance, including first-stage estimation effects. New historical-data simulations use \(p=20\), three to five signals, development samples of 5,000–20,000, and population \(R^2\) values of 8%–15%, with separate validation samples. They distinguish exact selection from retained variance reduction and apparent from externally evaluated performance. Existing diagnostic simulations are retained in an appendix. The target is null-reference variance reduction in a specified population and observation regime; its use for a future trial requires an explicit transportability argument.

Keywords: covariate adjustment; historical data; PATED; variance reduction; artificial grouping; influence function; adaptive lasso.

1. Introduction

1.1 Covariate adjustment and the planning problem

Randomization separates treatment assignment from patients’ baseline characteristics, but it does not remove the sampling variability associated with prognostic differences between participants. Baseline information can therefore improve the precision of a treatment comparison even when the unadjusted analysis is already valid. The FDA’s May 2023 guidance recommends considering prognostic baseline covariates in randomized, parallel-group trials, emphasizes prespecification of the analysis, and identifies previous studies as a source of information about useful covariates [11]. Its purpose is to improve randomized comparisons, rather than to resolve confounding in observational treatment comparisons. Historical data enter the present article as a resource for planning an eventual randomized analysis.

A useful distinction is between requiring a response model to be correct and using a working model to improve estimation. Ye, Shao, Yi and Zhao [12] organize their analysis around efficiency gain, applicability to the randomization design, and appropriate standard errors. Their results explain why a valid adjusted estimator need not automatically improve precision and why the particular adjustment matters. Reluga, Ye and Zhao [15] likewise compare regression adjustments without assuming a correct linear conditional-mean model. Thus, freedom from a correctly specified parametric outcome model is already an established feature of parts of the covariate-adjustment literature; it is not a novelty claim of the present work.

This perspective extends beyond linear working models. Bannick and colleagues [13] study augmented inverse probability weighted adjustment using generalized linear, nonparametric, and machine-learning working models, with explicit conditions for validity and efficiency under covariate-adaptive randomization. For survival endpoints, Ye, Shao and Yi [14] develop adjustment of the log-rank statistic itself. These developments suggest a principle that guides our construction: evaluate prognostic information in relation to the statistic that will actually be used for the treatment comparison. Predicting an outcome well under some convenient model and reducing the variance of a particular trial statistic are related tasks, but they are not interchangeable definitions.

Here we focus on PATED, whose name in multipleOutcomes is prognostic variables assisted treatment effect detection [9]. The package jointly analyzes a primary endpoint and baseline prognostic variables, each through its own marginal analysis against the same treatment label. Because genuinely baseline variables have no population treatment contrast under randomization, their estimated contrasts provide mean-zero adjustment terms. The primary estimate is corrected by a linear combination of those terms, with coefficients obtained from their joint covariance. This construction makes variance reduction a directly defined quantity even when the primary and auxiliary analyses use different types of statistics. Throughout this article we discuss estimator or statistic variance, rather than translating the results into power.

Two questions then arise before a trial is conducted:

1. Selection: among the available candidate variables, which should be retained as PATED covariates because they make an additional contribution to the planned adjustment?
2. Magnitude: for the selected set, how much of the primary statistic’s variance can that adjustment remove?

A historical trial cohort, a single-arm study, or a suitably defined real-world cohort can help answer these questions. Claims data may provide many candidate diagnoses, utilization measures, prior-treatment characteristics, and event histories. The first question is consequently a variable-selection problem. The second requires a performance measure tied to PATED, rather than a generic ranking of prognostic models. Li, Li and Luedtke [16] already provide a framework for estimating the efficiency gains of future covariate-adjusted analyses from external data. Our planning motivation is directly related to theirs; our specific object is the sparse projection induced by the marginal statistics available to PATED.

1.2 Why the variance target must match the analysis

It would be inaccurate to say that logistic or Cox regression has no \(R^2\). Generalized coefficients of determination and survival explained-variation measures have long been proposed [17,18]. Their definitions differ: they can summarize likelihood improvement, prediction error, or explained variation in a survival process. None is automatically the proportion by which a particular PATED treatment statistic’s variance is reduced. We therefore define \(R^2\) through the joint covariance of that primary statistic and its auxiliary contrasts. Its interpretation follows from the adjustment formula, without borrowing an interpretation from another regression model.

A second consideration concerns how different variables enter an analysis. Standard regression readily accepts numeric covariates, including a fully observed prior time-to-progression. The difficulty is more specific when the prior event time is censored: entering the observed minimum of event and censoring times as if it were the true event time changes the information being represented. For example, first-line time-to-progression may be prognostic for an endpoint in a later-line trial, yet some prior event times may be unresolved or administratively censored. An appropriate time-to-event marginal statistic can retain the observed time and event indicator together. In the PATED construction this statistic can serve as an auxiliary analysis, while the primary endpoint uses its own scale and observation model.

This flexibility does not make all variable encodings equivalent. It requires a clinically chosen summary, a valid observed-data score, and appropriate assumptions about censoring or missingness. Only information available before the future randomization can be used as a baseline covariate. The benefit is a common covariance-based adjustment space for heterogeneous summaries, without requiring all original measurements to enter one conditional response family.

1.3 From artificial labels to a historical-data analysis

A historical cohort need not contain the future treatment comparison. Instead, assign each person an independent artificial binary label and use that same label for every marginal analysis. The two artificial arms have identical population distributions, so all their population contrasts are zero. Nevertheless, the estimated contrasts fluctuate together: a split producing an unusual primary-endpoint statistic can also produce unusual statistics for related baseline variables. Their cross-covariances quantify precisely the linear adjustment opportunity that PATED uses.

We formalize that opportunity using influence scores, the subject-level contributions to an estimator’s first-order sampling fluctuation. The primary score is projected on the candidate covariate scores, and an adaptive lasso penalty selects additional contributions. Repeated artificial grouping is not needed: under the null factorization proved below, the required covariance is a scaled covariance of label-free scores computed from the historical observations. The resulting \(R^2\) describes variance reduction at this null reference, where a reduction remains meaningful despite a zero treatment contrast.

The historical target is not automatically a forecast for every future trial. The endpoint definition, baseline measurement rules, eligibility criteria, follow-up horizon, censoring mechanism, and relevant joint distribution must be comparable or explicitly transported. Our formal construction uses independent artificial allocation; it does not establish the validity of a trial analysis under every covariate-adaptive allocation scheme considered in [12–14]. Nor does a scalar null representation generally determine the covariance under a large nonnull treatment effect. These are boundaries of the planning interpretation, not reasons to replace the statistic-specific target by a generic regression score.

The contribution developed here is the connection between PATED’s covariance adjustment, historical-data selection, analytic removal of artificial randomization, and an oracle result that accounts for estimated scores. Selection and penalized control variates have substantial precedents [1–5], and external-data efficiency assessment is already studied [16]. We make no claim that the projection identity or the broad historical-data strategy is new.

2. PATED, its variance \(R^2\), and compatible marginal statistics

2.1 The adjustment is a linear combination

For a scalar primary statistic, let \(\hat{\theta}_Y\) estimate its target \(\theta_Y\) and let \(\hat{\delta}_X\) collect the \(p\) baseline-variable treatment contrasts, whose population vector is zero. Suppose

\[ \sqrt n \begin{pmatrix} \hat{\delta}_X\\ \hat{\theta}_Y-\theta_Y \end{pmatrix} \ \xrightarrow{d}\ N\left\{0, \begin{pmatrix} \Sigma_{XX}&\Sigma_{XY}\\ \Sigma_{YX}&\Sigma_{YY} \end{pmatrix}\right\}. \tag{P1} \]

The central requirement is joint asymptotic normality of the chosen marginal statistics, with a consistently estimable joint covariance. Separate normal limits for each marginal statistic do not by themselves establish (P1). A joint asymptotic linear expansion in aligned subject-level contributions, followed by a multivariate central limit theorem, is one sufficient route. The raw clinical measurements need not be jointly normal. We restrict the implementation discussion to analytic score-based adapters; KM and quantile adapters requiring the current package’s bootstrap route are excluded.

With \(\Sigma_{XX}\) positive definite, minimizing the first-order variance of \(\hat{\theta}_Y-c^\top\hat{\delta}_X\) gives

\[ c_*=\Sigma_{XX}^{-1}\Sigma_{XY}, \qquad \hat{\theta}_{\mathrm{PATED}} =\hat{\theta}_Y-\hat{c}_*^\top\hat{\delta}_X, \qquad V_{\mathrm{PATED}} =\Sigma_{YY}-\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY}. \tag{P2} \]

Here \(V_{\mathrm{PATED}}\) is the variance of the limiting \(\sqrt n\) fluctuation; the estimator variance is approximately \(V_{\mathrm{PATED}}/n\). A consistent \(\hat{c}_*\) has the same first-order variance, since \((\hat{c}_*-c_*)^\top\hat{\delta}_X=o_p(n^{-1/2})\). The code in pated.R uses the opposite coefficient sign, opt_c = -solve(mcov22) %*% mcov21, and then adds its linear combination. This is exactly (P2) on the working scale. The covariance matrix returned by vcov() is on the estimator scale, so its common \(1/n\) factor cancels in the solve.

Define

\[ R_{\mathrm{PATED}}^2 =\frac{\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY}}{\Sigma_{YY}}, \qquad \frac{V_{\mathrm{PATED}}}{\Sigma_{YY}}=1-R_{\mathrm{PATED}}^2. \tag{P3} \]

Thus \(R_{\mathrm{PATED}}^2=0.12\) means a 12% first-order variance reduction for the optimal linear adjustment in the stated population and statistic system. The standard-error ratio is \(\sqrt{0.88}\); the package’s scalar Rel. Eff. attribute is the inverse variance ratio \(1/(1-R_{\mathrm{PATED}}^2)\). Positive semidefiniteness gives \(0\le R_{\mathrm{PATED}}^2\le1\). A singular auxiliary block must be reduced or treated as a block-identification problem before applying this formula.

This is a projection of sampling fluctuations, not a regression of raw endpoints on treatment labels. It explains variance of the primary treatment statistic, not necessarily variance of the original clinical outcome. The zero population auxiliary contrast is essential: raw baseline means or arm-specific survival probabilities cannot simply be inserted in \(\hat{\delta}_X\).

2.2 From the PATED estimator to influence functions

The route to the previous projection formulation is an expansion of the PATED estimator itself. Suppose the marginal estimators admit the joint asymptotically linear representation

\[ \hat{\theta}_Y-\theta_Y =\frac1n\sum_{i=1}^n\phi_Y(O_i)+o_p(n^{-1/2}), \qquad \hat{\delta}_X =\frac1n\sum_{i=1}^n\phi_X(O_i)+o_p(n^{-1/2}). \tag{P4a} \]

The centered, finite-variance vector \((\phi_X^\top,\phi_Y)^\top\) is the complete influence function of the stacked marginal estimates. For consistent adjustment coefficients, (P2) consequently has influence function

\[ \phi_{\mathrm{PATED}}(O) =\phi_Y(O)-c_*^\top\phi_X(O). \tag{P4b} \]

Minimizing its variance is exactly the projection \(\min_c E\{(\phi_Y-c^\top\phi_X)^2\}\). This is the connection between the PATED estimator and the influence-function formulation below, rather than a separate justification based on a response regression model. Joint asymptotic linearity is a sufficient condition for joint asymptotic normality; joint normality alone does not imply that an arbitrary estimator has this observation-wise representation.

A likelihood score is not, in general, a special kind of influence function. For a regular finite-dimensional MLE, the estimator’s influence function is \(I^{-1}s\), where \(s\) is its likelihood score and \(I\) its information matrix. For a regular estimating equation it is \(-D^{-1}\psi\). These normalizations, and any nuisance contributions, are part of the influence function. With semiparametric estimators the corresponding efficient-score/information representation requires its own regularity conditions. We use influence function for the complete contribution and label-free influence component for the factor \(h\) below. The shorter phrase influence score elsewhere in this article always means this complete or explicitly factored contribution, not a raw likelihood score.

2.3 Analytic constructions and the package interfaces

The source audit below refers to commit 3bfc7d361d5bca858dce7f49511961f47319fd42 [9]. Write \(Z\) for one fully observed scalar measurement. Under a common independent artificial label \(A\) with probability \(\pi\), the representation sought is

\[ \sqrt n\,\hat{\delta} =\frac{1}{\sqrt n}\sum_{i=1}^n W_i h(O_i)+o_p(1), \qquad W_i=\frac{A_i-\pi}{\pi(1-\pi)}, \qquad E h(O)=0. \tag{P4} \]

The functions \(h\) below can be used for the primary endpoint or for a baseline variable, provided its complete record is available at the required time. Multiplying any one score by a nonzero constant changes its coefficient scale but not the population \(R^2\) or selected population support.

Measurement and contrast	Label-free null score \(h\), or construction	Implemented entry point and qualification
Continuous measurement, difference in means	\(Z-\mu\), \(\mu=E(Z)\)	glm_(Z ~ arm, family="gaussian"); GLMAdapter uses sandwich::estfun() and bread. Normality is not needed to define the mean contrast.
Binary measurement, difference in log odds	\((Z-\mu)/\{\mu(1-\mu)\}\), \(0<\mu<1\)	glm_(Z ~ arm, family="binomial"), with the logit link. Its score is proportional to the risk-difference score at the artificial null.
Count measurement, difference in log means	\((Z-\mu)/\mu\), \(\mu>0\)	glm_(Z ~ arm, family="poisson"), with the log link; use the joint sandwich rather than assuming the mean equals the variance.
Censored event time, treatment-only Cox coefficient	\(M(\infty)/d\), where \(d=E(\Delta)>0\) and \(M\) is defined below	coxph_(Surv(time,event) ~ arm); score residuals and inverse information are combined in CoxphAdapter. The displayed formula is the artificial-null treatment-only case.
Censored event time, log-rank statistic	A nonzero scalar multiple of \(M(\infty)\)	logrank_(); LogRankAdapter returns \(I_n^{-1/2}U_n\), an order-one standardized statistic. Convert to a compatible asymptotically linear scale before combining it with root-\(n\) contrasts.
Hierarchical net benefit or log win ratio	Pairwise Hájek projections in (P6)–(P7) below	netbenefit_(), winratio_(), with nb_continuous(), nb_binary(), and nb_tte().
Repeated or clustered measurements summarized by a coefficient	The desired coordinate or contrast of \(-D^{-1}\psi_i\), with one complete \(\psi_i\) per independent subject/cluster	gee_() and mmrm_(); the adapters use cluster/subject scores. The MMRM path requires the fitted engine to expose score_per_subject. Verify (P4) for the chosen contrast rather than assuming it for every longitudinal model.

The table describes analytic constructions compatible with the joint-normal requirement, not an assertion that every pooled null score is exposed by a public package function. The package already supplies analytic joint-covariance machinery for these adapters; exposing label-free pooled scores for historical-data integration remains an implementation step. A log-rank statistic has no treatment-effect estimation interpretation solely because a return column is named estimate. The implemented standardized log-rank statistic must first be put on the scale used in (P4).

For the survival rows, let \(\widetilde{T}=\min(T,C)\), \(\Delta=I(T\le C)\), \(N(u)=I(\widetilde{T}\le u,\Delta=1)\), and \(Y^{\mathrm{risk}}(u)=I(\widetilde{T}\ge u)\). Under independent right censoring in the population being represented,

\[ M(t)=N(t)-\int_0^tY^{\mathrm{risk}}(u)\,d\Lambda(u), \qquad y(u)=E\{Y^{\mathrm{risk}}(u)\}. \tag{P5} \]

The treatment-only Cox null formula follows because its treatment score has first-order contribution \((A-\pi)M(\infty)\) and its information per subject is \(\pi(1-\pi)d\). The same residual, up to normalization, represents log-rank fluctuations. This derivation does not assume a Cox model for the conditional relationship between the clinical outcome and the candidate covariates.

More generally, an estimating equation \(n^{-1}\sum_i\psi_i(\eta)=0\) gives the complete influence function \(-D^{-1}\psi_i\), where \(D=E(\partial\psi_i/\partial\eta^\top)\). For a contrast \(a^\top\eta\), retain \(-a^\top D^{-1}\psi_i\). In GEE, for example, \(\psi_i=\dot{\mu}_i^\top V_i^{-1}(Z_i-\mu_i)\) before the bread adjustment. A compatible MMRM mean-parameter score has the corresponding subject-level residual form, but the fitted engine’s score and bread should be used. Raw likelihood or estimating-equation scores alone generally omit the outer factors of the sandwich. If these complete subject-level contributions have finite second moments and jointly negligible remainders, their stacked sum gives joint asymptotic normality, even when endpoints have different types. The package aligns subjects by identifier when constructing cross-model covariance blocks [9]; rows from different endpoints must not be treated as independent observations.

2.4 Win-ratio-type endpoints and covariates

The implemented hierarchical comparison examines prespecified endpoint components in order, returning a win, loss, or tie at the first decisive comparison. The components can be continuous, binary, or time-to-event [9]. Let \(w(o,o')\) indicate that the first record wins over the second, let \(\ell(o,o')=w(o',o)\), and set \(k=w-\ell\). Endpoint priorities, directions, margins, and censoring rules are part of the statistic’s definition.

For treatment and control records with distributions \(P_1\) and \(P_0\), define \(p_W=Ew(O_1,O_0)\), \(p_L=E\ell(O_1,O_0)\), net benefit \(\theta_{\mathrm{NB}}=p_W-p_L\), and log win ratio \(\theta_{\mathrm{WR}}=\log(p_W/p_L)\). Write

\[ b_1(o)= \frac{E_0w(o,O_0)-p_W}{p_W} -\frac{E_0\ell(o,O_0)-p_L}{p_L}, \qquad b_0(o)= \frac{E_1w(O_1,o)-p_W}{p_W} -\frac{E_1\ell(O_1,o)-p_L}{p_L}. \tag{P6} \]

The full log-win-ratio influence function is \(A b_1(O)/\pi+(1-A)b_0(O)/(1-\pi)\). These are the row- and column-average pairwise contributions used by WinRatioAdapter, followed by the log-ratio delta method. Both win and loss probabilities must be positive. With fixed endpoint rules and a nondegenerate Hájek projection, the bounded pairwise kernels have a joint U-statistic central limit theorem; the smooth log-ratio map then gives joint asymptotic normality with the other regular marginal estimators. For net benefit, replace each difference of ratios by the centered conditional expectation of \(k\).

Under artificial grouping, \(P_1=P_0=P\), \(p_W=p_L=q>0\), and \(a(o)=E_P k(o,O')\) has mean zero. Antisymmetry gives

\[ h_{\mathrm{NB}}(o)=a(o), \qquad h_{\mathrm{WR}}(o)=a(o)/q, \qquad \mathrm{IF}_{\mathrm{WR}}(A,O)=W h_{\mathrm{WR}}(O). \tag{P7} \]

Consequently, the two statistics have identical population projection \(R^2\) at this null when they use the same comparison kernel. They need not have identical nonnull targets or finite-sample penalized fits. For one continuous endpoint, no ties, zero margin, and larger values preferred, \(a(z)=2F_Z(z)-1\) and \(q=1/2\). This special case connects the construction to centered rank scores without requiring rank correlation as its general definition.

A pooled empirical implementation estimates \(a(O_i)\) by averaging \(k(O_i,O_j)\) over \(j\ne i\), and estimates \(q\) by the corresponding win fraction. These are estimated, mutually dependent scores. Their Gram expansion must include this estimation step; Theorem 1 does not justify treating them as known independent scores. For CV, estimate the reference distribution from the training observations and compare a validation record only against that training reference.

The code’s nb_tte(censor_rule="informative") names a pair-comparison rule; it is not a correction for informative censoring. Under the default longer-is-better direction, an event/censoring pair is declared decisive only when the observed ordering establishes that result; unresolved comparisons pass to the next endpoint. Missing components also pass to the next endpoint in the inspected implementation. Thus the operational win statistic depends on the recorded observation process. Artificial labels reproduce that observed-data null; additional identification assumptions are needed to reinterpret it as a latent uncensored comparison.

2.5 The planning quantity and its validation

For a prespecified subset \(S\), replace the covariance blocks in (P3) by those for \(S\). This gives the optimal variance reduction \(R^2(S)\). Variable selection asks for a small useful subset; reporting only the full-candidate \(R^2\) does not answer that question.

For any fitted, frozen coefficient vector \(c\), its actually achieved population variance reduction is

\[ \mathcal{V}(c) =\frac{2c^\top g_{XY}-c^\top G_{XX}c}{G_{YY}} =R^2-\frac{(c-c_0)^\top G_{XX}(c-c_0)}{G_{YY}}, \tag{P8} \]

where \(G\) and \(c_0\) are defined below. Only the optimal coefficient gives \(R^2\) exactly. A noisy or over-shrunk coefficient can yield less reduction, and \(\mathcal{V}(c)\) can be negative. Apparent in-sample \(R^2\) is therefore insufficient to evaluate a selected adjustment.

Use an independent historical validation sample, or an outer validation scheme enclosing selection, to evaluate \(1-\widehat{\mathrm{Var}}_{\mathrm{val}}(h_Y-c^\top h_X)/ \widehat{\mathrm{Var}}_{\mathrm{val}}(h_Y)\). Keep the fitted score definitions and coefficients fixed when evaluating them. This ratio estimates achieved reduction for that fitted rule; it is not an exactly unbiased estimator of the population ratio. A refit on selected coordinates can be evaluated separately from the penalized coefficients. Outcome records in the validation cohort may be used to compute the validation scores, but not to tune the fitted rule.

For trial planning, the usual output is a selected set, its statistic definitions, and a validated estimate of variance reduction. The eventual trial can re-estimate PATED’s covariance coefficients on this prespecified set. If instead historical coefficients are frozen for use in the trial, (P8) must be evaluated under the trial covariance: population shifts can change its sign. Our historical-data oracle result concerns \(c_0\) and selection in the historical population; it is not a proof of transportability.

3. PATED results

4. Population construction

4.1 Marginal statistics and common artificial groups

Let \(O_1,\ldots,O_n\) be independent observations from a population distribution \(P\). Each candidate variable \(X_j\), \(j=1,\ldots,p\), and the response \(Y\) is represented by a prespecified scalar null-contrast score. For differences of regular marginal functionals \(T_j(P)\) and \(T_Y(P)\) these are their centered influence functions; for pairwise comparisons they are the arm-directional projections in (P7). Write

\[ h_j(O;P),\qquad h_Y(O;P),\qquad E_P h_j=E_P h_Y=0. \tag{1} \]

The notation permits observed-data records, such as an event time together with its censoring indicator. The score must be derived for that observation scheme. For fully observed marginal variables, it can depend only on that variable.

Generate the same artificial label \(A_i\sim\operatorname{Bernoulli}(\pi)\) for all variable-specific analyses of subject \(i\), independently of \(O_i\), with fixed \(0<\pi<1\). Let \(\hat{\beta}_j\) be the chosen two-group contrast; differences between groupwise estimators of \(T_j\) are one example. Assume the estimators have compatible, jointly valid asymptotic linear expansions. Because the group distributions coincide, all true contrasts are zero, and

\[ \sqrt n\,\hat{\beta}_j =\frac1{\sqrt n}\sum_{i=1}^n \frac{A_i-\pi}{\pi(1-\pi)}h_j(O_i;P)+o_p(1). \tag{2} \]

This expansion follows by subtracting the two groupwise expansions and using \(n_1/n\to_p\pi\). It also covers other regular contrast estimators when their null influence functions factor in the same way. It is a condition to verify for an adapter, not a guarantee for arbitrary model output. For win statistics, \(h_j\) is an arm-directional projection, not the influence function of the identically zero diagonal functional \(T(P,P)\); conflating these would incorrectly give a zero score. In the historical-data setting, (7) below equals the null-reference PATED quantity in (P3).

Stack \(h=(h_X^\top,h_Y)^\top\), where \(h_X=(h_1,\ldots,h_p)^\top\), and define

\[ G=E_P(hh^\top) =\begin{pmatrix}G_{XX}&g_{XY}\\g_{XY}^\top&G_{YY}\end{pmatrix}. \tag{3} \]

Why (3) already includes the sandwich adjustment. Throughout this article, an “influence score” means the complete influence function \(h\), not an unadjusted likelihood score or estimating function. For a regular estimator solving \(n^{-1}\sum_i\psi(O_i;\hat{\theta})=0\), write

\[ D=E\{\partial\psi(O;\theta_0)/\partial\theta^\top\},\qquad B=E\{\psi(O;\theta_0)\psi(O;\theta_0)^\top\}. \]

Its influence function is \(h=-D^{-1}\psi\), so

\[ G=E(hh^\top)=D^{-1}BD^{-\top}. \]

Thus the two outer factors of the sandwich are absorbed into \(h\). No information-matrix equality is assumed. For several marginal estimators, stack their complete influence functions, including any nuisance-estimation contributions; their cross-products supply the off-diagonal blocks. Substituting raw estimating functions for \(h\) would generally omit necessary adjustments.

Proposition 1 (artificial-group covariance). If the joint expansion (2) holds and \(E\|h\|^2<\infty\), the asymptotic covariance of \(\sqrt n(\hat{\beta}_X^\top,\hat{\beta}_Y)^\top\) is

\[ \Sigma=\frac{G}{\pi(1-\pi)}. \tag{4} \]

Proof. Independence of \(A\) and \(O\), together with \(E(A-\pi)^2=\pi(1-\pi)\), gives the covariance in (4); apply the multivariate central limit theorem to (2). \(\square\)

The same subject must receive the same artificial label in every variable-specific analysis. Labels may be independent across subjects. To see why sharing matters, let \(W=(A-\pi)/\{\pi(1-\pi)\}\). A common label gives \(E(W^2h_jh_k)=G_{jk}/\{\pi(1-\pi)\}\). Separate labels that are mutually independent and independent of the data instead give

\[ E(W_jW_kh_jh_k)=E(W_j)E(W_k)E(h_jh_k)=0,\qquad j\ne k. \]

Independent labels therefore remove the desired first-order cross-covariance of the contrasts, even when the original variables are associated. They do not remove association from the underlying data.

4.2 Why one realized vector of contrasts is insufficient

The equation \(\hat{\beta}_Y=c^\top\hat{\beta}_X\) does not identify a meaningful population coefficient from one grouping. For \(p>1\), it is typically one equation with many solutions. Even with \(p=1\), the ratio \(\hat{\beta}_Y/\hat{\beta}_X\) generally converges in distribution to a ratio of correlated normal variables rather than to a constant, since both contrasts have population value zero.

The well-defined target instead minimizes the variance of the unexplained fluctuation:

\[ c_0 =\arg\min_c E\{h_Y-c^\top h_X\}^2 =G_{XX}^{-1}g_{XY} =\Sigma_{XX}^{-1}\Sigma_{XY}. \tag{5} \]

Here and below \(G_{XX}\) is positive definite and \(G_{YY}>0\). Equation (5) makes the target independent of the artificial allocation proportion \(\pi\). If the loss is multiplied by the common covariance scale in (4), its penalty must be multiplied by the same scale to preserve a finite-sample penalized solution.

Let

\[ u=h_Y-c_0^\top h_X,\qquad E(h_Xu)=0, \qquad S_0=\{j:c_{0j}\ne0\}. \tag{6} \]

The orthogonality in (6) holds by projection. Neither \(E(u\mid h_X)=0\) nor homoskedasticity is assumed.

4.3 Association strength and additional contribution

For a subset \(S\) of predictors, write

\[ v(S)=\min_a E(h_Y-a^\top h_S)^2, \qquad R^2=1-\frac{v(\{1,\ldots,p\})}{G_{YY}} =\frac{g_{XY}^\top G_{XX}^{-1}g_{XY}}{G_{YY}}. \tag{7} \]

Define the additional contribution of variable \(j\), conditional on retaining all other candidate variables, as

\[ D_j=\frac{v(-j)-v(\mathrm{all})}{G_{YY}} =\frac{c_{0j}^2G_{j\mid-j}}{G_{YY}}, \qquad G_{j\mid-j}=G_{jj}-G_{j,-j}G_{-j,-j}^{-1}G_{-j,j}. \tag{8} \]

Proposition 2 (interpretation and invariance). Under the stated definiteness conditions, \(0\le R^2\le1\), \(0\le D_j\le1\), and \(D_j=0\) if and only if \(c_{0j}=0\). The support, \(R^2\), and each \(D_j\) are unchanged by separate nonzero rescalings of predictor scores and the response score.

Proof. Residualize \(h_j\) on \(h_{-j}\), obtaining \(r_j\) with variance \(G_{j\mid-j}>0\). Orthogonal projection gives \(v(-j)=v(\mathrm{all})+c_{0j}^2E(r_j^2)\), proving (8) and the bounds. Under \(h_j^*=a_jh_j\) and \(h_Y^*=bh_Y\), the coefficient becomes \(c_{0j}^*=bc_{0j}/a_j\), the residual predictor variance becomes \(a_j^2G_{j\mid-j}\), and response variance becomes \(b^2G_{YY}\). The factors cancel. \(\square\)

The \(D_j\) generally do not add to \(R^2\). They describe deletion of one variable from the full candidate set and therefore depend on that set. They are nonnegative strengths; the sign of \(c_{0j}\) can be reported separately on a prespecified score orientation.

A smooth re-expression of a marginal statistic, \(T_j^*=f_j(T_j)\) with \(f_j'(T_j)\ne0\), rescales its influence function. Proposition 2 therefore gives invariance to such reparameterizations. It does not give invariance to an arbitrary change in the underlying statistic.

Specifically, the influence-function chain rule gives \(h_j^*=f_j'(T_j)h_j\), a nonzero constant rescaling at the fixed population. For example, \(T=E(X)=\mu>0\) and \(T^*=\log\{E(X)\}\) have influence functions \(X-\mu\) and \((X-\mu)/\mu\). Their support, \(R^2\), and \(D_j\) agree, although coefficient magnitudes change. In contrast, replacing \(E(X)\) by \(E(\log X)\), when defined, gives \(\log X-E(\log X)\), generally not a constant multiple of \(X-\mu\). Replacing a mean by a tail probability similarly changes the statistic. The nonzero-derivative condition excludes transformations such as \(T^2\) at \(T=0\), which erase the first-order influence function. These are population invariance statements; they do not guarantee identical finite-sample selections after nonlinear refitting and tuning.

4.4 A distribution-free entry point

Regular model fitting is not mandatory. Choose a fixed measurable feature \(\phi_j\) and take

\[ T_j(P)=E_P\{\phi_j(X_j)\},\qquad h_j=\phi_j(X_j)-E_P\{\phi_j(X_j)\}. \tag{9} \]

This is a nonparametric marginal functional. Bounded features automatically satisfy the moment conditions needed below. Examples include category indicators, prespecified ordinal thresholds, and fixed bounded transformations of continuous variables. Several features for one original variable can be treated as a block, as developed in Section 9.

5. Integrating out artificial randomization

5.1 An exact identity for frozen scores

Fix the observed sample and any associated vectors \(\hat{h}_1,\ldots,\hat{h}_n\). They may have been estimated from the pooled data, but are now held fixed as assignments vary. Assign exactly \(n_1\) observations to group 1 uniformly over all such assignments, with \(n_0=n-n_1\). Let

\[ \Delta_A=\overline{\hat{h}}_{A=1}-\overline{\hat{h}}_{A=0}, \quad S_{\hat{h}}=\frac1{n-1}\sum_i(\hat{h}_i-\overline{\hat{h}}) (\hat{h}_i-\overline{\hat{h}})^\top. \]

Proposition 3 (finite-sample integration).

\[ \operatorname{Cov}_A(\Delta_A\mid O_{1:n},\hat{h}_{1:n}) =\left(\frac1{n_1}+\frac1{n_0}\right)S_{\hat{h}}. \tag{10} \]

Proof. Put \(z_i=\hat{h}_i-\overline{\hat{h}}\). Then \(\Delta_A=n(n_1n_0)^{-1}\sum_i A_iz_i\) and \(\sum_i z_i=0\). For complete randomization, \(\operatorname{Var}(A_i)=n_1n_0/n^2\) and \(\operatorname{Cov}(A_i,A_k)=-n_1n_0/\{n^2(n-1)\}\) for \(i\ne k\). Substitution, using \(\sum_{i\ne k}z_iz_k^\top=-\sum_i z_iz_i^\top\), yields (10). \(\square\)

Thus averaging the covariance over artificial groups is available analytically. One can estimate the association parameter from a centered pooled score covariance. Equation (10) is exact for these frozen-score contrasts. For an estimator that is refitted nonlinearly inside every artificial group, it is an influence-function approximation unless a separate exact identity is available.

Computation without resampling artificial groups. Estimate the influence functions once from the pooled observations and compute \(S_{\hat{h}}\). No artificial labels need to be generated. With \(\pi_n=n_1/n\), the exact conditional covariance of the frozen-score contrast is \(S_{\hat{h}}/\{n\pi_n(1-\pi_n)\}\). The corresponding estimator on the scale of (4) is

\[ \hat{\Sigma}_{\mathrm{int}} =n\left(\frac1{n_1}+\frac1{n_0}\right)S_{\hat{h}} =\frac{S_{\hat{h}}}{\pi_n(1-\pi_n)}. \]

The factor \(n\) matters: \(\Sigma\) is the asymptotic covariance of \(\sqrt n\,\hat{\beta}\), not of \(\hat{\beta}\) itself. Consistency for this population covariance additionally requires consistent estimation of the relevant influence-function second moments. The allocation factor cancels from \(\Sigma_{XX}^{-1}\Sigma_{XY}\), so the projection can be estimated directly from \(S_{\hat{h}}\). Equation (10) assumes uniform allocation with fixed group sizes; Bernoulli allocation gives the same identity conditional on nonempty realized group sizes. This identity does not assert that the finite-sample average of arbitrary covariance outputs from repeatedly refitted nonlinear models equals \(\hat{\Sigma}_{\mathrm{int}}\).

This also distinguishes three objects: a covariance over artificial assignments conditional on the sample; a sampling covariance of fitted contrasts; and the population score covariance \(G\). They are linked by the stated identities and expansions, but are not interchangeable finite-sample estimators.

5.2 How much noise can one grouping add?

For this subsection only, suppose the centered population scores \(h(O_i)\) are known. Generate independent Bernoulli assignments and put

\[ W_i=\frac{A_i-\pi}{\sqrt{\pi(1-\pi)}}. \]

Compare the integrated Gram estimator \(\hat{G}_0=n^{-1}\sum_i h_ih_i^\top\) with \(\hat{G}_A=n^{-1}\sum_i W_i^2h_ih_i^\top\) and their unpenalized projection estimators. Let \(V=G_{XX}^{-1}E(h_Xh_X^\top u^2)G_{XX}^{-1}\).

Proposition 4 (grouping variance factor). The integrated estimator has asymptotic covariance \(V\), whereas the one-assignment estimator has covariance \(\kappa(\pi)V\), where

\[ \kappa(\pi)=E(W^4) =\frac{(1-\pi)^2}{\pi}+\frac{\pi^2}{1-\pi} =1+\frac{(1-2\pi)^2}{\pi(1-\pi)}. \tag{11} \]

If the Gram matrices from \(M\) independent assignments of the same observations are averaged before projection, the factor is

\[ \kappa_M(\pi)=1+\frac{\kappa(\pi)-1}{M}. \tag{12} \]

Proof. Linearizing \(G_{XX}^{-1}g_{XY}\) gives influence \(G_{XX}^{-1}h_Xu\) for the integrated estimator and \(G_{XX}^{-1}W^2h_Xu\) for the randomized estimator. Since \(E(h_Xu)=0\) and \(W\) is independent of the scores, the covariance is multiplied by \(E(W^4)\). For the average \(Q=M^{-1}\sum_{m=1}^M W_m^2\), independence of assignments gives \(E(Q^2)=1+\{E(W^4)-1\}/M\). \(\square\)

At \(\pi=1/2\), \(W_i^2=1\) exactly: the two uncentered Gram matrices coincide. Centering the normalized contrast scores adds only a second-order difference when their mean is \(O_p(n^{-1/2})\). At \(\pi=0.2\), the factor is \(3.25\); with ten averaged assignments it is \(1.225\).

The calculation is specific to this known-score estimator and independent Bernoulli assignments. It is not a universal efficiency ratio for every implementation of a sandwich covariance from refitted models. Estimated-score derivatives can change the comparison. Repeating artificial assignments also does not multiply the number of independent subjects. Averaging already-computed coefficients, rather than averaging Gram matrices first, is a different finite-sample operation.

6. Sparse estimation and the oracle theorem

6.1 Objective and adaptive weights

From the original motivation to a loss function. The motivating question is whether \(c^\top\hat{\beta}_X(A)\) can represent \(\hat{\beta}_Y(A)\) across common artificial groupings. A single grouping generally cannot identify \(c\). We make “represent” precise by minimizing the mean squared representation error across groupings. This is a choice of association target, not a claim to detect every form of dependence.

For the frozen-score contrasts of Proposition 3, each contrast has conditional randomization mean zero. With \(k_n=1/n_1+1/n_0\), the exact identity is

\[ E_A\{(\Delta_{Y,A}-c^\top\Delta_{X,A})^2\mid O_{1:n},\hat{h}_{1:n}\} =k_n\{S_{\hat{h},YY}-2c^\top S_{\hat{h},XY}+c^\top S_{\hat{h},XX}c\}. \]

Divide by the positive allocation factor and by two, convert the sample covariance divisor from \(n-1\) to \(n\) if desired, and discard the term independent of \(c\). The resulting quadratic loss is \(\tfrac12c^\top\hat{G}_{XX}c-c^\top\hat{g}_{XY}\). Adding an adaptive lasso penalty gives (14). Any overall rescaling of the loss requires the same rescaling of \(\lambda_n\) to preserve the solution. For groupwise nonlinear refits, the connection uses expansion (2) and is first-order rather than a general finite-sample identity.

Equivalently, let \(\hat{h}_i^{\,c}=\hat{h}_i-\overline{\hat{h}}\). If \(\hat{G}=n^{-1}\sum_i\hat{h}_i^{\,c}(\hat{h}_i^{\,c})^\top\), the objective below, up to a constant, is

\[ \frac1{2n}\sum_i(\hat{h}_{Yi}^{\,c}-c^\top\hat{h}_{Xi}^{\,c})^2 +\lambda_n\sum_jw_{jn}|c_j|. \]

Thus an imagined regression across many artificial groupings can be implemented directly as a projection of pooled influence functions, without generating those groupings.

Let \(\hat{G}\) estimate (3), with blocks \(\hat{G}_{XX}\) and \(\hat{g}_{XY}\). In the simplest implementation it is the centered empirical covariance of pooled scores, using divisor \(n\). A pilot estimator is a preliminary or initial estimator used to construct the final estimator. Here it is the unpenalized full-model projection estimate of \(c_0\); its role is to set the adaptive penalty weights. It can use the same observations as the final fit and does not require a separate pilot study or sample. Define this initial estimator and its associated weights by

\[ \tilde{c}=\hat{G}_{XX}^{-1}\hat{g}_{XY},\qquad w_{jn}=(|\tilde{c}_j|\vee\epsilon_n)^{-\gamma},\quad \gamma>0. \tag{13} \]

The penalized estimate is

\[ \hat{c}=\arg\min_c\left\{ \frac12c^\top\hat{G}_{XX}c-c^\top\hat{g}_{XY} +\lambda_n\sum_{j=1}^p w_{jn}|c_j|\right\}. \tag{14} \]

Equivalently, its quadratic term is \(\tfrac12(c-\tilde{c})^\top\hat{G}_{XX}(c-\tilde{c})\) up to a constant. The weights must be based on an estimate of \(c_0\), not on the artificial contrasts \(\hat{\beta}_j\), whose population values are all zero.

This has the adaptive-lasso form introduced by Zou [1] and resembles the general quadratic strategy of Wang and Leng [2]. However, \(\hat{G}_{XX}\) here is the covariance of predictor scores. It is generally not the inverse sampling covariance of \(\tilde{c}\). Consequently, a likelihood-efficiency interpretation from another LSA application cannot simply be imported. The following result is proved for this particular projection objective.

6.2 Assumptions

The theorem is for fixed \(p\), which includes a fixed dimension such as 100 as sample size increases. A growing-dimensional sufficient condition is given separately in Section 9.

A1. Sampling and a fixed target. The sampling units are independent and identically distributed. The score system and its population limits define a deterministic \(G\). The true support \(S_0=\operatorname{supp}(c_0)\) is fixed, with \(1\le |S_0|<p\); write \(N_0=S_0^c\). These are population sets, distinct from the random selected set \(\hat{S}=\operatorname{supp}(\hat{c})\).

A2. Identifiability and separated signals. \(G_{XX}\) is positive definite, \(G_{YY}>0\), and \(\min_{j\in S_0}|c_{0j}|>0\). The nonzero coefficients do not approach zero with \(n\) in this theorem.

A3. A root-\(n\) expansion for the estimated Gram matrix. For a symmetric, mean-zero random matrix \(K(O)\) with \(E\|K\|_F^2<\infty\),

\[ \sqrt n(\hat{G}-G)=\frac1{\sqrt n}\sum_{i=1}^n K(O_i)+o_p(1). \tag{15} \]

The remainder is entrywise, equivalently in any matrix norm when \(p\) is fixed. The empirical objective is convex, and \(\hat{G}_{XX}\) is positive definite with probability approaching one. If artificial labels are retained, the independent sampling unit in (15) includes those labels.

A4. A root-\(n\) pilot and controlled weights. \(\sqrt n(\tilde{c}-c_0)=O_p(1)\), \(\gamma>0\), and \(0<\epsilon_n=o(n^{-1/2})\). The pilot in (13) satisfies the required rate under A2–A3. Another pilot may be substituted if its rate is established.

A5. Penalty separation. With the normalization in (14),

\[ \sqrt n\lambda_n\longrightarrow0, \qquad n^{(\gamma+1)/2}\lambda_n\longrightarrow\infty. \tag{16} \]

These limits can be in probability for a data-dependent penalty. For example, for \(\gamma=1\), \(\lambda_n=k\hat{G}_{YY}n^{-3/4}\) with fixed \(k>0\) satisfies (16) if \(\hat{G}_{YY}\to_pG_{YY}>0\); \(\epsilon_n=n^{-1}\) is admissible. The constant \(k\) still matters greatly in finite samples.

A6. Nondegenerate active uncertainty. For

\[ z(O)=K_{XY}(O)-K_{XX}(O)c_0, \tag{17} \]

define \(\Omega=\operatorname{Var}\{z(O)\}\). Its active block \(\Omega_{S_0S_0}=\operatorname{Var}\{z_{S_0}(O)\}\) is positive definite when nonsingular normal limits and Wald intervals are desired. A singular normal limit is possible without this last definiteness condition.

Assumption A3 is substantive. Consistency of a covariance estimate for artificial contrasts, by itself, is not enough to prove a root-\(n\) oracle result for a coefficient estimated from that covariance.

Covariance of the unpenalized pilot. On the event that \(\hat{G}_{XX}\) is invertible, the pilot in (13) satisfies the exact identity

\[ \tilde{c}-c_0 =\hat{G}_{XX}^{-1}\{(\hat{g}_{XY}-g_{XY})-(\hat{G}_{XX}-G_{XX})c_0\}. \]

A2–A3 therefore give

\[ \sqrt n(\tilde{c}-c_0) =\frac1{\sqrt n}\sum_iG_{XX}^{-1}z(O_i)+o_p(1) \ \overset{d}{\longrightarrow}\ N(0,V),\qquad V=G_{XX}^{-1}\Omega G_{XX}^{-1}. \]

The usual first-order covariance approximation for \(\tilde{c}\) is \(V/n\). Convergence of actual finite-sample covariance matrices requires additional moment control beyond convergence in distribution. In the fixed-feature case (25), \(z=h_Xu\) and \(\Omega=E(h_Xh_X^\top u^2)\). For estimated influence functions, use the full \(K\) from a justified expansion, including the contributions described in Section 8. In particular, \(\hat{G}_{XX}\) in (14) is not generally an inverse-covariance weight: if \(\Omega\) is nonsingular, \(V^{-1}=G_{XX}\Omega^{-1}G_{XX}\). Only under additional conditions such as \(E(u^2\mid h_X)=\sigma_u^2\) in the fixed-feature case does \(V=\sigma_u^2G_{XX}^{-1}\).

6.3 Theorem and proof

Theorem 1 (pointwise adaptive oracle property). Under A1–A6, the minimizer of (14) is unique with probability approaching one and satisfies

\[ P\{\operatorname{supp}(\hat{c})=S_0,\quad \operatorname{sign}(\hat{c}_{S_0})=\operatorname{sign}(c_{0,S_0})\}\longrightarrow1, \tag{18} \]

and

\[ \sqrt n(\hat{c}_{S_0}-c_{0,S_0}) \ \overset{d}{\longrightarrow} N\left(0,G_{S_0S_0}^{-1}\Omega_{S_0S_0}G_{S_0S_0}^{-1}\right). \tag{19} \]

Here \(S_0\) is the fixed true support. The limiting distribution in (19) coincides with that of the unpenalized estimator restricted to \(S_0\).

Proof. Put \(N_0=S_0^c\), \(s_0=\operatorname{sign}(c_{0,S_0})\) and \(e_n=\hat{g}_{XY}-\hat{G}_{XX}c_0\). Since \(g_{XY}=G_{XX}c_0\), A3 implies

\[ \sqrt n e_n=\frac1{\sqrt n}\sum_i z(O_i)+o_p(1), \qquad e_n=O_p(n^{-1/2}). \tag{20} \]

For \(j\in S_0\), A4 gives \(w_{jn}\to_p|c_{0j}|^{-\gamma}\), whence \(\lambda_nw_{S_0}=o_p(n^{-1/2})\) by A5. For the inactive coordinates,

\[ \sqrt n\lambda_n\min_{j\in N_0}w_{jn} =\frac{n^{(\gamma+1)/2}\lambda_n} {\{\max_{j\in N_0}(\sqrt n|\tilde{c}_j|\vee\sqrt n\epsilon_n)\}^{\gamma}} \ \longrightarrow_p\infty, \tag{21} \]

because the denominator is \(O_p(1)\) and the numerator diverges.

On the event \(\hat{G}_{XX}\succ0\), whose probability tends to one, define \(c^*_{N_0}=0\) and

\[ c_{S_0}^*=\hat{G}_{S_0S_0}^{-1} \{\hat{g}_{S_0Y}-\lambda_n w_{S_0}\odot s_0\}. \tag{22} \]

Since \(c_{0,N_0}=0\) and \(\hat{G}_{S_0S_0}^{-1}=O_p(1)\),

\[ c_{S_0}^*-c_{0,S_0} =\hat{G}_{S_0S_0}^{-1} \{e_{n,S_0}-\lambda_n w_{S_0}\odot s_0\} =O_p(n^{-1/2}). \tag{23} \]

By A2, \(\operatorname{sign}(c^*_{S_0})=s_0\) with probability tending to one; hence (22) satisfies the active Karush–Kuhn–Tucker equations. The inactive gradient is

\[ \hat{G}_{N_0S_0}c_{S_0}^*-\hat{g}_{N_0Y} =-e_{n,N_0}+\hat{G}_{N_0S_0}\hat{G}_{S_0S_0}^{-1} \{e_{n,S_0}-\lambda_n w_{S_0}\odot s_0\} =O_p(n^{-1/2}). \tag{24} \]

Together, (21) and (24) imply that every inactive gradient is smaller in absolute value than its penalty threshold \(\lambda_nw_{jn}\) with probability tending to one. The candidate therefore satisfies all Karush–Kuhn–Tucker conditions. Strict convexity gives \(\hat{c}=c^*\) with probability tending to one, proving (18).

Multiplying (23) by \(\sqrt n\) and applying (20), A5 and Slutsky’s theorem yields (19). Moreover, \(\hat{G}_{S_0S_0}^{-1}\hat{g}_{S_0Y}-c^*_{S_0} =\lambda_n\hat{G}_{S_0S_0}^{-1}(w_{S_0}\odot s_0)=o_p(n^{-1/2})\), which establishes equivalence to the unpenalized estimator restricted to the true support. \(\square\)

If \(S_0=\varnothing\), the inactive argument with \(c^*=0\) gives \(P(\hat{c}=0)\to1\) under the remaining assumptions; (19) is vacuous.

6.4 A transparent sufficient case

For fixed features as in (9), with \(E\|h\|^4<\infty\), a centered empirical covariance satisfies

\[ K(O)=h(O)h(O)^\top-G, \qquad z(O)=h_Xu. \tag{25} \]

Estimating the marginal means does not add a first-order term: the centering correction is the outer product of sample mean errors, of order \(O_p(n^{-1})\). Thus A3 follows directly from a central limit theorem for score products.

The active oracle covariance becomes

\[ V_{S_0}=G_{S_0S_0}^{-1}E(h_{S_0}h_{S_0}^\top u^2)G_{S_0S_0}^{-1}. \tag{26} \]

This sandwich form allows nonlinear conditional means and heteroskedastic residuals. Replacing it automatically by a residual variance times \(G_{S_0S_0}^{-1}\) would require additional conditions. The empirical sandwich, evaluated on the selected support, is pointwise consistent when support recovery holds. It does not supply a uniform post-selection confidence guarantee near zero coefficients.

6.5 A useful special case: orthogonal scores

If \(\hat{G}_{XX}\) is diagonal and the numerical floor is inactive, the solution is

\[ \hat{c}_j=\operatorname{sign}(\tilde{c}_j) \left(|\tilde{c}_j|-\frac{\lambda_n}{\hat{G}_{jj}|\tilde{c}_j|^\gamma}\right)_+. \tag{27} \]

Therefore a variable is selected exactly when

\[ |\tilde{c}_j|^{\gamma+1}>\lambda_n/\hat{G}_{jj}. \tag{28} \]

With \(\gamma=1\) and \(\lambda_n\asymp n^{-3/4}\), the selection threshold for \(|\tilde{c}_j|\) is of order \(n^{-3/8}\). Relative to root-\(n\) noise, this threshold increases only as \(n^{1/8}\). This slow separation explains why asymptotic support consistency can coexist with many false selections at realistic sample sizes when the penalty constant is small.

7. Tuning, interpretation, and inference

7.1 Cross-validation targets projection risk

For a training fold, estimate its scores, center them using training quantities, estimate its pilot and adaptive weights, and fit a penalty path. Evaluate held-out squared score residuals,

\[ \mathrm{CV}(\lambda)=\frac1n\sum_k\sum_{i\in I_k} \left\{\hat{h}_Y^{(-k)}(O_i) -\hat{c}_{\lambda}^{(-k)\top}\hat{h}_X^{(-k)}(O_i)\right\}^2. \tag{29} \]

Here the superscript means that nuisance quantities and centering values were obtained without the validation observations. A common grid of fractions of a fold-specific maximum penalty is also possible and is used in the simulation. In that case the tuning parameter in (29) is the fraction, and each fold converts it to its own absolute penalty.

This criterion is aligned with the defining projection loss. It estimates out-of-sample score prediction error rather than a response likelihood. It need not be interpreted as AIC or BIC. Training-only pilot weights are important: calling ordinary CV after computing adaptive weights once from the full sample leaks validation information into training.

Risk-optimal tuning and exact support recovery are different goals. Neither ordinary CV nor a one-standard-error rule automatically verifies (16). The oracle theorem applies to a data-selected penalty only after its required rates have been established. The simulation therefore reports the fixed-rate rule and CV separately.

7.2 Scaling and \(\gamma\)

The default examined here is \(\gamma=1\), a fixed positive exponent satisfying the theorem. It is not uniquely optimal. Larger \(\gamma\) increases the distinction between small and large pilot coefficients and can make selection more sensitive to a noisy pilot.

Without a floor, the penalty \(\sum_j|c_j|/|\tilde{c}_j|\) for \(\gamma=1\) is invariant to separate predictor rescalings. Under a response rescaling by \(b\), the quadratic loss is multiplied by \(b^2\), so the absolute penalty level must also be multiplied by \(b^2\). A rule proportional to \(\hat{G}_{YY}\) has this property. A numerical floor should be transformed with the corresponding coefficient scale if exact finite-sample invariance is required; an unscaled floor is only asymptotically negligible under A4.

7.3 What an oracle claim does and does not imply

Theorem 1 is pointwise at a fixed distribution with fixed nonzero signals. For a local coefficient \(c_{0j,n}=d/\sqrt n\), the pilot remains \(O_p(n^{-1/2})\). Under the strong penalty condition in (16), that local coefficient can be excluded with probability approaching one. There is no contradiction: A2 does not hold along that sequence. This illustrates the nonuniformity analyzed more generally by Pötscher and Schneider [7].

Consequently, a selected-variable confidence interval is not automatically a valid confidence statement for a variable chosen after looking at the data. The simulations of coverage concern fixed, prespecified active coordinates. The oracle theorem also provides neither false-discovery-rate control nor a conditional-independence interpretation.

If \(G_{XX}\) is singular, for example because two candidate scores are identical, individual coefficients and their support need not be identifiable. Combining redundant scores into a meaningful block, removing exact redundancy, or changing the estimand is necessary. Adding a fixed ridge penalty makes the optimization well behaved but changes the target and does not retain Theorem 1 unchanged.

7.4 Uncertainty in the estimated variance reduction

For a fixed covariate set and the joint Gram expansion in A3, differentiating (7) gives the first-order influence function

\[ \psi_{R^2}(O)= \frac{2c_0^\top K_{XY}(O)-c_0^\top K_{XX}(O)c_0-R^2K_{YY}(O)} {G_{YY}}. \tag{R1} \]

Thus the usual asymptotic variance is \(E(\psi_{R^2}^2)/n\), when nondegenerate. To obtain (R1), differentiate \(g^\top G_{XX}^{-1}g/G_{YY}\) and use \(d(G_{XX}^{-1})=-G_{XX}^{-1}(dG_{XX})G_{XX}^{-1}\). The same formula applies to a prespecified subset using its corresponding blocks. Estimated-score contributions enter through the complete \(K\), just as they do for coefficient uncertainty. At zero association the first derivative can vanish, so an ordinary nondegenerate normal interval is not justified by this formula alone.

For independent validation of a fitted rule, freeze its scores and coefficients conditional on the development sample. Let \(r_c=h_Y-c^\top h_X\), \(v_r=\operatorname{Var}(r_c)\) and \(v_Y=\operatorname{Var}(h_Y)\). The validation variance-reduction ratio has conditional influence function

\[ -\frac{(r_c-Er_c)^2-v_r}{v_Y} +\frac{v_r}{v_Y^2}\{(h_Y-Eh_Y)^2-v_Y\}. \tag{R2} \]

With finite fourth moments, its empirical variance divided by the validation sample size supplies a standard error for that frozen rule’s performance. This conditional interval does not measure variability from repeating model development, nor uncertainty about transport to a new population. These distinctions prevent a selected in-sample \(R^2\) from being reported with an unjustified fixed-model interval.

8. Estimated scores: the main additional statistical issue

8.1 Why covariance estimation needs its own expansion

For distribution-dependent scores \(h(O;P)\), the empirical covariance of estimated scores need not have influence function \(hh^\top-G\). Under a valid differentiable plug-in expansion, the \((j,k)\) entry has the form

\[ K_{jk}(o) =h_j(o)h_k(o)-G_{jk} +E\{\dot{h}_j(O;o)h_k(O)+h_j(O)\dot{h}_k(O;o)\}, \tag{30} \]

where \(\dot{h}_j(O;o)\) denotes the first-order change in the score function induced by contaminating \(P\) at \(o\). Equation (30) is a functional derivative formula under conditions permitting differentiation of the expectation. For nonsmooth scores a separate argument may be needed.

Theorem 1 still applies if the full \(K\) satisfies A3, but its covariance must use (17), including the additional terms. Merely observing that the original marginal estimators are asymptotically normal does not prove (15).

Cross-fitting can help control empirical-process dependencies. It does not by itself make the covariance moment insensitive to nuisance estimation. The often-used \(o_p(n^{-1/4})\) nuisance rate is not a generic solution here: it is useful when the relevant moment has the necessary first-order orthogonality and its remainder has a suitable product structure. For an arbitrary covariance of estimated influence functions, that orthogonality must be derived.

There are three concrete routes to a valid implementation: use fixed features for which (25) holds; derive the full plug-in expansion; or construct and justify an orthogonal estimator of the required Gram entries. The latter is an extension requiring score-specific work, rather than an assumption hidden inside cross-validation.

Corollary 1 (a constructive smooth-nuisance case). Suppose \(h(O;\eta)\) depends on a fixed-dimensional nuisance parameter, \(E h(O;\eta_0)=0\), and

\[ \sqrt n(\hat{\eta}-\eta_0)=n^{-1/2}\sum_i a(O_i)+o_p(1), \qquad E a=0,\quad E\|a\|^2<\infty. \]

Assume \(h\) is twice continuously differentiable near \(\eta_0\), \(E\|h(O;\eta_0)\|^4<\infty\), and integrable uniform envelopes for the first two derivatives of \(h\) and \(hh^\top\) justify a uniform law of large numbers and the Taylor remainders. Then the centered empirical covariance of \(h(O_i;\hat{\eta})\) satisfies A3 with

\[ K(O)=h(O;\eta_0)h(O;\eta_0)^\top-G +\sum_r D_r a_r(O), \quad D_r=E\left[\frac{\partial}{\partial\eta_r} \{h(O;\eta)h(O;\eta)^\top\}\bigg|_{\eta=\eta_0}\right]. \]

Proof. Taylor-expand the empirical uncentered second moment at \(\eta_0\). Its linear term is the empirical derivative multiplied by \(\hat{\eta}-\eta_0\); the empirical derivative converges to \(D_r\) and the quadratic remainder is \(O_p(n^{-1})\). Substituting the nuisance expansion gives the stated \(K\). The empirical mean of the estimated scores is \(O_p(n^{-1/2})\), so subtracting its outer product also changes the covariance by only \(O_p(n^{-1})\). The finite-moment conditions give the central limit theorem. \(\square\)

This corollary is one verifiable route to the estimated-score oracle theorem. It does not require independence between the nuisance estimate and the observations used in the Gram matrix. The cross-covariance between the two terms of \(K\) is retained automatically.

8.2 A constructive Gram expansion for win-based components

For fixed bounded antisymmetric kernels, the additional expansion required by the oracle theorem can be made explicit. Suppose coordinate \(j\) has \(a_j(o)=E\{k_j(o,O')\}\) and \(h_j=a_j/q_j\), with \(q_j>0\) either fixed or a win probability. Define \(B_{jk}=E(a_ja_k)\), so \(G_{jk}=B_{jk}/(q_jq_k)\). If \(q_j=E\{w_j(O,O')\}\), its influence function is

\[ r_j(o)=E\{w_j(o,O')+w_j(O',o)\}-2q_j; \tag{W1} \]

set \(r_j=0\) for a fixed denominator. Differentiating the outer distribution and the two inner reference distributions in \(B_{jk}\) gives

\[ L_{jk}(o) =a_j(o)a_k(o) +E\{a_k(Z)k_j(Z,o)+a_j(Z)k_k(Z,o)\} -3B_{jk}. \tag{W2} \]

Therefore the empirical pooled-score covariance has A3 with

\[ K_{jk}(o)= \frac{L_{jk}(o)}{q_jq_k} -G_{jk}\left\{\frac{r_j(o)}{q_j}+\frac{r_k(o)}{q_k}\right\}. \tag{W3} \]

Indeed, the numerator is a bounded statistic of triples of observations; its first-order Hoeffding projection is (W2). The denominator is a bounded statistic of pairs, giving (W1), and the delta method gives (W3). Diagonal and leave-one-out conventions differ by \(O_p(n^{-1})\). Antisymmetry centers the population components; empirical centering changes no first-order term. Thus a joint multivariate central limit theorem supplies A3. If the active variance is nondegenerate and the remaining assumptions hold, Theorem 1 applies to this win-component construction.

Centered bounded fixed features also fit this calculation by taking \(k_j(o,o')=f_j(o)-f_j(o')\) and \(q_j=1\). For comparisons involving unbounded features or estimated survival hazards, the corresponding moment and nuisance-expansion conditions must be verified separately. Joint normality of the marginal estimates is the starting point for PATED; (W3) verifies the stronger covariance expansion needed for selection in this particular nonparametric subclass.

9. Extensions with explicit guarantees

9.1 A finite-sample support certificate

The fixed-\(p\) theorem is not a theorem for arbitrary \(p\gg n\). A useful extension separates the algebra of selection from probabilistic assumptions on a growing-dimensional Gram estimator and pilot.

Let \(S_0=\operatorname{supp}(c_0)\), \(N_0=S_0^c\), \(e=\hat{g}-\hat{G}c_0\), and use the matrix maximum-row-sum norm. Suppose \(\hat{G}\) is positive semidefinite, \(\hat{G}_{S_0S_0}\) is positive definite, and

\[ \|e\|_\infty\le\delta,\quad \|\hat{G}_{S_0S_0}^{-1}\|_\infty\le M,\quad \|\hat{G}_{N_0S_0}\hat{G}_{S_0S_0}^{-1}\|_\infty\le L. \]

Let \(a=\max_{j\in S_0}w_j\), \(b=\min_{j\in N_0}w_j\), and \(\beta_{\min}=\min_{j\in S_0}|c_{0j}|\).

Proposition 5 (deterministic certificate). If

\[ \beta_{\min}>M(\delta+\lambda a), \qquad (1+L)\delta+\lambda La<\lambda b, \tag{31} \]

the adaptive quadratic estimator has a unique solution with exactly the correct support and signs.

Proof. Use the candidate (22). Its active maximum error is at most \(M(\delta+\lambda a)\), preserving the active signs. The inactive gradient in (24) is at most \((1+L)\delta+\lambda La\), giving strict inactive KKT inequalities. For a convex quadratic, these strict inequalities force any other minimizer to have zero inactive coordinates; positive definiteness of the active block then gives uniqueness. \(\square\)

This certificate remains meaningful when the full Gram matrix is singular. To turn it into a high-dimensional probability theorem, one must prove rates for \(\delta\), pilot-based \(a,b\), conditioning constants \(M,L\), and the minimum signal. A ridge or initial lasso pilot is a possibility, but no root-\(n\) or support guarantee for such a pilot is assumed without proof. The proposition shows exactly which additional work is needed.

9.2 Variable blocks and a multicomponent response

For a nominal variable with \(K\) categories, a scalar coding would impose an arbitrary ordering. Instead, represent it by \(K-1\) linearly independent centered indicators. The same principle permits several prespecified thresholds or features for a single variable. A nominal response can likewise be represented by a vector rather than an arbitrary numerical label.

Let \(H_j\in\mathbb R^{d_j}\) represent predictor variable \(j\) and \(V\in\mathbb R^q\) represent the response. Remove exact linear redundancies. For invariance to nonsingular within-variable recodings, whiten each marginal block and the response using their positive-definite covariance matrices. In the resulting coordinates, stack the predictor scores as \(H\) and consider

\[ C_0=\arg\min_C\frac12 E\|V-C^\top H\|_2^2, \tag{32} \]

with adaptive variable-block penalty

\[ \lambda_n\sum_j (\|\widetilde{C}_j\|_F\vee\epsilon_n)^{-\gamma}\|C_j\|_F. \tag{33} \]

Here \(C_j\) contains all predictor-feature rows for variable \(j\) and all response columns. Variable \(j\) is selected when \(C_j\ne0\).

Proposition 6 (fixed-block oracle extension). With a fixed number of blocks and fixed block sizes, positive definiteness, nonzero active block norms bounded away from zero, a joint root-\(n\) expansion for all estimated moments including whitening, a root-\(n\) pilot, and (16), the estimator has consistent block support and the same first-order distribution on active blocks as the known-active-block unpenalized estimator.

Proof. Vectorizing \(C\) gives a quadratic problem with Hessian \(I_q\otimes G_{HH}\), up to a coordinate permutation. On a nonzero block, the Frobenius norm is differentiable and its gradient has norm one, so the active penalty is \(o_p(n^{-1/2})\). For a zero block, the pilot norm is \(O_p(n^{-1/2})\), while its quadratic-loss gradient is \(O_p(n^{-1/2})\). The inactive group KKT condition compares that gradient’s norm with \(\lambda_nw_j\), which dominates by (16). Restricted consistency, the KKT argument, and the central limit theorem proceed as in Theorem 1. \(\square\)

After whitening, invertible recodings within a variable’s feature span induce orthogonal changes of coordinates. Both the squared loss and the Frobenius block penalties are invariant to those changes, provided the pilot and floor conventions are transformed coherently. Thus selection depends on the variable’s represented subspace rather than on an arbitrary basis. Adding genuinely new features changes that subspace and appropriately changes the estimand.

With \(\operatorname{Var}(V)=I_q\), a normalized variable contribution is

\[ D_j^{\mathrm{block}} =\frac1q\operatorname{tr} \{C_{0j}^\top G_{j\mid-j}C_{0j}\}, \tag{34} \]

the reduction in total residual score variance when the block is added to all other blocks, divided by the total response score variance. It is zero exactly when the block coefficient is zero under positive definiteness.

9.3 Bounded representations and observational data

For fixed bounded scalar features, the full unpenalized projection has influence \(G_{XX}^{-1}h_Xu\). If \(\|h_X\|\le B_X\) and \(|h_Y|\le B_Y\), then

\[ \|\operatorname{IF}_c(O)\| \le\|G_{XX}^{-1}\| B_X(B_Y+\|c_0\|B_X). \tag{35} \]

This gives bounded local contamination sensitivity when the Gram matrix is well conditioned. It is a concrete way to build a more stable association measure. It changes the representation and hence the estimand; it is not a free robustness correction for an unchanged unbounded-score target.

For claims data, the independent sampling unit should ordinarily be a patient rather than a claim line. Repeated records require an appropriate subject-level score construction or a separate cluster asymptotic argument. The present simulations use independent subjects. Censoring and missingness likewise require a specified observed-data score and identification assumptions; artificial assignment does not remove those requirements.

10. Historical-data simulation study

10.1 Design, estimands, and validation

The primary experiments now represent a large historical development cohort: \(n=5000,10000,20000\), \(p=20\), and three, four, or five active covariates. The three signal configurations have population \(R^2=0.08,0.12,0.15\), respectively. Each of these nine scenarios has 200 independent replications. Two additional scenarios at \(n=10000\), four signals, and \(R^2=0.12\) examine win-ratio scores and a censored baseline covariate. Each development cohort has an independently generated validation cohort of 5,000. In total there are 2,200 development–validation pairs, with seed 20261024.

These are newly generated mathematical experiments, not results from clinical records or from executing multipleOutcomes. They use Python 3.12.14, NumPy 2.3.5, SciPy 1.17.0, pandas 2.2.3, and scikit-learn 1.8.0. Scripts and replicate-level results are supplied. The old diagnostic studies are retained in Appendix A and are not pooled with these results. No runtime benchmark was conducted.

For the main mixed-variable design, \(X_1,\ldots,X_6\) are jointly standard normal with correlation \(0.5^{|j-k|}\). Remaining coordinates alternate independent standardized binary, Poisson, and normal measurements. The active set is the first \(s\) indices of \((1,7,8,9,12)\), with coefficients given by the first \(s\) entries of \((0.7,-0.6,0.5,0.4,-0.35)\). Let \(b\) denote this vector and \(G_X=\operatorname{Cov}(X)\). Generate

\[ Y=X^\top b+ a\left\{0.5X_{10}X_{11}+(0.8+0.3|X_1|)\varepsilon\right\}, \qquad \varepsilon\sim t_5\sqrt{3/5},\quad \varepsilon\perp X. \tag{H1} \]

The residual term in braces is orthogonal to every coordinate of \(X\) and has variance \(v_e=0.25+0.64+0.48\sqrt{2/\pi}+0.09\). For a requested population value \(r\), choose

\[ a^2=\frac{b^\top G_X b}{v_e}\left(\frac1r-1\right). \tag{H2} \]

This makes the population score-projection \(R^2\) exactly \(r\), rather than tuning it after inspecting simulations. Mean-contrast scores are used for the mixed measurements. The nonlinearity and heterogeneous residual variance mean that the linear score projection is not a correctly specified conditional-mean model. The construction isolates the planning problem; it is not intended as a full model of a claims database.

Fivefold CV uses 40 logarithmically spaced fractions from 1 to \(10^{-4}\) of each training fold’s maximum penalty. Score construction, centering, predictor scaling, pilot fitting, and adaptive weights are training-only. Adaptive weights use \(\gamma=1\) and floor \(n_{\mathrm{train}}^{-1}\) on the standardized-predictor coefficient scale. The compared rules are:

• Adaptive CV: minimum held-out squared score residual.
• Adaptive 1-SE: the largest penalty with mean fold loss no greater than the minimum plus one standard error of that minimum, estimated as the fold-loss standard deviation divided by \(\sqrt5\).
• Ordinary lasso CV: the same folds, score scaling, grid rule, and minimum-loss criterion with equal weights.
• Adaptive rate: \(\lambda_n=0.25\hat{G}_{YY}n^{-3/4}\), as a separate oracle-compatible rate, not a claim about CV.
• Full: all 20 covariates.
• Oracle: the true active set, with coefficients estimated from the development cohort.

Selection is determined by the original penalized fit, using \(10^{-8}\) as the nonzero threshold. For the primary variance assessment, selected coefficients are then refitted without a penalty on that selected set; the refit is frozen before external validation. The oracle is a benchmark, not an implementable method. For the full estimator, exact support is necessarily zero and FDP is large by construction; those entries describe its lack of selection, not a defect in its definition.

Exact support is the fraction of replications with \(\hat{S}=S_0\). TPR is \(|\hat{S}\cap S_0|/|S_0|\). FDP is \(|\hat{S}\setminus S_0|/\max(1,|\hat{S}|)\), averaged across replications. Excess risk is \((c-c_0)^\top G_{XX}(c-c_0)/G_{YY}\) here, normalized to the primary-score variance. By (P8), it is exactly the loss in population variance reduction relative to the optimal adjustment. The archive also records unnormalized excess risk for the older studies; the definitions must not be mixed.

Population gain is (P8) evaluated using the known data-generating covariance and the refitted coefficients. Validation gain is the independent-sample residual-variance ratio in Section 2.5. Apparent \(R^2\) is the refitted development-sample projection value. Percentages in the tables denote variance reductions; differences between them are percentage points. Reported \(\pm\) values are Monte Carlo standard errors of the simulation mean, not confidence intervals for a single historical study. For 200 replications the maximum Monte Carlo SE of an exact-support proportion is 0.035.

Direct KKT checks of the stated optimization objective have maximum violation \(3.7e-11\) across 28600 fitted paths/problems.

10.2 Selection and retained variance reduction

The first table shows the primary sparse rule, adaptive 1-SE followed by refitting, across all nine mixed-variable scenarios.

Development n	Signals	True R²	Exact support	TPR	Mean FDP	Population gain	Validation gain ± MC SE (pp)
5,000	3	8.00%	0.870	0.948	0.000	7.64%	7.55% ± 0.09
5,000	4	12.00%	0.745	0.926	0.000	11.43%	11.35% ± 0.11
5,000	5	15.00%	0.815	0.954	0.000	14.55%	14.53% ± 0.12
10,000	3	8.00%	0.980	0.992	0.000	7.92%	7.98% ± 0.06
10,000	4	12.00%	0.960	0.989	0.000	11.89%	11.99% ± 0.08
10,000	5	15.00%	0.970	0.993	0.000	14.91%	14.89% ± 0.08
20,000	3	8.00%	1.000	1.000	0.000	7.98%	7.95% ± 0.06
20,000	4	12.00%	1.000	1.000	0.000	11.98%	11.97% ± 0.06
20,000	5	15.00%	1.000	1.000	0.000	14.98%	14.99% ± 0.07

At \(n=5000\), the conservative rule sometimes misses a weaker active coordinate; the resulting population reduction can be below the requested 8%–15% optimum. At \(n=10000\), it usually retains the active set. In these simulations, all 200 replications in each \(n=20000\) scenario recover the support. This last observation is finite Monte Carlo evidence, not a guarantee: with zero observed failures in 200 replications, the approximate one-sided 95% upper bound for the failure probability is 1.5%.

The next table compares all methods at \(n=10000\), four active covariates, and population \(R^2=12\%\).

Method	Exact support	Mean FDP	Mean size	Population gain	Validation gain ± MC SE (pp)	Apparent R²
Adaptive 1-SE + refit	0.960	0.000	3.96	11.89%	11.99% ± 0.08	11.89%
Adaptive CV + refit	0.675	0.105	4.75	11.93%	12.03% ± 0.07	11.98%
Lasso CV + refit	0.010	0.523	9.28	11.85%	11.95% ± 0.07	12.06%
Adaptive rate + refit	0.140	0.263	5.64	11.90%	12.00% ± 0.07	12.01%
Full	0.000	0.800	20.00	11.82%	11.92% ± 0.07	12.09%
Oracle	1.000	0.000	4.00	11.96%	12.07% ± 0.07	11.95%

Here the adaptive 1-SE rule recovers the exact set in 96.0% of replications, compared with 67.5% for adaptive minimum-loss CV and 1.0% for ordinary lasso CV. Yet their refitted population variance reductions are all close to 12%. This distinction is central to the two planning questions: selecting unnecessary covariates and losing useful variance reduction are different errors. A method can estimate almost the same variance gain while selecting many more variables.

The rate-based rule illustrates the same limitation found in the earlier diagnostics: an asymptotically admissible penalty rate does not by itself give reliable finite-sample selection. Its exact recovery is 14.0% in this scenario. None of these comparisons establishes universal dominance of an adaptive rule or proves that the 1-SE choice satisfies the oracle theorem’s tuning-rate assumptions.

Figure 1. Adaptive 1-SE with selected-set refitting. Dotted lines show population optimal R²; error bars show 95% Monte Carlo intervals for mean validation gain, not patient-level uncertainty.

The refit matters. In the same scenario, adaptive 1-SE’s unrefitted coefficients achieve a mean population reduction of 10.31%, whereas the selected-set refit achieves 11.89%. Ordinary minimum-loss shrinkage can sometimes lower coefficient risk relative to its own refit, so refitting is not universally preferable. Both versions are retained in the result files. Apparent and validation summaries can differ in either direction in a finite simulation; external validation, rather than the sign of one observed difference, is what separates development from evaluation.

10.3 Statistic-specific examples

Win-ratio scores. Generate independent \(Z_1,\ldots,Z_{20},E\sim N(0,1)\) and a continuous primary endpoint \(Z_Y=\rho\sum_{j=1}^4Z_j+\sqrt{1-4\rho^2}E\). The representations are the standardized single-endpoint log-win-ratio null scores \(h_j=\sqrt3\{2\Phi(Z_j)-1\}\) and \(h_Y=\sqrt3\{2\Phi(Z_Y)-1\}\). For a bivariate normal pair, their covariance is \((6/\pi)\arcsin(\rho/2)\); consequently, choosing \(\rho=2\sin\{(\pi/6)\sqrt{0.12/4}\}\) gives \(R^2=0.12\). This identity follows by writing \(E\{\Phi(Z_j)\Phi(Z_Y)\}\) as a bivariate normal orthant probability with correlation \(\rho/2\).

The implementation estimates each reference CDF from the development data, or the training fold inside CV. Validation records are transformed with that frozen empirical CDF. Thus actual validation includes CDF-estimation effects. Population gains in this row are limiting-score benchmarks, while validation gains evaluate the fitted empirical-score rule. The experiment validates the uncensored, one-component win-ratio special case of (P7); it does not validate every hierarchical censoring rule.

Censored baseline event history. Generate an independent baseline event time \(T\sim\operatorname{Exp}(1)\) and censoring time \(C\sim\operatorname{Exp}(0.5)\). Use the standardized score \(h_1=\{\Delta-\widetilde{T}\}/\sqrt{2/3}\). This has mean zero and variance one, because the known-hazard martingale residual has variance \(E(\Delta)=2/3\). Use independent remaining covariate scores, four active coordinates \((1,7,8,9)\), and an additive standardized \(t_5\) residual scaled to give \(R^2=0.12\). This is an observed-score experiment with known nuisance quantities. It shows how a censored baseline record enters the projection; it does not establish an interpretation in terms of its unobserved latent event time or evaluate estimated-hazard uncertainty.

Score construction	Method	Exact support	Mean FDP	Population gain	Validation gain ± MC SE (pp)
Win ratio	Adaptive 1-SE + refit	1.000	0.000	11.96%	11.94% ± 0.06
Win ratio	Adaptive CV + refit	0.670	0.106	11.93%	11.91% ± 0.06
Win ratio	Lasso CV + refit	0.005	0.542	11.85%	11.82% ± 0.06
Win ratio	Full	0.000	0.800	11.82%	11.79% ± 0.06
Win ratio	Oracle	1.000	0.000	11.96%	11.94% ± 0.06
Censored baseline	Adaptive 1-SE + refit	0.925	0.000	11.85%	11.82% ± 0.07
Censored baseline	Adaptive CV + refit	0.660	0.121	11.92%	11.91% ± 0.07
Censored baseline	Lasso CV + refit	0.020	0.535	11.85%	11.83% ± 0.07
Censored baseline	Full	0.000	0.800	11.82%	11.81% ± 0.07
Censored baseline	Oracle	1.000	0.000	11.96%	11.94% ± 0.07

These examples support the score-construction principle while delimiting its validation: no current simulation establishes uniform performance over all endpoint types, missing-data mechanisms, or historical-to-trial shifts. The general estimated-score theorem requires its own Gram expansion for each proposed implementation.

11. Relationship to existing methods and a literature audit

11.1 Direct mathematical overlap

The most important identity for assessing novelty is already visible in (14). When the estimated scores are available as rows of a matrix, the procedure is exactly adaptive lasso applied to those scores. The variance-minimizing coefficient in (5) is also the usual optimal linear control-variate coefficient. The name or artificial-group interpretation does not change these identities.

Source	Directly relevant overlap	Difference from the present target
Zou (2006) [1]	Adaptive coefficient weights and an oracle argument based on different penalty behavior at zero and nonzero coefficients.	Here the target is a score projection and the Gram estimator may have additional first-stage influence terms.
Wang and Leng (2007) [2]	General penalization of a quadratic approximation around an initial estimator.	The curvature in this manuscript is a score Gram matrix; it need not be inverse coefficient-estimator covariance.
Leluc, Portier, and Segers (2021) [3]	Lasso selection of control variates, least-squares refitting, and support-recovery analysis. This is a close mathematical predecessor.	Their main estimand is a Monte Carlo integral. Here support and deletion contributions are themselves associational targets in an observed population.
South, Oates, Mira, and Drovandi (2023) [4]	Penalized regression for selecting control variates, including high-dimensional examples.	Controls are constructed for integration using derivatives of a target density; that construction is not required here.
Han, Shen, and Zubizarreta (2023), equation (9) [5]	Regression of estimated influence functions with a data-adaptive \(\ell_1\) penalty to choose information sources.	Their nonnegative ensemble weights and discrepancy-squared penalties address external-source bias; these are not inverse-pilot association-selection weights.
Yamada et al. (2014) [6]	Feature-wise kernelized lasso for nonlinear association-based feature selection.	HSIC Lasso approximates an outcome kernel matrix using predictor kernel matrices. It is not this projection of subject-level marginal influence scores.
Dai and Hughes (2012) [8]	Joint estimating-equation analysis of treatment-effect estimates for mixed outcome types.	Their target is trial-level surrogacy across actual trials, not sparse association induced by artificial groups in one population.
multipleOutcomes, pinned source [9]	Joint covariance estimation and an unpenalized variance-minimizing adjustment in pated().	The source examined does not implement this full artificial-group association definition, adaptive selection, and nuisance-aware oracle analysis as one workflow.
Li, Li and Luedtke (2023) [16]	External-data identification and estimation of future-trial covariate-adjustment efficiency, including partially observed outcomes.	Directly precedes the planning question; here the chosen target is PATED’s sparse marginal-score adjustment.
Zhou et al. (2026), preprint [19]	Historical-data prediction of martingale residuals for later log-rank/marginal-hazard-ratio adjustment, with a squared-correlation variance-reduction interpretation.	Close survival-specific overlap; it does not establish the full mixed-statistic adaptive-selection construction studied here.

In particular, Han et al. [5] make it inappropriate to claim that regressing one estimated influence function on others with a data-adaptive lasso penalty is itself an unstudied idea. Leluc et al. [3] likewise preclude claiming novelty for selecting variance-reducing controls by lasso.

11.2 What the search did and did not establish

The targeted search was conducted through 24 October 2026. It used both exact phrases and combinations involving artificial grouping, random grouping, influence-function association, adaptive lasso, score regression, control-variate selection, mixed outcomes, and covariance-based estimation. Primary paper PDFs, journal abstracts, author-hosted papers, and the supplied repository were examined. The reproducibility bundle contains a search log with queries, sources, and relevance assessments.

The additional search identified direct precedents for historical-data efficiency assessment [16] and a close 2026 survival preprint [19]. These substantially overlap the revised motivation. No exact published match was located for the complete combination of a common artificial group, a statistic-specific population association target, variable selection through inverse-pilot adaptive weights, analytic elimination of assignment randomness, and a Gram-expansion oracle result accounting for estimated scores. This is a bounded search finding, not evidence that no such publication exists. Some adjacent papers were available only as abstracts. The search was not a systematic citation-index review.

The plausible contribution is therefore a carefully specified synthesis: the association interpretation of artificial-group fluctuations; the reduction to an integrated score covariance; the allocation-noise calculation; a clear treatment of first-stage score estimation; and a coding-invariant block construction. Their individual proof ingredients are largely classical. The oracle theorem is an adaptation of established adaptive-penalty reasoning, not a new general principle. The elementary randomization and projection identities should also not be represented as unprecedented mathematical discoveries.

Any novelty claim must distinguish the mixed-statistic PATED selection workflow from external-data efficiency assessment, penalized control variates, and survival residual prediction. The present manuscript does not establish publication priority.

12. Implementation bridge to the supplied R framework

The supplied repository was examined at commit 3bfc7d361d5bca858dce7f49511961f47319fd42 [9]. In that source, pated() forms an optimal adjustment using a negative covariance-block solve. If its auxiliary contrast is \(\hat{\beta}_X\) and target contrast is \(\hat{\beta}_Y\), its sign convention corresponds to \(-c_0\) in the adjustment \(\hat{\beta}_Y-c_0^\top\hat{\beta}_X\).

asymptoticJointCovariance() builds covariance blocks from subject-aligned scores and model-specific bread matrices. The coefficient and covariance objects are useful for defining the original artificial-group route. The examined return object does not directly expose the complete pooled influence-score matrix required for the integrated route. Extracting or exposing that matrix, or implementing fixed-feature scores directly, is an additional implementation step.

An implementation should follow this order:

1. Specify the target. Record the population, independent sampling unit, marginal statistic or feature block for every variable, response representation, and observation assumptions.
2. Obtain compatible subject scores. Preserve subject alignment and a common coefficient scale. For example, a standardized log-rank statistic is not automatically on the same root-\(n\) coefficient scale as a fitted regression coefficient. Define its normalized score representation explicitly before combining covariance blocks.
3. Integrate grouping when possible. With pooled scores, use their centered covariance. If using refitted artificial-group models instead, establish their expansion and the covariance estimator’s expansion separately. If a package returns \(\widehat{V}_\beta\) estimating \(\operatorname{Cov}(\hat{\beta})\) for the common \(n\) subjects, normalize it as \(\hat{G}=n\pi(1-\pi)\widehat{V}_\beta\) under (2), before using the penalty normalization in (14). A covariance of contrasts is not inverse covariance of the projected coefficients.
4. Fit the adaptive projection. Use (13)–(14), or the block version (32)–(33), with a documented floor and scaling convention.
5. Tune without leakage. Recompute all estimated score ingredients and adaptive weights inside each training fold; validate with (29). Retain a separate theoretical statement for any fixed-rate penalty.
6. Answer both planning questions. Report the selected original variables and their statistic definitions. Refit on the selected set if that is the planned rule; freeze the fit and evaluate variance reduction on independent historical validation data. Report apparent \(R^2\), validated achieved reduction, and uncertainty separately. Prespecify the selected set for the trial and state whether PATED coefficients will be re-estimated there.
7. Use the appropriate uncertainty expansion. Apply (25) for fixed features and the full \(K\) for estimated scores. Do not interpret ordinary selected-model intervals as generally valid selective inference.

For scalar scores, glmnet can solve the required objective [10]. After centering, replace predictor-score column \(j\) by \(h_j/w_j\), fit squared-error lasso with coefficient \(\theta_j=w_jc_j\), and transform back. Use an explicit penalty grid, no additional standardization, and no intercept for already-centered pseudo-data. The family="gaussian" software label invokes a squared-error optimizer here; it does not assert a normal distribution for the original mixed-type variables. Using column rescaling also avoids ambiguity from internally rescaled penalty.factor values.

If only a positive-definite \(p\times p\) Gram matrix and \(\hat{g}\) are retained, write \(\hat{G}_{XX}=R^\top R\) and use \(p\) pseudo-observations with design \(\sqrt p\,R\) and response \(\sqrt p\,R\tilde{c}\), followed by the same column rescaling. Their squared-error objective is exactly (14) up to a constant. For CV, retain fold-specific score summaries and a valid held-out loss; a single full-sample covariance alone does not supply leakage-free CV.

Only the analytic joint-covariance adapters in Section 2 are included. This article supplies a mathematical and computational prototype, not a tested drop-in R extension of the repository.

12.1 Delivered R implementation

The companion pated_historical_selection.R implements the scalar workflow using glmnet and, for the Cox/log-rank component, survival. No artificial label is generated. Its public functions are:

Function	Role
mean_score(), binary_logit_score(), logmean_score()	Training-fitted analytic components for marginal mean, log-odds, and log-mean contrasts.
cox_null_score()	Pooled treatment-only null component from a fitted baseline cumulative hazard, on a prespecified horizon; select the Cox coefficient or proportional log-rank scale.
rank_win_score()	Standardized one-component win-ratio null component, with a training empirical reference distribution.
pairwise_score(), hierarchy_kernel()	Analytic pairwise projection for a specified win/loss kernel. The supplied hierarchy covers continuous/binary components and the source-audited longer-is-better time-to-event rule.
pated_select()	Training-only score construction and adaptive weights inside CV; minimum-loss or 1-SE selection, optional selected-set refitting, and retained CV summaries.
pated_validate()	Independent validation variance reduction and its conditional standard error (R2).
artificial_joint_covariance()	The component Gram matrix, asymptotic covariance, and estimator covariance with explicit allocation scaling.
pated_from_joint()	PATED’s covariance-block adjustment and \(R^2\) from compatible marginal estimates and their joint covariance.

The functions expect one row per independent participant, finite recorded values, and an identifiable scalar score system. They stop on unsupported missingness or degeneracy. A GEE/MMRM implementation can use the same builder interface once its subject-level null component is derived; no ready-made pooled GEE/MMRM builder is claimed in this delivery. The pairwise hierarchy helper is not a replacement for validating every comparison rule in the final trial specification.

run_historical_simulations.R implements all 11 primary data-generating scenarios and six analysis rules. Its default is 200 replications per scenario. R and NumPy use different random-number streams, so rerunning the full design in R should reproduce the population targets and Monte Carlo pattern, not identical table entries.

The executed same-input comparison used three development cohorts of size 5,000, separate validation cohorts, and identical fold assignments. Across adaptive minimum-loss CV, adaptive 1-SE, and ordinary lasso CV, all nine selected sets agreed between R and Python. The largest absolute penalized-coefficient difference was \(3.37\times10^{-9}\); refitted coefficients differed by at most \(2.53\times10^{-14}\), and validation gains by at most \(7.78\times10^{-16}\). The R objective’s maximum KKT discrepancy was \(2.06\times10^{-8}\). Additional checks verified the tied-rank and pairwise-component identities, the covariance-block PATED formula, and the pooled Cox martingale component against survival::residuals() to \(4.45\times10^{-16}\).

These checks used R 4.3.3, glmnet 4.1-8, and survival 3.5-8. All 11 R simulation scenarios were additionally executed with 20 replications each (220 development–validation pairs), with outputs retained in the archive. The table-producing 2,200-replication suite remains the Python run described in Section 10; the R implementation and its numerical validation are supplied alongside it. The complete multipleOutcomes package has not been installed or integration-tested here. Its audited code supplies the PATED construction and adapter mapping, while the delivered R functions implement the historical-data extension.

13. Discussion

The planning problem has two outputs: a covariate set compatible with the intended PATED analysis and a quantitative assessment of its variance reduction. Their common definition is the covariance of jointly asymptotically normal marginal-statistic fluctuations. Independent artificial labels make the historical-data target explicit, and the common null-score factorization removes the need to generate those labels in computation.

The fixed-dimensional oracle property concerns the support of this score projection. It requires identifiable nonvanishing signals, a suitable pilot, a joint root-\(n\) Gram expansion, and separated penalty rates. It neither detects every form of dependence nor establishes conditional independence or uniformly valid post-selection inference. In practical terms, sufficiently informative historical data can make the sparse fit behave as though the relevant score coordinates were known in advance, but only under the stated conditions. Joint normality of the original marginal estimates alone is not enough to establish this selection theorem; the additional covariance-estimation expansion is essential.

The large-sample simulations illustrate why selection and variance reduction should be evaluated separately. At moderate \(p\) and large \(n\), retaining unnecessary variables may have little effect on variance reduction even when exact selection is poor. Conservative selection can discard weaker useful variables, and strong shrinkage can reduce the achieved benefit unless coefficients are refitted. Estimated scores and independent validation remain relevant even when the historical sample is large.

PATED’s linear-combination formula supplies the interpretation of \(R^2\); it does not make the historical population identical to the future trial. Endpoint-specific observation assumptions, transportability, and the actual randomization design remain part of the planned analysis. The reserved PATED-results section can later document trial-specific properties without attributing them prematurely to the historical-data oracle theorem.

References

1. Zou, H. (2006). The adaptive lasso and its oracle properties. Journal of the American Statistical Association, 101(476), 1418–1429. DOI. Primary paper PDF.

2. Wang, H., and Leng, C. (2007). Unified LASSO estimation by least squares approximation. Journal of the American Statistical Association, 102(479), 1039–1048. DOI and journal abstract. Used as a conceptual quadratic-penalization predecessor; the theorem in this article is proved directly.

3. Leluc, R., Portier, F., and Segers, J. (2021). Control variate selection for Monte Carlo integration. Statistics and Computing, 31, 50. DOI. Author manuscript.

4. South, L. F., Oates, C. J., Mira, A., and Drovandi, C. (2023). Regularized zero-variance control variates. Bayesian Analysis, 18(3), 865–888. Author manuscript. Earlier versions used the longer title “Regularised zero-variance control variates for high-dimensional variance reduction.”

5. Han, L., Shen, Z., and Zubizarreta, J. R. (2023). Multiply robust federated estimation of targeted average treatment effects. Advances in Neural Information Processing Systems, 36. Proceedings paper, especially equation (9).

6. Yamada, M., Jitkrittum, W., Sigal, L., Xing, E. P., and Sugiyama, M. (2014). High-dimensional feature selection by feature-wise kernelized lasso. Neural Computation, 26(1), 185–207. DOI. Author manuscript.

7. Pötscher, B. M., and Schneider, U. (2009). On the distribution of the adaptive LASSO estimator. Journal of Statistical Planning and Inference, 139(8), 2775–2790. Author manuscript.

8. Dai, J. Y., and Hughes, J. P. (2012). A unified procedure for meta-analytic evaluation of surrogate end points in randomized clinical trials. Biostatistics, 13(4), 609–624. DOI. Journal abstract.

9. Zhang, H. multipleOutcomes source repository, examined at commit 3bfc7d361d5bca858dce7f49511961f47319fd42. Pinned repository; joint covariance implementation; PATED implementation.

10. Hastie, T., Qian, J., and Tay, K. An introduction to glmnet. Official documentation, accessed 21 October 2026. Used only for the implementation bridge.

11. U.S. Food and Drug Administration (2023). Adjusting for Covariates in Randomized Clinical Trials for Drugs and Biological Products: Guidance for Industry. May 2023. Official guidance PDF.

12. Ye, T., Shao, J., Yi, Y., and Zhao, Q. (2023). Toward better practice of covariate adjustment in analyzing randomized clinical trials. Journal of the American Statistical Association, 118(544), 2370–2382. DOI. Author manuscript.

13. Bannick, M. S., Shao, J., Liu, J., Du, Y., Yi, Y., and Ye, T. (2025). A general form of covariate adjustment in clinical trials under covariate-adaptive randomization. Biometrika, 112(3), asaf029. DOI. Author manuscript, with earlier title.

14. Ye, T., Shao, J., and Yi, Y. (2024). Covariate-adjusted log-rank test: guaranteed efficiency gain and universal applicability. Biometrika, 111(2), 691–705; published online in 2023. DOI. Author manuscript.

15. Reluga, K., Ye, T., and Zhao, Q. (2024). A unified analysis of regression adjustment in randomized experiments. Electronic Journal of Statistics, 18(1), 1436–1454. DOI. Author manuscript.

16. Li, X., Li, S., and Luedtke, A. (2023). Estimating the efficiency gain of covariate-adjusted analyses in future clinical trials using external data. Journal of the Royal Statistical Society, Series B, 85(2), 356–377. DOI. Author manuscript.

17. Nagelkerke, N. J. D. (1991). A note on a general definition of the coefficient of determination. Biometrika, 78(3), 691–692. DOI. Used to clarify that generalized regression coefficients of determination exist.

18. Schemper, M., and Henderson, R. (2000). Predictive accuracy and explained variation in Cox regression. Biometrics, 56(1), 249–255. DOI. Used to distinguish survival explained variation from PATED variance reduction.

19. Zhou, J., Liu, Q., Mo, M., and Xia, A. (2026). Improving power in randomized controlled trials with time-to-event endpoints: a risk-free approach. arXiv preprint, 2605.27711, version 1, 26 May 2026. Manuscript. Cited for its historical martingale-residual strategy and variance-reduction interpretation; not represented as a peer-reviewed publication.

Appendix A. Earlier diagnostic simulations

A.1 Reproducibility and analysis conventions

These earlier experiments are retained as diagnostics of the theory, rather than the primary historical-data design. Their signal strengths and dimensions differ from Section 10. All reported numbers were generated for this manuscript, not taken from published studies. The accompanying scripts and replicate-level CSV files reproduce the experiments. Across these six retained diagnostic studies, 6,110 independent data sets were generated; paired methods within a replication share the same data set. This earlier suite uses seed 20261021, with independent deterministic streams for each experiment and scenario. The simulations use Python 3.12.14, NumPy 2.3.5, SciPy 1.17.0, pandas 2.2.3, and scikit-learn 1.8.0. They evaluate the mathematical estimators; they do not execute or validate the R package multipleOutcomes. No runtime benchmark was conducted.

For scalar adaptive fits, the initial rule fixed before inspecting results was \(\gamma=1\), \(\epsilon_n=n^{-1}\), and \(\lambda_n=0.25\hat{G}_{YY}n^{-3/4}\). A coefficient with absolute value greater than \(10^{-8}\) is counted as selected. A weighted lasso was solved through an algebraically equivalent pseudo-data problem. Direct checks of the original objective’s KKT conditions gave maximum violation \(3.9e-11\) over 8280 fitted paths/problems in the primary suite. An exploratory sensitivity analysis is explicitly separated below.

Exact recovery means that the selected set equals the entire true support. TPR is the fraction of true signals selected. FDP is the fraction of selected variables that are null, defined as zero for an empty selection. Coefficient MSE means the squared Euclidean coefficient error, \(\|\hat{c}-c_0\|_2^2\), averaged over replications, rather than divided by the number of coefficients. Excess risk is \(E\{h_X^\top(\hat{c}-c_0)\}^2\), evaluated using the known population predictor covariance. This is an unconditional population projection criterion; it does not require a correctly specified conditional response model.

For a Monte Carlo proportion \(\hat{p}\) based on \(R\) replicates, the reported standard error is \(\{\hat{p}(1-\hat{p})/R\}^{1/2}\). Its maximum is 0.025 for \(R=400\), 0.029 for \(R=300\), 0.035 for \(R=200\), 0.040 for \(R=160\), and 0.014 for \(R=1200\). A displayed empirical rate of 1.000 does not mean its true probability is exactly one; for zero observed failures, the approximate one-sided 95% upper bound on the failure probability is \(3/R\).

A.2 Mixed predictors, nonlinear residual structure, and oracle behavior

The first six predictors are standard normal with correlation \(0.5^{|j-k|}\). The remaining predictors alternate between standardized binary, standardized Poisson with mean 2, and normal variables, independently of the first block and of one another. The active coordinates are \(1,7,8,9\), with coefficients \((0.7,-0.6,0.5,0.4)\). We generate

\[ Y=X^\top c_0+0.5X_{10}X_{11} +(0.8+0.3|X_1|)\varepsilon, \quad \varepsilon\sim t_5\sqrt{3/5},\quad \varepsilon\perp X. \tag{36} \]

The interaction term is orthogonal to every individual predictor, so the stated \(c_0\) is the population projection coefficient. The conditional mean is nonlinear and the residual variance is heterogeneous. The scores are centered raw variables, a special case of (9); the \(t_5\) errors have the moments needed here. The oracle comparator knows the four active coordinates but estimates their coefficients from each sample. The full estimator uses every predictor.

n	p	Method	Replicates	Exact support	Mean FDP	Excess risk
250	30	adaptive	400	0.010	0.495	0.05180
250	30	full	400	0.000	0.867	0.18599
250	30	oracle	400	1.000	0.000	0.02581
1000	30	adaptive	400	0.060	0.369	0.00960
1000	30	full	400	0.000	0.867	0.04308
1000	30	oracle	400	1.000	0.000	0.00582
4000	30	adaptive	400	0.240	0.217	0.00191
4000	30	full	400	0.000	0.867	0.01043
4000	30	oracle	400	1.000	0.000	0.00153
2000	100	adaptive	200	0.005	0.626	0.00725
2000	100	full	200	0.000	0.960	0.07241
2000	100	oracle	200	1.000	0.000	0.00294

At \(p=30,n=4000\), adaptive shrinkage substantially lowers excess risk relative to the full projection and approaches the oracle risk. However, with the initial small penalty constant, exact support recovery is only 24.0%. At \(p=100,n=2000\), it is 0.5%. This is a limitation of this finite-sample tuning rule, not evidence of reliable variable discovery merely because its rate satisfies the theorem.

The standard deviation of the first adaptive coefficient relative to the oracle standard deviation was 1.036, 1.017, 1.004 across \(n=250,1000,4000\) at \(p=30\). The corresponding pointwise 95% sandwich coverage was 90.0%, 95.0%, 94.8%, compared with 92.2%, 96.0%, 95.8% for the oracle estimator. These are prespecified active-coordinate coverages; no claim about selective coverage or simultaneous inference is made.

A.3 Exploratory penalty sensitivity and larger samples

After inspecting the preceding results, we added an explicitly exploratory comparison of constants \(k=0.25\) and \(k=1\) in \(\lambda_n=k\hat{G}_{YY}n^{-3/4}\). The original results were retained. This extension also increased \(n\) while keeping \(p=30\) and the data-generating process unchanged.

n	Penalty constant k	Replicates	Exact support	Mean null selections	Excess risk
1000	0.25	300	0.063	2.673	0.009329
1000	1.00	300	0.967	0.033	0.009736
4000	0.25	300	0.250	1.327	0.002010
4000	1.00	300	0.993	0.007	0.001991
16000	0.25	200	0.555	0.635	0.000417
16000	1.00	200	1.000	0.000	0.000427
64000	0.25	150	0.807	0.207	0.000095
64000	1.00	150	1.000	0.000	0.000098

Both rules have the same oracle-compatible rate. The larger constant gives much better support recovery here, at the cost of additional shrinkage bias. At \(n=1000\), the average first coefficient is 0.6970 for \(k=0.25\) and 0.6815 for \(k=1\), against truth 0.7. With \(k=0.25\), support recovery improves slowly as predicted by (28), reaching 80.7% at \(n=64000\). This experiment demonstrates sensitivity, not a universal recommendation that \(k=1\) is optimal.

Figure A1. Exact support recovery and excess projection risk in the exploratory sensitivity analysis. Error bars in the left panel are approximate 95% Monte Carlo intervals, clipped to \([0,1]\); all-success cells use a rule-of-three lower endpoint. The oracle has the true support by construction.

A.4 Cross-validated adaptive lasso versus ordinary lasso

Using (36) with \(n=1000,p=30\), we conducted 160 fresh replications. Both CV procedures used the same five subject folds within each replication. Each training fold recomputed its centering, Gram matrix, pilot, and adaptive weights. The grid contained 35 geometric fractions from 1 to \(10^{-4}\) of the training-specific maximum penalty \(\lambda_{\max}=\max_j|\hat{g}_j|/w_j\). The fraction minimizing total held-out squared residuals was refitted on the full sample. Ordinary lasso used \(w_j=1\). The fixed adaptive rule was included on the same data.

Method	Exact support	MC SE	TPR	Mean FDP	Excess risk
Adaptive, fivefold CV	0.738	0.035	1.000	0.087	0.00878
Ordinary lasso, fivefold CV	0.006	0.006	1.000	0.588	0.01881
Adaptive, fixed k = 0.25	0.037	0.015	1.000	0.382	0.00935

In this experiment, adaptive CV both reduced excess risk and selected fewer null variables than ordinary-lasso CV. It also improved support recovery over the initial fixed adaptive constant. The comparison establishes an empirical advantage in this specified setting, not universal dominance or a proof that CV satisfies (16). Its finite-sample success should not be conflated with a tuning-selection theorem.

A.5 Artificial randomization noise

We generated \(n=800\) observations with eight independent standard-normal predictors and \(Y=0.7X_1-0.5X_2+\varepsilon\), with independent standard-normal noise. Known, population-centered scores were used to isolate Proposition 4. Across 1200 replications, each method used the same original observations. The randomized methods used the \(W^2\)-weighted, uncentered score Gram matrix; they did not refit a nonlinear first-stage model.

Method	Predicted ratio	Empirical ratio	95% MC interval
Integrated	1.000	1.000	[1.000, 1.000]
One assignment, pi = 0.5	1.000	1.000	[1.000, 1.000]
One assignment, pi = 0.2	3.250	3.206	[2.942, 3.536]
Ten assignments, pi = 0.2	1.225	1.236	[1.176, 1.300]
One assignment, pi = 0.1	8.111	8.671	[7.838, 9.635]

Variance ratios refer to the first unpenalized projection coefficient, relative to integration. The intervals are paired bootstrap intervals for the Monte Carlo variance ratio, using 1000 resamples of complete replications. The balanced estimate is numerically identical to the integrated estimate. The observed ratios track (11)–(12); at \(\pi=0.1\), the larger finite-sample ratio is compatible with the asymptotic prediction within its Monte Carlo interval.

An additional enumeration of all 70 assignments with \(n=8,n_1=4\) verified (10) to maximum absolute error \(1.1e-16\). This is a numerical check of an exact identity, separate from the sampling simulation.

A.6 A censored outcome with mixed predictor types

We generated 20 independent standardized predictors alternating among binary, uniform, and three-level ordinal variables. The first three coefficients were \((0.10,0.08,0.06)\) and the remaining coefficients were zero. Let

\[ q(X)=0.45+0.10X_1+0.08X_2+0.06X_3, \qquad P(T\le1\mid X)=q(X). \tag{37} \]

The predictor ranges ensure \(0<q(X)<1\). Event times were generated exponentially with rate \(-\log\{1-q(X)\}\). Independent censoring times were exponential with hazard \(a_C\in\{0.35,1.2\}\). The analysis used the inverse-censoring-weighted observed score

\[ Z=I(T\le C,T\le1)\exp(a_CT). \tag{38} \]

The censoring survival function is treated as known in this experiment. Since \(E(Z\mid X)=P(T\le1\mid X)=q(X)\), projecting centered \(Z\) on the predictors has exactly the target coefficients in (37). The analysis does not fit the event-time generating distribution. A naive comparator treated subjects censored before an observed event as non-events; an unavailable-in-practice complete-outcome comparator used \(I(T\le1)\). All used the initial fixed adaptive penalty. There were 300 replications per scenario.

n	Censoring hazard	Method	Mean first coefficient (truth 0.10)	Coefficient MSE	Exact support
1000	0.35	IPCW	0.0953	0.00231	0.007
1000	0.35	Naive	0.0841	0.00226	0.000
1000	0.35	Complete outcome	0.0965	0.00167	0.003
1000	1.2	IPCW	0.0921	0.00468	0.013
1000	1.2	Naive	0.0631	0.00440	0.003
1000	1.2	Complete outcome	0.0969	0.00164	0.010
4000	0.35	IPCW	0.0979	0.00046	0.063
4000	0.35	Naive	0.0864	0.00074	0.063
4000	0.35	Complete outcome	0.0985	0.00032	0.077
4000	1.2	IPCW	0.0975	0.00092	0.077
4000	1.2	Naive	0.0659	0.00273	0.087
4000	1.2	Complete outcome	0.0991	0.00031	0.097

The censoring fraction before the event or horizon was approximately 22.5% at \(a_C=0.35\) and 55.4% at \(a_C=1.2\). Weighting largely preserves the intended coefficient as sample size grows, whereas the naive coefficient has appreciable attenuation, especially under heavier censoring. The IPCW method still has poor exact support recovery with the small fixed penalty constant, consistent with the earlier tuning limitation. More censoring increases variability; weighting is not a promise of lower finite-sample MSE in every setting.

This experiment establishes the behavior of the specified known-censoring-score construction. Estimating the censoring distribution would require its contribution to the Gram expansion. Moreover, although the raw coefficients in this example retain their full-data target, \(G_{YY}\) and normalized association strengths can depend on the censoring mechanism because the observed score has different variance. The example is not a proof that every survival score is observation-mechanism invariant.

A.7 Local alternatives and the limit of the oracle claim

For 20 independent standard-normal predictors, we set \(c_{01}=0.7\), \(c_{02}=2/\sqrt n\), all other coefficients to zero, and generated \(Y=X^\top c_0+\varepsilon\) with independent standard-normal noise. There were 400 replications per sample size, using the initial fixed adaptive penalty.

n	Local coefficient	Strong-signal selection	Local-signal selection	MC SE (local)
500	0.0894	1.000	0.755	0.022
2000	0.0447	1.000	0.682	0.023
8000	0.0224	1.000	0.545	0.025

The fixed strong signal was selected in every replication. Selection of the local signal declined as sample size increased, even though its standardized signal-to-noise scale stayed constant. This agrees with the threshold argument and shows why oracle properties cannot be read as uniform sensitivity to weak associations.

Reproducibility and rendering

The companion archive includes the R implementation, R simulation runner, R/Python verification script and shared fixtures, the new historical_simulations.py, all historical replicate-level and summary CSVs, earlier diagnostic scripts and results, document and figure builders, and an updated literature search log. The new study uses 2,200 development cohorts plus 2,200 independent validation cohorts. The original suite’s 8,510 data sets remain archived for reproducibility; the 6,110 retained diagnostic data sets are reported in Appendix A. The excluded median-score experiment is not used in the revised argument. No real patient records were used. The code validates the mathematical score estimators and optimization, not an installed R implementation of multipleOutcomes.

The Markdown and R Markdown files contain embedded figures. The YAML explicitly selects MathJax 3.2.2 with SVG output. In addition, the supplied HTML has every formula rendered at build time with MathJax 3.2.2 into embedded SVG, so reading it requires neither downloaded math fonts nor a runtime MathJax script. This is the reviewed version for RPubs preparation. The R Markdown file is provided as an editable source; it was not knitted in R in this environment.

Accents attach only to their base: for example, \hat{X}_{t}, \tilde{c}_{S_0}, and \bar{X}_{t}. Subscripts follow the accented base, rather than being included inside the accent. Composite objects use explicit braces. A default MathJax setting alone does not establish the cause of an observed alignment problem; font, CSS, and renderer interactions can also matter. The explicit SVG build and visual checks address the delivered document’s actual rendering. See the MathJax 3.2 configuration documentation. Nothing has been published to RPubs by this revision.