StepReg: A Comprehensive and Intuitive R Package for Stepwise Regression Analysis

Stepwise regression remains widely used for model selection, yet the field lacks a comprehensive, well-documented tool that supports diverse model families and selection strategies, implements multiple information criteria, and addresses overfitting and post-selection inference. We present StepReg, an R package that unifies stepwise selection across linear, generalized linear (e.g., logistic, Poisson, Gamma, negative binomial), and Cox models, supporting forward, backward, bidirectional, and best-subset search under various information criteria. StepReg also supports multivariate multiple linear stepwise regression, enabling simultaneous modeling of multiple dependent variables. Users can explore multiple strategies and information criteria within a single function call, and optionally combine them for more efficient and flexible model selection. To enhance robustness, StepReg provides an optional randomized forward selection mode to mitigate overfitting and a data-splitting workflow to improve the reliability of post-selection inference. The package further provides logging and visualization of the selection path, along with exporting results in common formats. A companion R package, StepRegShiny, has also been developed to provide a Shiny-based GUI for point-and-click analysis. Accuracy was assessed on public datasets by cross checking results against SAS. Together, these features provide a transparent, extensible framework that streamlines stepwise regression while promoting best practices.

1 Introduction

Stepwise regression remains a practical, widely used technique for model building and predictor screening across many applied domains (De Caires et al. 2022; Makepeace et al. 2022; Brandão-Dias et al. 2023; Caires et al. 2023; Seter et al. 2023; Tadros et al. 2023; Yang et al. 2023; Zhu et al. 2023; Ebrahimi et al. 2024; Ho et al. 2024; Naei et al. 2024; Yang et al. 2024a; Yang et al. 2024b), valued for its interpretability and ease of use when the candidate set is large. However, it carries well-recognized limitations, including a greedy, path-dependent search, susceptibility to overfitting, and challenges in conducting valid post-selection inference (Lewis-Beck 1978; Austin and Tu 2004; Whittingham et al. 2006; Flom and Cassell 2007; Smith 2018). These concerns, while important, can be mitigated through strategies such as principled stopping criteria (e.g., Akaike Information Criterion (AIC) (Darlington 1968; Judge et al. 1985; Hurvich and Tsai 1989; Al-Subaihi 2002) and Bayesian Information Criterion (BIC) (Sawa 1978; Judge et al. 1985)), cross-validation, data splitting with hold-out refitting to support more reliable inference, penalized regression as a complementary check, bootstrap or post-selection inference adjustments, and randomized variants of forward selection to reduce overfitting and stabilize variable selection (Stone 1974; Mitchell and Beauchamp 1988; Efron and Tibshirani 1993; George and McCulloch 1993; Tibshirani 1996; Hastie et al. 2009,; O’Hara and Sillanpää 2009). With these safeguards in place, stepwise procedures can preserve their practical utility and interpretability while mitigating key limitations.

Stepwise regression encompasses a variety of methodologies (termed “regression families”); each regression family is applicable to different types of regression models with its own complexities and considerations. For instance, linear stepwise regression (Draper and Smith 1998a; Burnham and Anderson 2002) is used for predicting continuous outcomes, selecting predictors based on information criteria like AIC, BIC, and p-values. Logistic stepwise regression (Menard 2002; Hosmer et al. 2013), on the other hand, is used for categorical outcomes, requiring adjustments to information criteria due to the non-linear nature of logistic models. The Cox proportional-hazards model (Cox 1972; Therneau and Grambsch 2000) is utilized for survival analysis, focusing on the time until an event occurs and requiring careful handling of both the proportional hazards assumption and censored data. Each of these methods necessitates careful considerations to balance model complexity, accuracy, and interpretability.

In terms of variable selection strategies, stepwise regression typically employs three widely used methods—forward selection, backward elimination, and bidirectional selection—along with the best-subset selection (see Supplementary Table 1). All strategies operate within a user-defined model scope and are influenced by predictor collinearity; they differ chiefly in path dependence and computational cost. Forward selection is a greedy, one-step-ahead procedure that can overfit if unchecked, yet it often attains comparable predictive performance against best-subset at substantially lower computational cost (Hastie et al. 2020). Backward elimination begins with a full model and is most reliable when the number of predictors (p) is modest relative to the number of observations (n). Bidirectional selection permits both additions and deletions, offering flexibility but remaining path dependent with the same sensitivities. Best subset selection exhaustively evaluates all subsets of the specified predictors under a chosen criterion. This approach is path-independent but computationally expensive (scaling roughly as \(2^p\)) and susceptible to overfitting without proper validation. The inherent complexity of model selection motivates a unified framework that streamlines applying and comparing tailored or combined strategies across datasets of varying size and type.

To evaluate model performance, various information criteria have been proposed (see Supplementary Table 2 and Supplementary Table 3 for details). These include information criteria-based metrics, such as AIC and BIC, as well as significance level-based (SL) methods like the F-test or the Wald test. Each metric offers distinct advantages and limitations, with the penalty coefficient size influencing model complexity: too small a penalty can lead to overfitting, while too large a penalty can cause underfitting. Researchers often need to combine different variable selection strategies and information criteria to identify the most optimal model.

Currently, ten open-source tools and four SAS procedures implement stepwise regression or best-subset selection (see Figure 1 and Supplementary Table 4). Although each has its strengths, capabilities are uneven and often narrow. Notably, to our knowledge, no existing tool implements randomized forward selection, a procedure that evaluates a random subset of remaining predictors at each step to reduce overfitting (Mentch and Zhou 2020; Chen et al. 2025), and most lack built-in data-splitting workflows to support more reliable post-selection inference. Coverage of regression families and information criteria is frequently limited. These constraints substantially restrict the model space that can be explored within any single tool, complicating comprehensive stepwise regression analyses. We summarize supported regression families, selection strategies, information criteria, and additional features by tool in Figure 1 and Supplementary Table 4.

Figure 1: Information criteria/metrics supported by 14 stepwise regression tools for each regression family

Of the 14 existing stepwise/subset selection tools we surveyed, eight are implemented in R (R Core Team 2025). The most widely used procedures—stats::step() (R Core Team 2025) and MASS::stepAIC() (Venables and Ripley 2002)—provide forward, backward, and bidirectional search but are oriented primarily toward AIC-type criteria. Other R packages cover narrower scopes. For example, leaps (Fortran code by Alan Miller 2024) implements best-subset selection for linear models with Mallows’ Cp (Cp) (Mallows 1973; Hocking 1976), R-squared, and adjusted R-squared (Darlington 1968; Judge et al. 1985). bestglm::bestglm() (McLeod et al. 2020) extends best-subset ideas to general linear model (GLM) families (e.g., logistic, Poisson, and Gamma), but remains focused on information criteria such as AIC, Schwarz Bayesian Information Criterion (SBC) (Schwarz 1978; Judge et al. 1985; Hurvich and Tsai 1989; Al-Subaihi 2002) and related variants, along with mean square error evaluated via cross-validation(Shao 1993, 1997; Davison and Hinkley 1997; Hastie et al. 2009). olsrr (Hebbali 2024) and blorr (Hebbali 2020) offer visualization utilities for linear and logistic regression workflows, respectively, while supporting a restricted set of information criteria; in our test environment, their Shiny graphical user interface (GUI) did not launch due to unresolved dependencies. My.stepwise (Hu 2017) prioritizes SL entry/removal and does not support alternative information criteria for GLM and Cox families.

Beyond R, two Python packages, stepwise-regression (Vijayakumar 2019) and openturns (Baudin et al. 2016), implement stepwise procedures but are limited to linear regression. stepwise-regression provides forward and backward selection using SL thresholds and has relatively limited documentation. openturns offers broader functionalities and up-to-date documentation; for stepwise routines, its manual indicates native support primarily for AIC and SBC, with other criteria potentially emulated via the penalty parameter, but users need to explore this further. Experienced users can use functions or modules from statsmodels (Seabold and Perktold 2010) or sklearn (Pedregosa et al. 2011) for stepwise regression. To our knowledge, no R or Python package currently implements randomized forward selection or offers built‑in data‑splitting workflows to mitigate overfitting and enable more reliable post‑selection inference.

In SAS, multiple procedures support stepwise selection across various model families, variable selection strategies, and information criteria. Logistic and Cox regression are fit via PROC LOGISTIC (SAS Institute Inc. 2019b) and PROC PHREG (SAS Institute Inc. 2019c), respectively. PROC HPGENSELECT (SAS Institute Inc. 2019a) supports generalized linear models, including Poisson, Gamma, logistic, and negative binomial (NB) regression. Selection strategies typically include forward, backward, and stepwise variants; by default, these use SL entry/removal thresholds, with information criterion and resampling options available in some procedures. While data-splitting or cross-validation can be configured, we are not aware of randomized forward selection being available natively in SAS. SAS is commercial software. Although graphical interfaces such as SAS Studio/Enterprise Guide exist, stepwise selection capabilities are primarily accessed through programmatic procedures, which may pose a barrier for some users.

Unlike univariate regression, which focuses on one outcome, multivariate regression looks at how several outcomes interact with each other and with the same set of predictors (Izenman 2013). This method can reveal complex interactions and correlations that are not apparent in univariate analyses, leading to more informed decision-making and better insights. This analysis is useful in fields like economics (Tintner 1964; Seyedaliakbar et al. 2019), where multiple economic metrics are analyzed, or medicine (Matsui et al. 2007; Peng et al. 2010; Bonnini and Borghesi 2022), where multiple health outcomes are studied. Multivariate regression uses techniques such as multivariate analysis of variance, canonical correlation analysis, and structural equation modeling to model and interpret these complex relationships. Notably, existing stepwise regression tools lack support for multivariate regression.

To address gaps in existing stepwise tools, we developed StepReg, a unified, versatile, and user-friendly R package for comprehensive stepwise regression analyses. StepReg supports six commonly used regression families —linear, logistic, Cox, Gamma, Poisson, and negative binomial regression—implements forward, backward, bidirectional, and best subset selection under various information criteria (e.g., AIC and BIC), and adds capabilities not found in other software tools, including multivariate regression (see Figure 1 and Supplementary Table 4). To enhance robustness, it provides optional randomized forward selection mode and a data-splitting workflow to mitigate overfitting and improve reliability of post-selection inferences. For purely predictive modeling, where the emphasis is on out-of-sample performance rather than formal inference, these inferential concerns are less central. We verify numerical correctness by cross-checking StepReg’s results against SAS on public datasets. To enhance accessibility, we provide an interactive Shiny-based GUI in the companion package StepRegShiny. Together, these features make StepReg a practical and comprehensive tool for stepwise regression analyses.

2 Materials and methods

2.1 Implementation

StepReg implements six commonly used stepwise regression families, including linear, logistic, Cox, Gamma, Poisson, and negative binomial regression in R. To streamline parameter acquisition for model comparison, StepReg utilizes the function lm() in stats for linear models, and function glm() in stats for logistic, Poisson, and Gamma regression models. For Cox proportional hazards models and negative binomial generalized linear models, StepReg employs function coxph() in survival (Therneau and Grambsch 2000; Therneau 2024) and function glm.nb() in MASS (Venables and Ripley 2002), respectively. Additionally, StepReg supports multivariate multiple linear stepwise regression by utilizing function lm() in stats to fit multivariate regression models and function anova() in stats to approximate F-statistics for model selection. Visualization of the selection process and output generation are implemented using ggplot2 (Wickham 2016) and flextable (Gohel and Skintzos 2024). The GUI is developed using the shiny (Chang et al. 2024) package and its extensions (Chang 2021; Attali 2022; Attali et al. 2024; Xie et al. 2024).

2.2 Selection strategies and information criterion formulas

StepReg implements four widely used variable selection strategies: forward selection, backward elimination, bidirectional selection, and the best subset selection (see Supplementary Table 1). Candidate models are evaluated using standard selection metrics. For information criteria such as AIC, BIC, and SBC, a lower value suggests a better model when adding or removing variables, whereas a higher adjusted R-squared reflects stronger model fit. Alternatively, for the SL, a predictor is included if its p-value falls below the entry threshold (sle) and removed if it exceeds the stay threshold (sls). Because definitions of some information criteria—particularly for linear models and small sample corrections—differ across sources and software (Darlington 1968; Judge et al. 1985; Hurvich and Tsai 1989; Hurvich et al. 1998), we documented the exact formulas implemented in StepReg in Supplementary Table 3 to ensure transparency and reproducibility.

2.3 Benchmark and datasets

We benchmarked StepReg against SAS using publicly available datasets: the mtcars dataset from the stats R package (R Core Team 2025), and four datasets bundled with the StepReg R package—remission (Lee 1974), lung (Therneau and Grambsch 2000; Therneau 2024), affairs (Kleiber and Zeileis 2008), and tobacco (Draper and Smith 1998b). To compare the outputs of multivariate regression, we utilized the SAS/IML script from a previous study (Al-Subaihi 2002) and modified the information criterion formulas. Specifically, for the SBC metric, we adjusted the penalty coefficients by multiplying them by the number of dependent variables to mitigate overfitting issues, as suggested by previous work (Anderson 2003). For other metrics, we made adjustments by either adding or multiplying a constant to align with the definition used in StepReg. The parameters used in the comparison are listed in Supplementary Table 5, together with the corresponding benchmark R script (2.3_Benchmark_StepReg.R) and SAS script (2.3_Benchmark_SAS.SAS) available in the Supplemental Materials.

3 Results

3.1 Workflow overview

StepReg provides four core functions, stepwise(), plot(), performance(), and report(), while the companion package StepRegShiny offers a Shiny-based GUI launcher, StepRegGUI(). Figure 2 summarizes the workflow. Analyses begin with stepwise(), which runs the selection procedure and returns a StepReg object recording the selection path and final model(s). This object can then be passed to plot() to visualize the path, to performance() to evaluate model performance on training and test data, and to report() to export results in common formats. Detailed arguments and default settings for each function are provided in Section 3.2 (Core functions and parameter definitions).

Figure 2: StepReg workflow overview

3.2 Core functions and parameter definitions

Below is a concise reference for the main arguments of the four core functions. See Table 1 for guidance on multiple dependent variables, survival outcomes, and predictor specification.

3.2.1 stepwise()

• Purpose

Runs the stepwise selection procedure and returns a "StepReg" object containing the selection path and final model(s).

• Key arguments

Model and data specification

formula: Model formula specifying a single response, multiple dependent variables, or survival outcomes on the left and predictors on the right (see Table 1).

data: data.frame containing the variables referenced in formula.

type (linear): Regression family. Options include "linear" (linear/Gaussian), "logit" (logistic/binomial), "poisson" (Poisson), "cox" (Cox proportional hazards), "gamma" (Gamma), and "negbin" (negative binomial).

Variable selection strategy

strategy: Variable selection strategy. Options include "forward" (forward selection), "backward" (backward elimination), "bidirection" (bidirectional selection), and "subset" (best-subset selection).

Selection criterion

metric: Model selection criterion. Supported values are "AIC", "AICc", "BIC", "SBC", "CP", "HQ", "adjRsq", "SL", "IC(3/2)", and "IC(1)". Note: Definitions can differ across sources; the exact formulas used in StepReg are enumerated in Supplementary Table 3.

Entry/stay controls (for SL-based metric)

sle (0.15): Significance level threshold used for forward and bidirectional selection. A predictor is added to the model if its p-value is below this cutoff.

sls (0.15): Significance level threshold used for backward and bidirectional selection. A predictor remains in the model if its p-value is below this cutoff.

Hypothesis test options

test_method_linear: For multivariate linear models: "Pillai"(default), "Wilks", "Hotelling-Lawley", "Roy". For univariate linear regression, the F statistic is used automatically.

test_method_glm: For GLMs: "Rao" (default) or "LRT". (Only "Rao" is available with the "subset" strategy.)

test_method_cox: Ties handling for Cox: "efron" (default), "breslow", "exact".

Constraints and robustness

include (NULL): Character vector of predictors that must be included in all models, namely, forced-in terms.

tolerance (1e-07): Threshold for multicollinearity checks; smaller values impose stricter screening.

weight (NULL): Optional numeric vector of observation weights ranging from 0 to 1.

Randomized forward selection, data splitting, and reproducibility

features_ratio (1): Proportion of remaining predictors sampled at random at each forward step (randomized forward selection). Set to 1 for standard forward selection.

test_ratio (0): Proportion of data reserved as a test set (e.g., 0.30). When > 0, performance() reports training/test performance metrics to support out-of-sample evaluation.

seed (123): Random seed used when test_ratio > 0 and/or features_ratio < 1.

Output control

best_n (3): Maximum number of top models displayed at each model size in best-subset search.

num_digits (6): Number of decimal places for printed results.

Table 1: Example formulas implemented in StepReg

Formulas	Descriptions
y ~ x1 + x2	Specify multiple predictors x1 and x2
y ~ .	Specify all predictors
y ~ . - x1	Specify all predictors except x1
y ~ x1 * x2	Specify main effects x1, x2, and their interaction
y ~ x1:x2	Specify the continuous-nested-within-class effects x1:x2
cbind(y1, y2) ~ .	Specify multiple response variables y1 and y2
y ~ . + 0 or y ~ . - 1	Fits model without an intercept (The Cox model does not use an explicit intercept, adding + 0, - 1, or + 1 to the formula will not change the results)
Surv(time, status) ~ . + strata(x1)	Cox regression with stratification variable x1

• Return value of stepwise()

stepwise() returns an S3 object of class "StepReg", a named list that records the analysis specification and results. The key components are:

argument

A data.frame summarizing all user-supplied settings (formula, type, strategy, metric, seeds, etc.).

variable

A data.frame describing all variables referenced in the model, including names, storage types, and roles (Dependent/Independent).

performance

A data.frame containing model-level performance metrics for each strategy-by-criterion. The model evaluation metrics vary by regression type:

For linear, Poisson, Gamma, and negative binomial regression models, the key performance indicators for both training and test data sets include adjusted R-squared (adjR-squared_train / adjR-squared_test applicable only to linear regression), mean squared error (mse_train / mse_test), and mean absolute error (mae_train / mae_test). A higher adjusted R-squared and lower MSE or MAE generally indicate a better model fit, while a large discrepancy between training and test metrics often signals overfitting, suggesting the model may not generalize well to unseen data. Note: If the adjusted R-squared for the test set exceeds 1, this is typically due to the small size of the test sample, which can produce unstable or uninterpretable values for this statistic.

For logistic regression models, the key performance indicators for both training and test data sets include accuracy (accuracy_train / accuracy_test), area under the receiver operating characteristic (ROC) curve (AUC) (auc_train / auc_test), and log loss (log_loss_train / log_loss_test). Higher accuracy and AUC values, combined with lower log loss, generally indicate stronger predictive performance. As with other models, substantial gaps between training and testing results suggest potential overfitting and warrant further investigation.

For Cox proportional hazards regression models, the primary performance indicators are concordance index (c_index_train / c_index_test) and time-dependent area under the curve (auc_sh). Higher values indicate stronger discrimination between risk groups. The c-index ranges from 0.5 (random prediction) to 1.0 (perfect prediction). Values below 0.5 typically suggest reversed risk direction or miscoding of events. In such cases, verify event coding and ensure that the scoring system aligns with the interpretation that higher scores correspond to higher risk.

overview

A nested list (by strategy and metric) giving the step-by-step selection path: which variables entered/left at each step, along with the corresponding criterion values (e.g., AIC, BIC, SBC).

detail

A nested list (by strategy and metric) with candidate step diagnostics: variables evaluated at each step, test statistics, p values, and selection decisions.

fitted models (within the strategy-metric hierarchy)

For each combination of selection strategy (e.g., “forward”, “backward”, “bidirection”, “subset”) and metric (e.g., “AIC”, “BIC”, “SBC”), the object stores the corresponding fitted model. These can be analyzed with standard S3 generics, summary(), anova(), plot(), or coef(), which dispatch to the underlying model class (e.g., coxph(), lm(), glm()). Users can also access components directly, e.g., result$forward$AIC$coefficients. Note: "coxph" objects in the survival package include detailed summaries (hazard ratios, standard errors, test statistics), whereas lm() / glm() display only basic coefficients by default. Use summary() to obtain standard errors, t/z statistics, p values, and R-squared.

• Printing behavior and downstream use

By default, a "StepReg" object prints a concise summary (e.g., selected model(s), criterion values, path length) while suppressing verbose components—argument, variable, performance, overview, and detail—from the console. These elements remain in the object for programmatic access and are utilized by the downstream functions such as plot(), performance(), and report().

3.2.2 plot()

• Purpose

Creates visualizations of the variable selection process from a "StepReg" object.

• Key arguments

x: A "StepReg" object.

strategy: Which strategy to display (e.g., "forward", "backward", "bidirection", "subset")

process: Selection view, one of "detail" (step-by-step path) or "overview" (summary of selected variables).

Outputs clear graphics for either the detailed path or a concise overview of included terms.

3.2.3 performance()

• Purpose

Produces a summary table of model performance across strategy-by-criterion combinations.

• Key arguments

x: A "StepReg" object.

If a hold-out split was specified, both training and test metrics are reported.

3.2.4 report()

• Purpose

Generates a formatted report from a "StepReg" object.

• Key arguments

x: A "StepReg" object.

name: Base file name for the report.

format: One or more output formats: "html", "docx", "pptx", and "rtf".

Reports include selection summaries, performance tables, and (optionally) figures.

3.2.5 StepRegGUI() in the companion R package StepRegShiny

It launches the Shiny GUI for point and click use. An online version is available at: https://junhuili1017.shinyapps.io/StepRegShiny/.

Together, these interfaces provide a concise print-view for readability and full programmatic access for analysis, visualization, benchmarking, and reporting. Comprehensive documentation, including vignettes and a user manual, is available with the package on CRAN: https://CRAN.R-project.org/package=StepReg.

3.3 Example usage

3.3.1 Data and preprocessing

This example used the lung dataset from the StepReg R package. We analyzed complete cases (n = 167) from patients with advanced lung cancer. The variables are:

• time: Overall survival time (days).

• status: Event status (1 = censored, 2 = death).

• age: Age (years).

• sex: Biological gender (1 = male, 2 = female).

• ph.ecog: ECOG performance score (0 = asymptomatic; 1 = symptomatic but completely ambulatory; 2 = in bed <50% of the day; 3 = in bed > 50% of the day but not bedbound; 4 = bedbound).

• ph.karno: Karnofsky performance score rated by physician (0 = worst, 100 = best).

• pat.karno: Karnofsky performance score rated by patient (0 = worst, 100 = best).

• meal.cal: Calories consumed at meals (kcal).

• wt.loss: Weight loss in last six months (pounds).

• inst: Institution code.

Preprocessing: We converted sex to a factor with levels “male” and “female” and excluded observations with missing data.

library(StepReg)
data(lung)
lung$sex <- factor(lung$sex, levels = c(1, 2), labels = c("male", "female"))
lung <- na.omit(lung)
str(lung)

'data.frame':   167 obs. of  10 variables:
 $ inst     : num  3 5 12 7 11 1 7 6 12 22 ...
 $ time     : num  455 210 1022 310 361 ...
 $ status   : num  2 2 1 2 2 2 2 2 2 2 ...
 $ age      : num  68 57 74 68 71 53 61 57 57 70 ...
 $ sex      : Factor w/ 2 levels "male","female": 1 1 1 2 2 1 1 1 1 1 ...
 $ ph.ecog  : num  0 1 1 2 2 1 2 1 1 1 ...
 $ ph.karno : num  90 90 50 70 60 70 70 80 80 90 ...
 $ pat.karno: num  90 60 80 60 80 80 70 80 70 100 ...
 $ meal.cal : num  1225 1150 513 384 538 ...
 $ wt.loss  : num  15 11 0 10 1 16 34 27 60 -5 ...
 - attr(*, "na.action")= 'omit' Named int [1:61] 1 3 5 12 13 14 16 20 23 25 ...
  ..- attr(*, "names")= chr [1:61] "1" "3" "5" "12" ...

3.3.2 Model specification and selection settings

We fit a Cox model using forward selection under two information criteria, AICc and SL with entry threshold sle = 0.05. Data are randomly split into 70% training / 30% test (test_ratio = 0.3). We constructed the survival response as Surv(time, status) and treated sex as a labeled factor; all remaining variables were considered as candidate predictors. See Table 1 for formula guidance and Section 3.2 (Core functions and parameter definitions) for argument details (e.g., strategy, metric, features_ratio for randomized forward selection, and test_ratio for train–test splitting).

3.3.3 Procedure

Following the StepReg workflow in Figure 2, we called stepwise() to run forward selection and obtain a "StepReg" object, employed plot() to visualize the selection path, used performance() to compute training/test metrics for the 0.7/0.3 split, and optionally generated a formatted report with report().

result <- stepwise(formula = Surv(time, status) ~ .,
                   data = lung,
                   type = "cox",
                   strategy = "forward",
                   metric = c("AICc", "SL"),
                   sle = 0.05,
                   test_ratio = 0.3)
result

$forward
$forward$AICc
Call:
coxph(formula = Surv(time, status) ~ ph.ecog + sex + wt.loss + 
    inst + ph.karno, data = data, weights = NULL, method = "efron")

               coef exp(coef)  se(coef)      z        p
ph.ecog    1.024028  2.784387  0.287734  3.559 0.000372
sexfemale -0.589911  0.554377  0.239437 -2.464 0.013749
wt.loss   -0.017024  0.983120  0.008828 -1.928 0.053799
inst      -0.030890  0.969582  0.015778 -1.958 0.050256
ph.karno   0.022345  1.022597  0.014420  1.550 0.121227

Likelihood ratio test=23.03  on 5 df, p=0.0003326
n= 117, number of events= 86 

$forward$SL
Call:
coxph(formula = Surv(time, status) ~ ph.ecog + sex, data = data, 
    weights = NULL, method = "efron")

             coef exp(coef) se(coef)      z       p
ph.ecog    0.4463    1.5625   0.1490  2.996 0.00274
sexfemale -0.5212    0.5938   0.2354 -2.214 0.02681

Likelihood ratio test=13.81  on 2 df, p=0.001001
n= 117, number of events= 86 

3.3.4 Outputs

The "StepReg" object consists of a fitted Cox model for each strategy-by-criterion combination. With sex coded as a factor, treatment contrasts yield a dummy variable (e.g., sexfemale) representing the effect of the female level relative to the male reference. For each selected model, stored results include coefficient estimates (coef), hazard ratios (exp(coef)), standard errors, z-statistics, p-values, and likelihood-ratio tests, enabling step-by-step review of the selection path and the final model summary. Standard S3 generics (e.g., summary(), coef(), anova()) work on the fitted objects. Individual components are also accessible via the $ operator.

coef(result$forward$SL)

   ph.ecog  sexfemale 
 0.4462914 -0.5212152 

result$forward$AICc$coefficients

    ph.ecog   sexfemale     wt.loss        inst    ph.karno 
 1.02402770 -0.58991071 -0.01702409 -0.03088999  0.02234546 

3.3.5 Performance

performance() produces a comparison table across strategy–criterion pairs. For Cox models, it reports training and test concordance (c-index) and time-dependent AUC, enabling side-by-side evaluation.

performance(x=result)

                                                           model
1 Surv(time, status) ~ ph.ecog + sex + wt.loss + inst + ph.karno
2                             Surv(time, status) ~ ph.ecog + sex
  strategy:metric c-index_train c-index_test auc_sh
1    forward:AICc        0.6429       0.4355 0.6676
2      forward:SL        0.6384       0.4113 0.6126

Practical notes: As a quality check, we verified event coding (death = 1, censored = 0) and confirmed that discrimination used the Cox linear predictor (higher = higher hazard); reversing either can spuriously drive the test-set C-index below 0.5. This example is for illustrative purposes only—we do not interpret test-set metrics further. In real analyses, if the test set C-index ≈ 0.5 (no discrimination) or < 0.5, we recommend rechecking event coding and risk-score direction; quantifying uncertainty with a bootstrap confidence interval (CI); assessing covariate and event-time distribution shift between training and test sets; running k-fold cross-validation or bootstrap to evaluate stability/overfitting; reviewing proportional-hazards and calibration diagnostics; and, under heavy censoring, also examining time-dependent AUC or the censoring-adjusted C-statistic (Uno et al. 2011). If issues persist, consider model simplification or penalization.

3.3.6 Visualization

To illustrate the selection process, we used the StepReg::plot() and arranged panels with cowplot::plot_grid() to produce two views (Figure 3): (i) a step-by-step path highlighting the top candidate at each step under both AICc and SL (selected entries shaded in green), and (ii) an overview summarizing the variables ultimately retained, with AICc and p-values (SL) tracked across steps. In the forward strategy, at each step the remaining predictors are evaluated, and the one that yields the lowest AICc or the smallest p-value for SL is considered for inclusion. A variable is accepted only if the AICc value decreases relative to the previous step or if the p-value for SL is below the entry threshold sle. If no candidate meets these conditions, the procedure terminates and the current model is retained.

library(cowplot)
plot_list <- setNames(
    lapply(c("forward"), function(i){
    setNames(
      lapply(c("detail", "overview"), function(j){
        plot(result, strategy=i, process=j)
      }),
        c("detail", "overview")
    )
    }),
    c("forward")
)
cowplot::plot_grid(plotlist = plot_list$forward, ncol = 1, rel_heights = c(2, 1))

Figure 3: Visualization of variable selection: detailed and overview perspectives

3.3.7 Reporting

The report() function can export results in multiple formats such as “rtf”, “docx”, “html”, and “pptx”. For example,

report(x=result, name = "results", format = c("html", "docx"))

creates results.html and results.docx (file extensions are added automatically). Replace "results" with any base name as needed.

3.4 Benchmark using SAS

We evaluated StepReg’s performance by comparing its outputs with those of SAS using multiple public datasets. For most cases, the information criterion values calculated at each step and the final models generated by StepReg were consistent with those from SAS across different strategies and metrics. However, it is important to note that p-value calculations for variable inclusion in Cox models may differ slightly between StepReg and SAS for variable entry. This discrepancy arises because survival::coxph() utilizes the likelihood ratio test (LRT), whereas SAS employs the score chi-square method (SAS Institute Inc. 2019b). Despite this difference, the final models remained consistent between StepReg and SAS across the test datasets (see Supplementary Table 5). These results confirm StepReg’s accuracy and reliability as a tool for stepwise regression analysis.

4 Discussion

4.1 StepReg enables comprehensive stepwise regression analyses

StepReg provides a unified and flexible framework for stepwise regression. It supports six commonly used regression families (type), four variable selection strategies (forward, backward, bidirectional, best-subset), and a broad set of information criteria (metrics). Unlike tools that run a single strategy–criterion combination at a time, StepReg lets users specify multiple strategies and criteria in one call, enabling efficient, side-by-side comparisons and reducing boilerplate code. To improve robustness, it offers data splitting via test_ratio for out-of-sample evaluation and more reliable post-selection summaries and an optional randomized forward selection mode via features_ratio to mitigate overfitting and stabilize choices.

For interpretability and diagnostics, StepReg performs basic multicollinearity checks (tolerance) and records the full selection path. The plot() function visualizes the path and selected variables; the performance() function summarizes train/test performance across all strategy-by-criterion combinations; and the report() function produces HTML/DOCX/RTF/PPTX reports for easy sharing.

To enhance accessibility, an interactive Shiny-based GUI is accessible through the companion package StepRegShiny, enabling point-and-click analyses for non-programmers.

In summary, by integrating broad model-family coverage, multiple strategies and criteria within a unified workflow, along with features such as data splitting, randomized forward selection, visualization, and automated reporting, StepReg makes stepwise regression more efficient, transparent, and accessible than existing tools.

4.2 Best practices and considerations for using StepReg

4.2.1 Scope of use

StepReg provides a flexible framework for stepwise regression suitable for both exploratory analysis and predictive modeling. Effective use depends on appropriate parameterization and careful interpretation.

4.2.2 Parameterization

• Variable selection strategy. Users should choose the strategy to match data and goals. Forward selection is efficient when starting from many candidate predictors; backward elimination is appropriate when a well-justified, estimable full model is available and the number of predictors is modest relative to the number of observations; bidirectional selection mixes additions and deletions but remains path-dependent; best-subset evaluates all subsets within the specified scope and is computationally intensive (roughly \(2^p\)).

• Information criteria. Information criteria (e.g., AIC, BIC, and SBC) generally balance fit and complexity for prediction; SL rules align with hypothesis-driven screening but warrant caution in interpretation. StepReg allows multiple strategies and criteria in a single call (strategy, metric) to enable side-by-side comparisons.

• A priori inclusion. Use include to retain domain-critical predictors irrespective of selection statistics.

4.2.3 Limitations and mitigation

Models obtained via stepwise procedures are not suitable for classical post-selection inference without adjustment; the search can induce overfitting, spurious associations, and biased estimates. StepReg provides two built-in mitigations: data splitting (test_ratio) to select on a training subset and refit/evaluate on a hold-out set, and randomized forward selection (features_ratio) to sample a subset of candidates at each step, which can stabilize choices and reduce overfitting. When formal inference is the goal, complementary approaches—penalized regression, bootstrap or post-selection inference methods, and cross-validated regularization—should be considered.

4.2.4 Performance and scalability

StepReg runs efficiently on small to moderate problems; wall time increases with the number of candidate predictors p, the complexity of the model family (e.g., Cox, negative binomial), and the number of strategy-by-criterion combinations evaluated. We omit parallel execution to avoid extra dependencies and variability, though it could reduce wall time in high-throughput settings. Tree-based and boosting methods are out of scope for this release; future versions may expose them through connectors to higher-level ecosystems (e.g., tidymodels).

4.2.5 Summary

StepReg is well suited for exploratory analysis, enabling rapid screening and ranking of candidate predictors, as well as for predictive modeling, where the priority is out-of-sample performance rather than coefficient interpretation. In the latter scenario, it provides a practical, unified workflow for selecting and comparing high-performing models. With thoughtful choices of strategy and criterion, principled inclusion of a priori predictors, and use of mitigations such as data splitting and randomized forward selection, StepReg can be applied effectively and responsibly across diverse analytic contexts. Although very large or highly complex problems may be computationally demanding, the current design balances functionality and usability while remaining readily extensible. While the current release standardizes stepwise selection for six likelihood-based families, future development will emphasize interoperability—exploring a tidymodels integration to expose alternative model engines while retaining StepReg’s standardized interface (formula syntax, information criteria, visualization, reporting).

4.3 Adoption and impact

StepReg has been publicly available on CRAN since 2019 and has been cited by over 80 research studies across diverse fields. It has been especially valuable in medicine and health sciences, where Cox and logistic regression are widely used for prognostic or predictive modeling and biomarker discovery. For example, Flammia et al. (2024) used StepReg to develop a Cox model predicting chronic kidney disease upstaging after robot-assisted partial nephrectomy, and Pigeyre et al. (2023) applied forward-selection logistic regression to identify 25 serum biomarkers for subtyping type 2 diabetes. Beyond clinical applications, StepReg has been used in other domains. Yoon et al. (2022) employed stepwise linear regression to assess how landscape factors reduce PM2.5 levels. Overall, StepReg’s comprehensive regression toolkit supports a wide range of scientific inquiries, making it a practical resource for advanced analyses across disciplines.

5 Acknowledgments

We thank the Department of Molecular, Cell, and Cancer Biology (MCCB) at UMass Chan Medical School for internal funding and Dr. Haibo Liu for his helpful suggestions during manuscript revision.

5.1 Supplementary materials

Supplementary materials are available in addition to this article. It can be downloaded at RJ-2026-005.zip

5.2 CRAN packages used

MASS, leaps, bestglm, olsrr, blorr, My.stepwise, StepReg, StepRegShiny, survival, ggplot2, flextable, shiny, cowplot, tidymodels

5.3 CRAN Task Views implied by cited packages

ChemPhys, ClinicalTrials, Databases, Distributions, Econometrics, Environmetrics, MachineLearning, MixedModels, NetworkAnalysis, NumericalMathematics, Phylogenetics, Psychometrics, ReproducibleResearch, Robust, Spatial, Survival, TeachingStatistics, WebTechnologies