Answer

Is stepwise regression a good way to choose predictors?

Inspired by a question on Cross Validated · October 31, 2026

The short answer

Usually not. Stepwise search picks terms using the same data you then test them on, so the final model's p-values are too small, its R-squared too high and its coefficients exaggerated. On pure noise, backward elimination still ended with a significant overall test in 80% of simulated datasets. Choose predictors from subject knowledge, or use a penalized method such as the lasso checked by cross-validation.

What stepwise selection does

Stepwise regression is a family of automatic search rules for deciding which of many candidate predictors to keep in a linear model. Forward selection starts with no predictors and adds, one at a time, whichever improves the fit most. Backward elimination starts with every candidate and removes the least useful one at each step. Bidirectional search does both, so a term added early can be dropped later. Each move is judged by a rule: a p-value threshold for the t or F test, or an information criterion such as AIC. The search stops when no single addition or removal improves that rule.

The appeal is obvious: it is fast, it needs no judgement, and it hands back one tidy model. The trouble is what happens when you then read that model's summary table as if you had chosen it in advance.

Why the final model looks better than it is

Every test in a regression summary assumes the set of predictors was fixed before the data were seen. Stepwise selection breaks that assumption, because the terms that survive are exactly the ones that happened to fit this sample well, partly by chance. The consequences are well documented, for example in Frank Harrell's Regression Modeling Strategies:

R-squared is biased upward, because the search keeps whatever chance correlations it finds.
p-values are too small and confidence intervals too narrow. The reported t and F statistics no longer follow the distributions used to compute them.
Coefficients are biased away from zero. A predictor tends to be kept in the samples where its estimate came out large, so the surviving estimates overstate the effects.
The chosen set is unstable. With correlated predictors, a few changed observations can swap one variable for another, and different directions (forward versus backward) can end in different models.
Noise can masquerade as signal. With many candidates and modest sample sizes, the search often keeps several predictors that have no real relationship with the outcome. This is sometimes called Freedman's paradox, after a 1983 paper that demonstrated it.

Using an information criterion instead of p-values does not fix this. Dropping or adding one term by AIC is equivalent to a likelihood-ratio test at roughly α = .157, so a stepwise AIC search is still a sequence of significance tests, just with a looser threshold. For when AIC is the right tool, which is comparing a small set of models you chose in advance, see AIC vs BIC: which should you use?

See it in R and Python: selection on pure noise

The example simulates 100 observations of an outcome and 20 candidate predictors that are all independent random noise, so the correct model contains no predictors at all. It runs backward elimination by AIC with base R's step(), records what the final model's summary would claim, and repeats the whole experiment on 200 fresh datasets. For comparison it also records the overall F test of the full 20-predictor model, which was specified in advance and should be significant only about 5% of the time. The Python tab codes the same backward search by hand with NumPy and SciPy, using the same criterion as R.

set.seed(20836)
n <- 100; p <- 20

# Pure noise: y is unrelated to all 20 candidate predictors
make_data <- function() {
  d <- as.data.frame(matrix(rnorm(n * p), n, p))
  d$y <- rnorm(n)
  d
}

# Backward stepwise selection by AIC, starting from all 20 predictors
select_once <- function(d) {
  full  <- lm(y ~ ., data = d)
  final <- step(full, direction = "backward", trace = 0)
  s  <- summary(final)
  k  <- length(coef(final)) - 1
  pt <- if (k > 0) coef(s)[-1, 4] else numeric(0)
  fs <- if (k > 0) s$fstatistic else c(0, 1, 1)
  c(kept = k, sig_05 = sum(pt < 0.05), r2 = s$r.squared,
    f_p = unname(pf(fs[1], fs[2], fs[3], lower.tail = FALSE)),
    full_f_p = unname(pf(summary(full)$fstatistic[1], p, n - p - 1,
                         lower.tail = FALSE)))
}

# One dataset
one <- select_once(make_data())
round(one, 3)

# 200 datasets: how often does the selected model look "significant"?
res <- replicate(200, select_once(make_data()))
round(c(mean_kept = mean(res["kept", ]),
        mean_sig_terms = mean(res["sig_05", ]),
        share_any_sig_term = mean(res["sig_05", ] > 0),
        share_final_F_sig = mean(res["f_p", ] < 0.05),
        share_full_F_sig = mean(res["full_f_p", ] < 0.05),
        mean_r2_final = mean(res["r2", ])), 3)

Python

import numpy as np
from scipy import stats

rng = np.random.default_rng(20836)
n, p = 100, 20

def fit(X, y):
    # OLS with intercept: RSS, coefficient p-values, R-squared, overall F p-value
    Z = np.column_stack([np.ones(len(y)), X])
    beta, *_ = np.linalg.lstsq(Z, y, rcond=None)
    rss = np.sum((y - Z @ beta) ** 2)
    k, df = X.shape[1], len(y) - Z.shape[1]
    se = np.sqrt(np.diag(np.linalg.inv(Z.T @ Z)) * rss / df)
    pt = 2 * stats.t.sf(np.abs(beta / se), df)[1:]
    r2 = 1 - rss / np.sum((y - y.mean()) ** 2)
    f_p = stats.f.sf(r2 / k * df / (1 - r2), k, df) if k > 0 else 1.0
    return rss, pt, r2, f_p

def aic(X, y):
    # Same criterion as R's extractAIC for lm: n log(RSS/n) + 2 * coefficients
    return n * np.log(fit(X, y)[0] / n) + 2 * (X.shape[1] + 1)

def select_once():
    # Pure noise: y is unrelated to all 20 candidate predictors
    X, y = rng.normal(size=(n, p)), rng.normal(size=n)
    keep = list(range(p))
    while keep:  # backward elimination: drop the term whose removal lowers AIC most
        current = aic(X[:, keep], y)
        drops = [aic(X[:, [j for j in keep if j != i]], y) for i in keep]
        if min(drops) >= current:
            break
        keep.pop(int(np.argmin(drops)))
    _, pt, r2, f_p = fit(X[:, keep], y)
    full_f_p = fit(X, y)[3]
    return dict(kept=len(keep), sig_05=int(np.sum(pt < 0.05)), r2=r2,
                f_p=f_p, full_f_p=full_f_p)

# One dataset
print({k: round(float(v), 3) for k, v in select_once().items()})

# 200 datasets: how often does the selected model look "significant"?
res = [select_once() for _ in range(200)]
col = lambda k: np.array([r[k] for r in res])
summary = dict(mean_kept=col("kept").mean(), mean_sig_terms=col("sig_05").mean(),
               share_any_sig_term=np.mean(col("sig_05") > 0),
               share_final_F_sig=np.mean(col("f_p") < 0.05),
               share_full_F_sig=np.mean(col("full_f_p") < 0.05),
               mean_r2_final=col("r2").mean())
print({k: round(float(v), 3) for k, v in summary.items()})

    kept   sig_05       r2      f_p full_f_p 
   5.000    0.000    0.128    0.023    0.653 
         mean_kept     mean_sig_terms share_any_sig_term  share_final_F_sig 
             3.840              1.370              0.760              0.800 
  share_full_F_sig      mean_r2_final 
             0.030              0.125

The figures come from one seeded run of the R code. Python's random numbers differ from R's, so its figures differ slightly (for example, 82% of its datasets ended with at least one significant term, against 76% in R), but they show the same pattern.

What the simulation shows

In the first dataset, the honest full model gives no hint of a relationship (overall p = .653). Backward elimination keeps 5 of the 20 noise variables, and the final model reports R² = .128 with an overall F test of p = .023. Nothing in the printed summary warns that this "finding" was manufactured by the search.
Across 200 datasets, the search kept 3.8 noise predictors on average. In 76% of datasets at least one kept predictor had p < .05, and in 80% the final model's overall F test was significant, although every true coefficient was zero.
The prespecified full model behaved as it should. Its overall F test was significant in only 3% of datasets, close to the nominal 5%.
R-squared was inflated. The selected models averaged R² = .125 on data where the true value is zero.

Real data are less extreme than pure noise, but the same mechanism operates whenever weak or useless candidates sit alongside real ones: their chance associations get rewarded, and the summary table overstates the evidence. For a refresher on what those p-values are supposed to mean, see What do p-values and t-values mean?

Better ways to choose predictors

Decide from subject knowledge. Specify the model before looking at outcomes, based on theory, earlier studies and the design. Its tests are then valid, and a non-significant predictor can simply stay in.
Reduce dimensions without the outcome. Combine related predictors into a score, or summarise them with principal components, before fitting. Because the outcome plays no part, this does not bias the tests.
Use penalized regression for prediction. The lasso (which can set coefficients exactly to zero), ridge regression and the elastic net shrink estimates toward zero instead of making all-or-nothing choices, and the amount of shrinkage is chosen by cross-validation. They usually predict new data better than stepwise models.
Validate whatever you do. Judge predictive performance on held-out data or by cross-validation, and repeat the full selection procedure inside each fold or bootstrap sample. Bootstrapping the selection also shows how often each variable is chosen, which exposes how unstable the result is.
If you must report a stepwise model, present it as exploratory, give the full candidate list and the rules used, and do not treat its p-values or confidence intervals as confirmatory.

In Python, scikit-learn's SequentialFeatureSelector performs forward or backward selection using a cross-validated score rather than p-values, which helps with prediction but still leaves the final model's in-sample tests invalid. For more plain-language guides to building and reporting models, see the DASS blog.

How to report a stepwise regression in APA style (7th edition)

If a stepwise model appears in a paper, the method section should make the search reproducible and its limits clear: list every candidate predictor, the direction of the search, the selection rule and the software. In the results, report the final model's fit with its degrees of freedom and flag that the tests are not adjusted for selection. Using the first dataset above (the F value is a placeholder, since the example printed only its p-value):

"Twenty candidate predictors were entered into a backward stepwise linear regression using AIC as the selection criterion (the step function in R). The selected model retained five predictors, R² = .13, F(5, 94) = x.xx, p = .023. Because the predictors were chosen from the same data, these tests do not account for the selection process and the results should be regarded as exploratory."

In a table, give each retained predictor's coefficient (b), SE and p-value, and add a note naming the candidates that were dropped. Where possible, also report performance on held-out data, which is the figure readers can trust.

Related tools and guides

Working with your own data?

General answers only go so far. Send us your situation and we'll reply by email within two business days.

Ask your question

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.