Why Selection Can Change When the Data Barely Change
A regularized model can be useful for prediction even when its list of selected variables is sensitive to small changes in the data. This distinction matters in practice: a variable included in one fitted model may be excluded in a nearby fit, while the resulting predictions remain close. Conversely, a small-looking change in the data can sometimes have a noticeable predictive effect. Neither possibility should be assumed without checking.
Here, a perturbation means a specified small change to the observations or predictors, or a resampling of the observed cases. To assess stability meaningfully, state what changes and what is held fixed. For example, a comparison might hold the penalty parameter fixed and alter only the response, or refit the penalty parameter in each resample and thereby include tuning variability.
The examples use a Gaussian linear model with an intercept removed by centering. Predictors are centered and scaled so that each column has mean zero and squared length \(n\), where \(n\) is the number of observations. The intercept, if present, is not penalized. For a fixed penalty parameter \(\lambda>0\), the lasso minimizes
The selected set depends on the data, the predictor representation and scaling, the penalty, and the optimization specification. Stability comparisons are interpretable only when these choices are made explicit.
A Threshold Rule for Seeing Instability
A simple special case makes the mechanism precise. Suppose the columns of \(X\) are mutually orthogonal and each has squared length \(n\). Write \(z_j=x_j^\mathsf{T}y/n\). The lasso estimate for predictor \(j\) is the soft-thresholded value
Thus, predictor \(j\) is selected exactly when \(|z_j|>\lambda\). A value just below the threshold becomes selected after a sufficiently small change in the data. This does not mean the algorithm has malfunctioned: the objective itself changes its minimizer as the data change.
The formula follows by expanding the squared-error term. Under orthogonality, \(\lVert y-X\beta\rVert_2^2\), apart from a term not depending on \(\beta\), separates into a sum of terms \(n(\beta_j-z_j)^2\). For one coefficient, minimize \(\frac12(\beta_j-z_j)^2+\lambda|\beta_j|\). If \(\beta_j>0\), the derivative is \(\beta_j-z_j+\lambda\), giving \(\beta_j=z_j-\lambda\), valid when \(z_j>\lambda\). If \(\beta_j<0\), the derivative is \(\beta_j-z_j-\lambda\), giving \(\beta_j=z_j+\lambda\), valid when \(z_j<-\lambda\). At zero, the minimum condition is \(|z_j|\le\lambda\). Combining the three cases gives the threshold rule. This derivation relies on orthogonal, equally scaled columns; with correlated predictors, the coefficient equations are coupled and the rule does not apply one variable at a time.
Worked Example: A Response Perturbation Crosses the Lasso Threshold
Let \(n=4\), \(x_1=(1,1,-1,-1)^\mathsf{T}\), and \(x_2=(1,-1,1,-1)^\mathsf{T}\). Both columns are centered, have squared length \(4\), and are orthogonal. Set \(\lambda=0.50\), and initially let \(y=x_1+0.49x_2\). Determine the lasso selected set and compare it with the fit after replacing \(y\) by \(y'=x_1+0.51x_2\). Use the same predictors and penalty in both fits.
Method. Compute \(z_j=x_j^\mathsf{T}y/4\) for each response and apply the orthogonal-design threshold rule. The explicit vectors are \(y=(1.49,0.51,-0.51,-1.49)^\mathsf{T}\) and \(y'=(1.51,0.49,-0.49,-1.51)^\mathsf{T}\).
For the initial response, \(x_1^\mathsf{T}y=1.49+0.51+0.51+1.49=4.00\), so \(z_1=4.00/4=1.00\). Also, \(x_2^\mathsf{T}y=1.49-0.51-0.51+1.49=1.96\), so \(z_2=1.96/4=0.49\). Therefore \(\hat{\beta}_1=1.00-0.50=0.50\) and \(\hat{\beta}_2=0\), giving \(S=\{1\}\).
For the perturbed response, \(x_1^\mathsf{T}y'=1.51+0.49+0.49+1.51=4.00\), so \(z'_1=1.00\). Further, \(x_2^\mathsf{T}y'=1.51-0.49-0.49+1.51=2.04\), so \(z'_2=2.04/4=0.51\). The estimates are now \(\hat{\beta}'_1=0.50\) and \(\hat{\beta}'_2=0.51-0.50=0.01\), giving \(S'=\{1,2\}\).
The response change is \(y'-y=0.02x_2\). The fitted-value change is \(0.01x_2\), whose root mean squared value across the four observations is \(0.01\), because \(x_2\) has entries of magnitude one. Thus membership changes, but the fitted values change only slightly. Common mistake avoided: a coefficient exactly at or below the threshold is zero; the initial \(z_2=0.49\) is not selected merely because it is close to \(0.50\).
Correlated Predictors: Membership Can Move Between Variables
With correlated predictors, a small data perturbation can shift which variable receives credit for a shared signal. The following calculation uses standardized predictors with Gram matrix \(X^\mathsf{T}X/n\) and response correlations \(z=X^\mathsf{T}y/n\). The stated Gram matrix and correlations are sufficient for the lasso calculations; any component of \(y\) orthogonal to both columns does not affect them.
Worked Example: Two Correlated Predictors Exchange Lasso Selection
Suppose \(X^\mathsf{T}X/n=\begin{pmatrix}1&0.98\\0.98&1\end{pmatrix}\), and initially \(z=(0.80,0.78)^\mathsf{T}\). Fit the lasso with \(\lambda=0.10\). A small response perturbation changes the correlations to \(z'=(0.78,0.80)^\mathsf{T}\), leaving the predictors and \(\lambda\) fixed. Find both selected sets and the root mean squared difference between the two fitted-value vectors.
Method. For the initial fit, test a one-variable solution. If only variable 1 is active and its coefficient is positive, its coefficient is \(z_1-\lambda\), since its diagonal Gram entry is one. This gives \(\hat{\beta}_1=0.80-0.10=0.70\). The inactive-variable condition is that the absolute residual correlation be no greater than \(\lambda\). That correlation is \(z_2-0.98\hat{\beta}_1=0.78-0.98(0.70)=0.78-0.686=0.094\). Since \(0.094\le0.10\), the solution \((0.70,0)\) satisfies the lasso conditions, so \(S=\{1\}\).
For \(z'=(0.78,0.80)^\mathsf{T}\), the corresponding candidate is \((0,0.70)\). The inactive-variable residual correlation is \(0.78-0.98(0.70)=0.094\), again no greater than \(0.10\). Thus the perturbed fit has \(S'=\{2\}\).
The fitted-value difference is \(0.70(x_1-x_2)\). Since \(x_1^\mathsf{T}x_1/n=x_2^\mathsf{T}x_2/n=1\) and \(x_1^\mathsf{T}x_2/n=0.98\), its mean squared value is \(0.70^2(1+1-2(0.98))=0.49(0.04)=0.0196\). Taking the square root gives \(0.14\). The support has switched completely, while the predictions remain relatively close because the predictors are highly correlated. Common mistake avoided: high correlation does not make the selected variable uniquely determined; it can make different supports produce similar fitted values.
How the Elastic Net Changes the Comparison
For the same Gaussian setting, define the elastic net objective using \(0\le\alpha\le1\):
The lasso is the case \(\alpha=1\); when \(0<\alpha<1\), the added squared-coefficient penalty discourages extreme differences among correlated predictors, but does not guarantee identical selected sets. Under orthogonal columns, the coefficient is the soft-thresholded \(z_j\) divided by \(1+\lambda(1-\alpha)\). With correlated columns, coefficients must instead be found jointly.
Worked Example: Elastic Net Retains Both Correlated Predictors
Use the Gram matrix and initial correlations from the preceding example, but fit an elastic net with \(\lambda=0.10\) and \(\alpha=0.80\). Find the coefficients when both are positive, verify that this condition holds, and calculate the root mean squared difference in fitted values when the correlations are swapped.
Method. For positive coefficients, the first-order equations are \((X^\mathsf{T}X/n)\hat{\beta}+ \lambda(1-\alpha)\hat{\beta}+\lambda\alpha(1,1)^\mathsf{T}=z\). Here \(\lambda(1-\alpha)=0.02\), \(\lambda\alpha=0.08\), so the coefficient matrix has diagonal \(1.02\) and off-diagonal \(0.98\), with right-hand side \(z-(0.08,0.08)^\mathsf{T}=(0.72,0.70)^\mathsf{T}\).
Adding the two equations gives \(2(\hat{\beta}_1+\hat{\beta}_2)=1.42\), hence their sum is \(0.71\). Subtracting the second equation from the first gives \(0.04(\hat{\beta}_1-\hat{\beta}_2)=0.02\), hence their difference is \(0.50\). Therefore \(\hat{\beta}_1=(0.71+0.50)/2=0.605\) and \(\hat{\beta}_2=(0.71-0.50)/2=0.105\). Both are positive, so the assumed positive-coefficient equations are valid.
Swapping the two correlations swaps the coefficients to \((0.105,0.605)\). Their difference in fitted values is \(0.50(x_1-x_2)\), with mean squared value \(0.50^2(0.04)=0.0100\), and root mean squared value \(0.10\). In this calculation both predictors remain selected, although their coefficients change. Common mistake avoided: the elastic net does not automatically make correlated variables’ coefficients equal; the equations still reflect the difference in their response correlations.
Measuring Stability Across Perturbations
A stability assessment should report more than one selected set. One useful summary is the selection frequency for variable \(j\): the proportion of perturbation fits in which its coefficient is nonzero. Another is the Jaccard similarity for two selected sets \(S_a\) and \(S_b\), defined as the size of their intersection divided by the size of their union. A value of one means identical sets; a value of zero means disjoint sets. These summaries describe the chosen perturbation exercise, not a universal probability that a variable is truly relevant.
Worked Example: Summarizing Five Perturbation Fits
Suppose five fits, made under a stated perturbation procedure with fixed preprocessing and penalty, produce \(S_1=\{A,B\}\), \(S_2=\{A\}\), \(S_3=\{A,C\}\), \(S_4=\{A,B\}\), and \(S_5=\{A,C\}\). Calculate selection frequencies and the average pairwise Jaccard similarity.
Method. Variable \(A\) appears in all five sets, so its frequency is \(5/5=1.00\). Variable \(B\) appears in \(S_1,S_4\), giving \(2/5=0.40\). Variable \(C\) appears in \(S_3,S_5\), also giving \(2/5=0.40\).
There are \(\binom{5}{2}=10\) distinct pairs. In order, the Jaccard similarities for pairs \((1,2),(1,3),(1,4),(1,5),(2,3),(2,4),(2,5),(3,4),(3,5),(4,5)\) are \(1/2,1/3,1,1/3,1/2,1/2,1/2,1/3,1,1/3\). Their sum is \(2+4/3+2=16/3\), so their average is \((16/3)/10=8/15\), or \(0.5333\) to four decimal places. Variable \(A\) is consistently selected, but the alternatives \(B\) and \(C\) appear intermittently. Common mistake avoided: a high frequency for one variable does not mean the full selected set is stable; set-level similarity also reflects changing membership.
Perturbations can be designed to answer different questions. Resampling observations assesses sensitivity to which cases are represented. Small changes to measurements assess sensitivity to predictor or response values. Repeating the analysis with different random splits assesses sensitivity to a partition, if that is relevant to the intended workflow. State the perturbation mechanism, its size, and whether the penalty parameter is fixed or re-estimated. If preprocessing or variable screening is part of model fitting, say whether those steps are repeated within each perturbation.
For a prediction-focused assessment, compare predictions on the same cases across fits, using an appropriate scale such as root mean squared prediction difference. Do not compare predictions from different sets of evaluation cases as if the differences arose only from model perturbation. Selection frequency and prediction differences answer complementary questions: one concerns membership, the other concerns fitted outputs.
Common Traps and Practical Interpretation
- Treating one selected set as definitive. A single fit hides how close coefficients or residual correlations were to selection boundaries.
- Changing several choices without recording them. If penalty tuning, scaling, or factor coding changes as well as the observations, the resulting variability cannot be attributed to one cause.
- Equating unstable support with useless prediction. Correlated predictors can substitute for one another while producing similar fitted values, as the worked calculation showed.
- Equating stable support with a correct model. Repeated selection under a particular perturbation scheme is evidence of reproducibility under that scheme, not proof that the model contains every relevant predictor.
- Ignoring the size and direction of a perturbation. A stability claim is incomplete unless it states what was perturbed and by how much or by what resampling method.
Check Your Understanding
For each problem, distinguish what changes in the selected set from what can be concluded about the fitted values.
- Two orthogonal, standardized predictors have \(z=(0.62,0.39)\) and the lasso uses \(\lambda=0.40\). A response perturbation changes the correlations to \((0.62,0.41)\). Give both selected sets and the nonzero coefficient values. (Answer: initially \(S=\{1\}\), \(\hat{\beta}_1=0.22\); after perturbation \(S'=\{1,2\}\), \(\hat{\beta}'_1=0.22,\hat{\beta}'_2=0.01\).)
- For standardized predictors with correlation \(0.95\), suppose the lasso fit selects only predictor 1 with coefficient \(0.60\). A perturbation makes the symmetric fit select only predictor 2 with coefficient \(0.60\). Calculate the root mean squared difference between fitted values, using the same cases. (Answer: \(\sqrt{0.60^2(2-2(0.95))}=\sqrt{0.036}=0.1897\), to four decimals.)
- Across eight perturbation fits, a variable is selected in six. What is its selection frequency, and what does that number describe? (Answer: \(6/8=0.75\); it is the proportion selected under the specified eight-fit perturbation exercise, not a universal probability of relevance.)
- Two selected sets are \(S_a=\{A,B,D\}\) and \(S_b=\{A,C,D\}\). Calculate their Jaccard similarity and interpret it. (Answer: intersection size \(2\), union size \(4\), so \(2/4=0.50\); half the variables in their union are shared.)
- With the orthogonal-design elastic net, \(\lambda=0.20\), \(\alpha=0.75\), and \(z=0.30\) for a predictor, calculate its coefficient. (Answer: \((0.30-0.20(0.75))/(1+0.20(0.25))=0.15/1.05=0.1429\), to four decimals; it is selected because \(0.30>0.15\).)