Two-Step Sampling Procedure

• Two-step sampling procedure is a statistical framework that decouples variable selection from function estimation to achieve near-oracle bias-variance trade-off in high dimensions.
• The first step employs group Lasso for component selection, effectively identifying active predictors while managing sparsity challenges.
• The second step uses penalized least squares with Sobolev penalties to estimate smooth functions, optimizing estimation accuracy and computational efficiency.

A two-step sampling procedure is a statistical or algorithmic approach in which two distinct stages are executed sequentially, with each phase performing a specific complementary task. In high-dimensional models, survey sampling, and design-based or Bayesian inference, two-step procedures frequently partition variable screening, sampling, or selection from subsequent estimation, adjustment, or uncertainty propagation. This paradigm is especially prominent in high-dimensional nonparametric regression, sparse modeling, and post-selection inference, as well as in optimal survey allocation and semiparametric methodology.

1. Model Structure and Motivation in the Two-Step Paradigm

The central setting motivating the two-step procedure is the high-dimensional additive regression model: $$y_i = c^* + g^*(z_i) + u_i,\quad g^*(z_i) = \sum_{j=1}^d g^*_j(z_{ij}),$$ where $g^*_j$ lies in a class of smooth functions (e.g., a Sobolev space of order $v$). The primary challenge arises when the number of candidate predictors $d$ is large (even larger than the sample size $n$), but only a small subset $s^*$ of the $g^*_j$ are nonzero. One seeks both interpretability (by identifying active components) and minimax prediction efficiency. The two-step procedure capitalizes on the sparsity of active components by decoupling feature selection from function estimation (Kato, 2012).

2. Step One: Group Lasso-Based Component Selection

In the first step, variable selection is conducted via the group Lasso. Each function $g_j$ is represented using a basis expansion: $g_j(z) \approx \sum_{k=1}^m \beta_{jk} B_k(z),$ with $m$ basis functions per component. The high-dimensional model thus becomes

$$y_i \approx c + \sum_{j=1}^d \sum_{k=1}^m \beta_{jk} B_k(z_{ij}) + u_i,$$

or, in vector-matrix notation, $y_i \approx c + \sum_{j=1}^d x_{ij}'\beta_j$.

The group Lasso solves

$$\min_{c,\,\beta_1,\dots,\beta_d} \left\{ \frac{1}{n} \sum_{i=1}^n (y_i - c - \sum_{j=1}^d x_{ij}'\beta_j)^2 + \lambda_1 \sum_{j=1}^d \|\beta_j\|_2 \right\},$$

where the block-$\ell_2$ penalty $\|\beta_j\|_2$ induces sparsity at the group (component function) level.

This yields a selected set $T\subset\{1,\ldots,d\}$ of indices corresponding to non-zero estimated group coefficients. Note that selection may include false positives (redundant or inactive components) and false negatives (missed true components), particularly when signals are weak or predictors highly correlated.

The first step thus identifies the likely active additive functions but does not shrink or smooth within selected components, nor control the functional bias due to overpenalization.

3. Step Two: Penalized Least Squares with Sobolev Penalties

Given the selected set $T$, the second step restricts estimation to the reduced additive model and employs penalized least squares with Sobolev-type roughness penalties. The optimization problem is: $$\min_{c,\,\{g_j\}_{j\in T}} \Bigg\{ \frac{1}{n} \sum_{i=1}^n \Big( y_i - c - \sum_{j\in T} g_j(z_{ij}) \Big)^2 + \sum_{j\in T} \lambda_2^2 I(g_j)^2 \Bigg\},$$ where $$I(g_j)^2 = \int_0^1 \left[ g_j^{(v)}(z) \right]^2 dz$$ is the Sobolev norm or another roughness functional associated with the smoothness order $v$. The tuning parameter $\lambda_2$ can be selected by cross-validation or other methods.

The output is the function

$$\hat{g}(z) = \hat{c} + \sum_{j\in T} \hat{g}_j(z_j),$$

where each $\hat{g}_j$ is supported only on the selected subset and optimized for bias-variance tradeoff via the Sobolev penalty.

Critically, the complexity of this step (and the choice of $\lambda_2$) depends on $|T|$ rather than $d$, freeing smoothness regularization from the "curse of dimensionality."

4. Oracle Properties, Error Decomposition, and Adaptivity

A central theoretical result is an oracle inequality that decomposes the mean squared error: $$\|\hat{g} - g^*\|_2^2 + \sum_{j\in T} I(\hat{g}_j)^2 \lesssim s^* n^{-2v/(2v+1)} + (\text{error due to extraneous components}) + (\text{error due to missed components}),$$ where $s^*$ is the effective support size of $g^*$. The leading term is the minimax-optimal rate for estimation in $s^*$ univariate Sobolev classes, i.e., $s^* n^{-2v/(2v+1)}$.

Additional terms (the size and magnitude of $T\setminus T^*$ and of $T^*\setminus T$) quantify the explicit cost of model selection error in Step 1. If the selected set is near-minimal ($T\approx T^*$), the procedure attains the oracle risk. Notably, the error bound does not require perfect variable selection, but is adaptive to the selection quality.

Unlike approaches using simultaneous double penalization for sparsity and smoothness, here the $\lambda_2$ tuning in Step 2 does not depend on the ambient dimension $d$—removing log $d$-scale inflation often required by joint penalization.

5. Bias-Variance Trade-off and Implementation Considerations

The separation of variable selection and function estimation in the two-step approach mitigates the bias induced by shrinkage, which is typically unavoidable when a single regularizer controls both sparsity and roughness. Step one enforces group sparsity (i.e., feature selection), while step two induces smoothness and allows each selected function to be estimated with (nearly) optimal bias-variance efficiency.

Both group Lasso selection and penalized least squares with Sobolev-type penalties are highly tractable. Standard solvers (e.g., block coordinate descent, alternating direction methods) exist for group Lasso and smoothing spline estimation. The decoupled scheme considerably simplifies computation in large $d$ settings.

The method is robust to the overselection and underselection of variables—errors can be directly interpreted in the final error bound. This framework accommodates highly nonorthogonal designs and function bases, and enables aggressive tuning of the selection penalty to avoid excessive shrinkage.

6. Applications and Extensions

Two-step estimation procedures are broadly applicable across fields such as genomics, signal processing, and economics, where high-dimensional additive structure and sparse effects are prevalent. The improved interpretability (by explicitly selecting functions) and statistical efficiency (by decoupling model selection from functional smoothness) make such procedures especially suitable in modern "large $d$, small $n$" regimes.

The explicit adaptivity of the theoretical risk bound to imperfect selection motivates further research in adaptive selection methods, handling correlated predictors, and more aggressive thresholding. Extensions could target more general penalizations, further decoupling, or refined second-step function estimation (e.g., local polynomial methods). The structure also lends itself naturally to subsequent post-selection inference.

7. Summary Table: Two-Step Additive Model Estimation

Step	Method	Objective
1: Variable Selection	Group Lasso (block ℓ₂)	Identify active functions, enforce sparsity
2: Function Estimation	Penalized LS (Sobolev)	Estimate selected functions, control smoothness

The two-step approach thus represents a principled and computationally efficient solution for high-dimensional, sparsely supported additive modeling, achieving near-oracle risk rates and reducing shrinkage bias by explicitly separating feature selection from function estimation (Kato, 2012).