Subsampling and Selection Stability

Why Examine Selection Stability?

A predictor-selection procedure can return different models when the data change slightly. This can happen even when the same predictors, response, and selection rule are used: a modest change in the observed rows may alter which variables meet the rule for inclusion. Repeating the procedure on many subsamples provides a way to describe how sensitive its choices are to such perturbations.

The resulting summaries are about the observed data and the specified procedure. A predictor selected often is not thereby proven to be genuinely important, and a predictor selected inconsistently is not thereby proven irrelevant. Stability is useful descriptive information alongside the selection result, not a substitute for scientific reasoning or appropriate inference.

Definition: In repeated subsampling, each run draws \(m\) distinct observational units from the \(n\) available units, without replacement within that run, applies the same predictor-selection procedure, and records the selected model. Repeating this process \(B\) times gives \(B\) selected models. Different runs may contain some of the same units.

Unlike bootstrap resampling, a subsample does not include repeated copies of a row within a run. The rows are sampled without replacement inside each subsample. The repeated subsamples themselves are not required to be disjoint.

Define the Procedure Before Repeating It

A stability summary is meaningful only when the selection procedure is held fixed across runs. Specify the response, candidate predictors, any predictors that are forced in or excluded, the selection rule and its stopping conditions, and how the data are prepared. Apply the same choices to every subsample. A frequency for a forced-in predictor is guaranteed by design, so it does not measure whether the search would independently select that predictor. An excluded predictor is not assessed by the search.

The subsample size must also permit the planned procedure to run. For example, if a regression fit cannot be calculated for a particular subsample, the selection method needs a pre-specified way to handle that situation; silently omitting unsuccessful runs can distort the summaries. The observational units and the way they are sampled should respect the data structure. If rows are clustered or ordered in time, simply sampling individual rows may not preserve the relevant structure.

These are descriptive perturbation analyses, not a new inferential model. Their interpretation does not require treating the \(B\) selected models as independent samples from a population of models. In fact, subsamples commonly overlap, so their results are dependent.

Fix the analysis choices.
Specify the candidate predictors, data preparation, selection procedure, stopping rule, subsample size \(m\), and number of repetitions \(B\).

Draw and analyze each subsample.
For each run, sample \(m\) distinct observational units and apply the unchanged selection procedure.

Record each selected set.
Represent the selected predictors in run \(b\) by a set \(S_b\), for \(b=1,\ldots,B\).

Summarize consistency.
Calculate predictor inclusion rates and, when useful, pairwise model overlap. Interpret them in light of \(m\) and the procedure used.

Predictor Inclusion Rates

For a candidate predictor \(j\), define an indicator \(I_{bj}\) that equals 1 if predictor \(j\) is selected in run \(b\), and 0 otherwise. Its subsample inclusion rate is the proportion of runs in which it appears:

\hat{\pi}_j=\frac{1}{B}\sum_{b=1}^{B} I_{bj}

The rate ranges from 0 to 1. For example, \(\hat{\pi}_j=0.70\) means that the predictor appeared in 70% of the \(B\) models produced under this particular repeated-subsampling setup. It does not mean there is a 70% probability that the predictor is truly associated with the response. The rate can change if the data, candidate set, selection rule, or subsample size changes.

Worked Example: Calculating Inclusion Rates

Consider a hypothetical data set of 30 patients, with a continuous recovery score as the response and hours of rehabilitation, age, and treatment dosage as candidate predictors. A fixed selection procedure is run on \(B=5\) subsamples, each containing \(m=20\) patients. Assume all five fits can be carried out and the same selection rule is used each time. The selected predictor sets are:

\(S_1=\{\text{hours, age}\}\), \(S_2=\{\text{hours, dosage}\}\), \(S_3=\{\text{hours, age, dosage}\}\), \(S_4=\{\text{hours}\}\), and \(S_5=\{\text{hours, dosage}\}\).

Hours appears in all five models, age in two, and dosage in three. Thus:

\hat{\pi}_{\text{hours}}=\frac{5}{5}=1.00,\qquad \hat{\pi}_{\text{age}}=\frac{2}{5}=0.400,\qquad \hat{\pi}_{\text{dosage}}=\frac{3}{5}=0.600

Under this setup, hours was selected in every run, while age was selected in two of the five and dosage in three of the five. With only five runs, the rates can take only the values 0, 0.200, 0.400, 0.600, 0.800, or 1.00. These results describe variation across the five subsamples; they do not establish that hours is more scientifically important than the other predictors.

Measuring Overlap Between Models

Inclusion rates describe predictors one at a time. To summarize how similar two whole models are, use the Jaccard index. For selected sets \(S_r\) and \(S_s\), it is the number of shared predictors divided by the number of distinct predictors appearing in either model:

J(S_r,S_s)=\frac{|S_r\cap S_s|}{|S_r\cup S_s|}

The index is between 0 and 1. It is 1 when two nonempty models contain exactly the same predictors and 0 when they have no predictors in common. If both models are empty, define their overlap as 1 because they are identical; if only one is empty, their overlap is 0. These conventions make the measure interpretable even in those boundary cases.

To see the calculation, let \(a\) be the number of shared predictors, and let \(b\) and \(c\) be the numbers appearing only in the first and second models. The union contains \(a+b+c\) predictors and the intersection contains \(a\), so \(J=a/(a+b+c)\). This divides shared predictors by the total distinct predictors, rather than by the size of just one model. Therefore, exchanging the two models does not change the index.

Worked Example: Calculating Average Pairwise Overlap

Suppose a hypothetical analysis produces four selected models from four subsamples: \(S_1=\{A,B\}\), \(S_2=\{A,C\}\), \(S_3=\{A,B,C\}\), and \(S_4=\{B,C\}\). There are six distinct pairs of models. For the pair \(S_1,S_2\), the intersection is \(\{A\}\) and the union is \(\{A,B,C\}\), so \(J(S_1,S_2)=1/3\).

For \(S_1,S_3\), the intersection is \(\{A,B\}\) and the union has three predictors, giving \(J(S_1,S_3)=2/3\). For \(S_1,S_4\), the intersection is \(\{B\}\), giving \(1/3\). The remaining pairs give \(J(S_2,S_3)=2/3\), \(J(S_2,S_4)=1/3\), and \(J(S_3,S_4)=2/3\).

The average pairwise Jaccard index is therefore:

\overline{J} =\frac{1}{6}\left(\frac{1}{3}+\frac{2}{3}+\frac{1}{3} +\frac{2}{3}+\frac{1}{3}+\frac{2}{3}\right) =\frac{3}{6}=0.500

The average overlap in this set of runs is 0.500. It summarizes whole-model similarity, but it does not say which predictors are driving the overlap; the inclusion rates answer that different question. Reporting both types of summaries can therefore be more informative than reporting either alone.

How Subsample Size Affects Stability

The choice of \(m\) changes the perturbation being examined. A smaller subsample contains less of the observed data in each run and may lead to more variable selected models. A larger subsample retains more rows in each run and may yield selections closer to those from the full data. Neither pattern is guaranteed: the effect depends on the data, candidate predictors, selection procedure, and stopping rule.

Consequently, an inclusion rate has no full interpretation without its subsample size. A rate of 0.70 from half-sized subsamples is not automatically comparable to 0.70 from subsamples containing 90% of the data. The procedures are being challenged by different perturbations. When comparing analyses across choices of \(m\), report the sizes and treat the results as summaries under those separate settings, not as if they came from one identical experiment.

Worked Example: Interpreting Two Subsample Sizes

A hypothetical data set contains \(n=40\) observations. The same selection rule and the same candidate predictors are used in two separate exercises, each with \(B=10\) runs. In the first exercise, \(m=20\); predictor \(X\) appears in 8 selected models. In the second, \(m=32\); \(X\) appears in 5 selected models. Assume every run is successfully fitted and the only planned change is the subsample size.

The two inclusion rates are:

\hat{\pi}_{X,\,m=20}=\frac{8}{10}=0.800,\qquad \hat{\pi}_{X,\,m=32}=\frac{5}{10}=0.500

In these runs, \(X\) was selected more often with subsamples of 20 than with subsamples of 32. The correct conclusion is not that the smaller sample gives a more reliable model, or that the predictor has a 0.800 chance of being relevant. Rather, the procedure's selections for \(X\) differed under two perturbation sizes. Further investigation could examine the other selected predictors, the model overlaps, and whether a particular pattern in the data explains the difference.

With ten runs, each inclusion rate changes in increments of 0.100. Increasing \(B\) would make the empirical summaries less dependent on the particular set of random subsamples used, but it would not add new observational units or remove uncertainty about the underlying relationships. Because runs overlap, a simple binomial standard error based on treating all \(B\) outcomes as independent is not generally justified.

Practical Interpretation and Common Pitfalls

Stability summaries matter because they reveal whether a reported selected model is sensitive to modest changes in which observations are analyzed. They can prompt a careful account of competing predictors or a review of the modeling choices. They should be reported with enough detail for a reader to understand what was perturbed and what was held fixed.

Confusing frequency with truth. A high inclusion rate means a predictor was selected often by this procedure under this setup; it is not a probability that the predictor is truly relevant.
Omitting the subsample size. An inclusion rate cannot be interpreted properly without saying how many observations were used per run.
Changing the rules between runs. If the candidate set, preprocessing, or stopping rule changes, variation in selected models no longer reflects only the repeated subsamples.
Treating overlapping runs as independent. Subsamples can share observations, so the selected models can be dependent. Do not automatically use independent-run calculations to attach a standard error to a frequency.
Equating predictor frequency with model similarity. Predictor rates summarize one variable at a time; pairwise overlap summarizes whole selected sets. The two summaries answer different questions.

Key takeaway: Repeated subsampling applies one fixed selection procedure to many subsets of the observed data. Inclusion rates and model-overlap measures describe how its choices vary, conditional on the data, procedure, and subsample size.

Check Your Understanding

Answer each question using the definitions and interpretations in this tutorial.

A predictor appears in 18 of 30 selected models. Calculate its inclusion rate and state what the rate describes.
Two selected models are \(\{A,B\}\) and \(\{B,C\}\). Calculate their Jaccard index.
Why can repeated subsamples overlap even though sampling is without replacement within each subsample?
What is one reason an inclusion rate from subsamples of size 20 should not be treated as directly equivalent to the same rate from subsamples of size 90?
Why does a high inclusion rate not establish that a predictor is truly associated with the response?