Artificial Grouping and Sparse Influence-Score Association

Oracle theory, integration of randomization, and simulation evidence

Abstract

We study a statistic-specific association framework motivated by the joint covariance of estimators for heterogeneous outcomes. An artificial binary group is assigned independently of each subject’s observed variables. Although all population group contrasts are zero, their joint sampling fluctuations retain information about association. We define a multivariable association parameter by projecting the influence score for a response statistic onto the influence scores for candidate predictor statistics, and use an adaptive lasso penalty to select nonzero projection contributions. No parametric model for the conditional distribution of the response given the predictors is required. The framework applies to a single-arm sample or an observational population, provided that the sampling and observation assumptions needed for the chosen scores hold.

We establish the population interpretation, an exact finite-sample randomization identity for frozen scores, an asymptotic variance formula for unnecessary grouping noise, and a fixed-dimensional oracle theorem allowing estimated scores. We give a deterministic support-recovery condition, a variable-block extension for nominal and other multicomponent variables, and a special case showing why first-stage estimation can change standard errors. Seven simulation studies examine support recovery, prediction risk, randomization, nuisance estimation, censoring, and local alternatives. They support the analytic results while showing that an oracle-compatible tuning rate alone can yield poor finite-sample selection. The mathematical ingredients overlap substantially with adaptive lasso, least-squares approximation, and penalized control-variate methods. The intended contribution is their association interpretation and careful integration, rather than a claim that penalized influence-function regression is new.

Keywords: artificial grouping; mixed variable types; influence function; adaptive lasso; oracle property; control variates; observational data.

1. Scientific question and scope

The question is: which variables make an additional contribution to a specified measure of association with a response, after accounting for the other candidate variables? The measure need not coincide with Pearson correlation, an odds ratio, or a regression coefficient from a conventional response model. Its meaning is fixed by the statistics used to represent the variables.

For example, one variable may be represented by its event probability, another by a location statistic, and a censored outcome by a prespecified survival summary. Their subject-level influence scores place them in a common numerical space. Sparse projection in that space provides one coherent definition of selection across variable types. The resulting linear projection is a definition of an association parameter; it does not assert that the original response has a linear conditional mean.

Artificial grouping is a device for understanding this construction. Imagine repeatedly splitting the same population into two groups independently of all measured variables. Subjects with unusual values can push several group-contrast estimates in related directions. The covariance of those fluctuations reflects association in the underlying population. Once the relevant influence scores are available, the artificial grouping can be integrated out, so repeated random splitting is unnecessary for estimation.

The sampling population can be the patients in a single treatment arm or the population represented by a real-world data source. There is no requirement for a genuine treatment assignment. The target remains an associational quantity in the specified population. Artificial grouping does not itself identify a causal effect or correct nonrepresentative sampling.

Two distinctions are essential throughout:

1. The selected set is the support of a joint projection, not a list of all marginally associated variables. Redundant predictors may have zero additional contribution.
2. The oracle result concerns that projection and the chosen covariance estimator. It is not a claim of universal efficiency, uniformly valid post-selection inference, or detection of every possible form of dependence.

2. Population construction

2.1 Marginal statistics and common artificial groups

Let \(O_1,\ldots,O_n\) be independent observations from a population distribution \(P\). Initially, each candidate variable \(X_j\), \(j=1,\ldots,p\), and the response \(Y\) is represented by one regular scalar functional \(T_j(P)\) or \(T_Y(P)\). Let their centered influence functions be

\[ 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 difference between regular estimators of \(T_j\) in the two artificial groups. 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.

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} \]

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\)

Using independent artificial labels for different variables would generally eliminate the desired cross-covariances. A common label is part of the construction.

2.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.

2.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.

2.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 7.

As another example, for a continuous variable with positive density at its median \(m\), the median influence function is proportional to \(\operatorname{sign}(X-m)\). Its density factor is irrelevant to the scale-invariant quantities in (7)–(8), although estimating \(m\) can affect their sampling uncertainty. Section 6 makes that distinction explicit.

3. Integrating out artificial randomization

3.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.

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.

3.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.

4. Sparse estimation and the oracle theorem

4.1 Objective and adaptive weights

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\). Define a pilot and 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.

4.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 7.

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\) is fixed, with \(1\le |S_0|<p\).

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 \(\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} \]

the active covariance \(\Omega_{SS}=\operatorname{Var}\{z_S(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.

4.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)=\operatorname{sign}(c_{0S})\}\longrightarrow1, \tag{18} \]

and

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

This is the same first-order distribution as the unpenalized projection estimator using the same \(\hat{G}\) with the true support supplied in advance.

Proof. Write \(S=S_0\), \(N=S_0^c\), and \(e_n=\hat{g}_{XY}-\hat{G}_{XX}c_0\). By (15)–(17),

\[ \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\), \(w_{jn}\to_p|c_{0j}|^{-\gamma}\). For the fixed collection \(j\in N\), root-\(n\) consistency of the pilot and the floor condition imply

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

The denominator is \(O_p(1)\); a denominator near zero only strengthens the conclusion.

Construct a candidate \(c^*\) with \(c_N^*=0\) and

\[ c_S^*=\hat{G}_{SS}^{-1} \{\hat{g}_{SY}-\lambda_n w_S\mathbin{\odot}\operatorname{sign}(c_{0S})\}. \tag{22} \]

Using \(c_{0N}=0\), its error equals

\[ c_S^*-c_{0S} =\hat{G}_{SS}^{-1} \{e_{nS}-\lambda_n w_S\mathbin{\odot}\operatorname{sign}(c_{0S})\} =O_p(n^{-1/2}). \tag{23} \]

Positive definiteness, nonzero active signals, and (16) imply that the candidate has the correct active signs with probability approaching one. Thus (22) satisfies the active Karush–Kuhn–Tucker equations.

Its inactive gradient is

\[ \hat{G}_{NS}c_S^*-\hat{g}_{NY} =-e_{nN}+\hat{G}_{NS}\hat{G}_{SS}^{-1} \{e_{nS}-\lambda_n w_S\mathbin{\odot}\operatorname{sign}(c_{0S})\} =O_p(n^{-1/2}). \tag{24} \]

Equation (21) makes this gradient strictly smaller in absolute value, coordinate by coordinate, than the inactive thresholds \(\lambda_n w_N\), with probability approaching one. Therefore the inactive KKT inequalities also hold. The strictly convex objective then has the unique minimizer \(\hat{c}=c^*\), proving (18).

Finally, multiply (23) by \(\sqrt n\). The active penalty disappears by (16); (20), the multivariate central limit theorem, and Slutsky’s theorem yield (19). The true-support unpenalized estimate is \(\hat{G}_{SS}^{-1}\hat{g}_{SY}\), differing from (22) by \(o_p(n^{-1/2})\). This proves the oracle comparison. \(\square\)

In plain terms, pilot estimates of genuine signals stay away from zero, so their penalties become negligible. Pilot estimates of zero coefficients shrink toward zero, making their adaptive penalties large relative to estimation noise. The two requirements in (16) make these behaviors compatible.

The all-null case \(S_0=\varnothing\) follows from the same inactive KKT argument with candidate \(c^*=0\). The nonzero-signal requirement and active normal limit are then vacuous, and the probability of selecting no variables tends to one under the remaining conditions.

4.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=G_{SS}^{-1}E(h_Sh_S^\top u^2)G_{SS}^{-1}. \tag{26} \]

This sandwich form allows nonlinear conditional means and heteroskedastic residuals. Replacing it automatically by a residual variance times \(G_{SS}^{-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.

4.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.

5. Tuning, interpretation, and inference

5.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.

5.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.

5.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.

6. Estimated scores: the main additional statistical issue

6.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.

6.2 An analytic example with an estimated median

Let \(X\) have a unique median \(m\) and a continuous positive density in a neighborhood of \(m\), and assume \(E(Y^2)<\infty\). Let \(s(X)=\operatorname{sign}(X-m)\), and suppose \(r(x)=E(Y\mid X=x)\) is continuous at \(m\). Since \(E(s)=0\) and \(s^2=1\) almost surely,

\[ c(P)=E\{s(X)Y\} \]

is the projection coefficient of \(Y-EY\) on the normalized median influence score. Estimate it using the sample median and sample covariance. For even \(n\) and no ties, the sample signs have mean zero and variance one.

Proposition 5 (median-score projection influence). Under the preceding conditions, using the usual sample median,

\[ \operatorname{IF}_c(X,Y) =s(X)\{Y-r(m)\}-c(P). \tag{31} \]

Proof. For a fixed threshold, differentiating \(E\{\operatorname{sign}(X-m)Y\}\) with respect to the distribution gives \(s(X)Y-c\). Its derivative with respect to \(m\) is \(-2f_X(m)r(m)\). The median has influence \(s(X)/\{2f_X(m)\}\). Adding the threshold-estimation contribution yields (31). To justify replacing the threshold in the empirical term, the weighted indicator class \(\{YI(X\le t):t\text{ near }m\}\) has a square-integrable envelope and the usual VC-class stochastic equicontinuity; its \(L_2(P)\) increments vanish as \(t\to m\). Consequently its empirical-process difference at the consistent sample median is \(o_p(n^{-1/2})\). For odd sample sizes, the single central sign or centering correction is also negligible at first order under the stated moment condition. \(\square\)

If one incorrectly treats the estimated signs as fixed centered scores, the apparent influence is \(s(X)(Y-EY)-c\). The difference is \(\{EY-r(m)\}s(X)\), which need not vanish.

For a concrete case, let \(X\sim N(0,1)\) and

\[ Y=X^2+aX+\varepsilon,\qquad \varepsilon\perp X,\quad E\varepsilon=0,\quad \operatorname{Var}(\varepsilon)=\sigma^2. \tag{32} \]

Then \(m=0\), \(EY=1\), \(r(0)=0\), and \(c=a\sqrt{2/\pi}\). The correct and naive asymptotic variances are

\[ V_{\mathrm{correct}}=3+a^2+\sigma^2-c^2, \qquad V_{\mathrm{naive}}=2+a^2+\sigma^2-c^2. \tag{33} \]

At \(a=0.5\) and \(\sigma^2=1\), these equal \(4.09085\) and \(3.09085\). The simulation estimates \(r(m)\) nonparametrically for the corrected standard error. This example is a derivation for the present manuscript; the literature search did not establish whether this exact illustrative example has appeared elsewhere.

7. Extensions with explicit guarantees

7.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=\operatorname{supp}(c_0)\), \(N=S^c\), \(e=\hat{g}-\hat{G}c_0\), and use the matrix maximum-row-sum norm. Suppose \(\hat{G}\) is positive semidefinite, \(\hat{G}_{SS}\) is positive definite, and

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

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

Proposition 6 (deterministic certificate). If

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

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.

7.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 quantiles, time horizons, or prespecified 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{35} \]

with adaptive variable-block penalty

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

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

Proposition 7 (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{37} \]

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.

7.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{38} \]

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.

8. Simulation study

8.1 Reproducibility and analysis conventions

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 the seven studies, 8,510 independent data sets were generated; paired methods within a replication share the same data set. The primary 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\).

8.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{39} \]

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.

8.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 1. 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.

8.4 Cross-validated adaptive lasso versus ordinary lasso

Using (39) 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.

8.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.

8.6 Estimated median scores and standard-error correction

For (32), we used \(a=0.5\), \(\varepsilon\sim N(0,1)\), and 1200 replications at each sample size. The target was \(c=0.398942\). All methods used the same coefficient estimate based on the sample median. Only their standard errors differed. The corrected method estimated \(r(m)\) using Gaussian-kernel local averaging with bandwidth \(1.06\hat{\sigma}_Xn^{-1/5}\), and substituted it into (31). The known-boundary comparator used the true \(r(0)=0\).

n	Standard error	Empirical SD	Mean estimated SE	95% coverage	MC SE
500	Estimated signs treated as fixed	0.0921	0.0784	0.900	0.009
500	Corrected, estimated boundary mean	0.0921	0.0885	0.944	0.007
500	Corrected, known boundary mean	0.0921	0.0904	0.952	0.006
2000	Estimated signs treated as fixed	0.0447	0.0392	0.919	0.008
2000	Corrected, estimated boundary mean	0.0447	0.0446	0.953	0.006
2000	Corrected, known boundary mean	0.0447	0.0452	0.958	0.006

Ignoring median estimation understates uncertainty, as predicted by (33). The corrected intervals approach nominal coverage. This example shows that the projected coefficient can be consistent while the usual score-regression sandwich is inconsistent when estimated features are treated as fixed.

Figure 2. Left: empirical variance inflation from retaining artificial assignment noise, compared with the analytic factors. Right: coverage of 95% intervals for the median-score projection coefficient. The dashed line is nominal coverage.

8.7 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{40} \]

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{41} \]

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 (40). 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.

8.8 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.

9. Relationship to existing methods and a literature audit

9.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.

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.

9.2 What the search did and did not establish

The targeted search was conducted through 21 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.

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.

A claim of methodological novelty would need a focused comparison against the closest papers and a compelling application demonstrating why this particular association target answers a useful scientific question. The present results support a technically developed research direction, but do not alone establish publication priority or guarantee acceptance.

10. 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 (35)–(36), 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. Report the chosen estimand. Provide the selected original variables, coefficient orientations, \(R^2\), and deletion contributions, together with the exact score definitions. Report sensitivity to plausible tuning and representation choices.
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.

The current analytic engine does not provide every statistic through the same route: for example, the inspected quantile and Kaplan–Meier wrappers use a bootstrap route rather than the same analytic joint-score engine. Their use therefore requires a separate score or covariance-estimation justification. This article supplies a mathematical and computational prototype, not a tested drop-in R extension of the repository.

11. Discussion

Artificial grouping can provide a coherent interpretation of mixed-type association without requiring actual treatment arms. Its useful information is the joint covariance of contrast fluctuations, not the realized null contrasts themselves. Under a common influence-score representation, the construction reduces to a projection parameter that is independent of allocation proportion. For frozen scores, randomization can be integrated exactly; for known-score sampling estimators, its avoidable variance cost is explicit.

The adaptive estimator has a pointwise oracle property under identifiable fixed signals, a root-\(n\) pilot, a root-\(n\) Gram expansion, and separated penalty rates. The estimated-score condition is central. The median example and simulation show that omitting first-stage contributions can produce misleading uncertainty even when the point estimator behaves well.

The empirical advantages in this study are lower projection risk than the full fit, better support recovery than ordinary-lasso CV in the chosen mixed-predictor scenario, and removal of avoidable assignment noise. The observed defects are equally material: a formally admissible small penalty constant can select many null variables; weak local signals are not covered by pointwise oracle claims; heavy censoring reduces information; and coefficient meaning depends on the chosen score representation. These findings support using the framework as an explicitly defined association analysis, with its assumptions and tuning choices made visible.

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.

Reproducibility note

The companion archive contains simulations.py, sensitivity.py, replicate-level and summarized CSV results, numerical-check metadata, figure-generation and document-build code, this article’s source template, and literature_search_log.md. The sensitivity extension is identified as exploratory in both its script and the article. No real patient records were used. The Markdown and R Markdown versions contain embedded figures so they can be rendered without separate image files; a rendered HTML copy is provided for review and subsequent RPubs upload.