Stepwise Regression

Stepwise Regression to Select Appropriate Models

stepwiselm creates a linear model and automatically adds to or trims the model. To create a small model, start from a constant model. To create a large model, start with a model containing many terms. A large model usually has lower error as measured by the fit to the original data, but might not have any advantage in predicting new data.

stepwiselm can use all the name-value options from fitlm, with additional options relating to the starting and bounding models. In particular:

• For a small model, start with the default lower bounding model: 'constant' (a model that has no predictor terms).

• The default upper bounding model has linear terms and interaction terms (products of pairs of predictors). For an upper bounding model that also includes squared terms, set the Upper name-value pair to 'quadratic'.

Compare large and small stepwise models

This example shows how to compare models that stepwiselm returns starting from a constant model and starting from a full interaction model.

Load the carbig data and create a table from some of the data.

load carbig
tbl = table(Acceleration,Displacement,Horsepower,Weight,MPG);

Create a mileage model stepwise starting from the constant model.

mdl1 = stepwiselm(tbl,'constant','ResponseVar','MPG')

1. Adding Weight, FStat = 888.8507, pValue = 2.9728e-103
2. Adding Horsepower, FStat = 3.8217, pValue = 0.00049608
3. Adding Horsepower:Weight, FStat = 64.8709, pValue = 9.93362e-15

mdl1 = 
Linear regression model:
    MPG ~ 1 + Horsepower*Weight

Estimated Coefficients:
                          Estimate         SE         tStat       pValue  
                         __________    __________    _______    __________

    (Intercept)              63.558        2.3429     27.127    1.2343e-91
    Horsepower             -0.25084      0.027279    -9.1952    2.3226e-18
    Weight                -0.010772    0.00077381    -13.921    5.1372e-36
    Horsepower:Weight    5.3554e-05    6.6491e-06     8.0542    9.9336e-15


Number of observations: 392, Error degrees of freedom: 388
Root Mean Squared Error: 3.93
R-squared: 0.748,  Adjusted R-Squared: 0.746
F-statistic vs. constant model: 385, p-value = 7.26e-116

Create a mileage model stepwise starting from the full interaction model.

mdl2 = stepwiselm(tbl,'interactions','ResponseVar','MPG')

1. Removing Acceleration:Displacement, FStat = 0.024186, pValue = 0.8765
2. Removing Displacement:Weight, FStat = 0.33103, pValue = 0.56539
3. Removing Acceleration:Horsepower, FStat = 1.7334, pValue = 0.18876
4. Removing Acceleration:Weight, FStat = 0.93269, pValue = 0.33477
5. Removing Horsepower:Weight, FStat = 0.64486, pValue = 0.42245

mdl2 = 
Linear regression model:
    MPG ~ 1 + Acceleration + Weight + Displacement*Horsepower

Estimated Coefficients:
                                Estimate         SE         tStat       pValue  
                               __________    __________    _______    __________

    (Intercept)                    61.285        2.8052     21.847    1.8593e-69
    Acceleration                 -0.34401       0.11862       -2.9     0.0039445
    Displacement                -0.081198      0.010071    -8.0623    9.5014e-15
    Horsepower                   -0.24313      0.026068    -9.3265    8.6556e-19
    Weight                     -0.0014367    0.00084041    -1.7095      0.088166
    Displacement:Horsepower    0.00054236    5.7987e-05     9.3531    7.0527e-19


Number of observations: 392, Error degrees of freedom: 386
Root Mean Squared Error: 3.84
R-squared: 0.761,  Adjusted R-Squared: 0.758
F-statistic vs. constant model: 246, p-value = 1.32e-117

Notice that:

• mdl1 has four coefficients (the Estimate column), and mdl2 has six coefficients.

• The adjusted R-squared of mdl1 is 0.746, which is slightly less (worse) than that of mdl2, 0.758.

Create a mileage model stepwise with a full quadratic model as the upper bound, starting from the full quadratic model:

mdl3 = stepwiselm(tbl,'quadratic','ResponseVar','MPG','Upper','quadratic');

1. Removing Acceleration:Horsepower, FStat = 0.075209, pValue = 0.78405
2. Removing Acceleration:Weight, FStat = 0.072756, pValue = 0.78751
3. Removing Horsepower:Weight, FStat = 0.12569, pValue = 0.72314
4. Removing Weight^2, FStat = 1.194, pValue = 0.27521
5. Removing Displacement:Weight, FStat = 1.2839, pValue = 0.25789
6. Removing Displacement^2, FStat = 2.069, pValue = 0.15114
7. Removing Horsepower^2, FStat = 0.74063, pValue = 0.39

Compare the three model complexities by examining their formulas.

mdl1.Formula

ans = 
MPG ~ 1 + Horsepower*Weight

mdl2.Formula

ans = 
MPG ~ 1 + Acceleration + Weight + Displacement*Horsepower

mdl3.Formula

ans = 
MPG ~ 1 + Weight + Acceleration*Displacement + Displacement*Horsepower
         + Acceleration^2

The adjusted $$R^2$$ values improve slightly as the models become more complex:

RSquared = [mdl1.Rsquared.Adjusted, ...
    mdl2.Rsquared.Adjusted, mdl3.Rsquared.Adjusted]

RSquared = 1×3

    0.7465    0.7580    0.7599

Compare residual plots of the three models.

subplot(3,1,1)
plotResiduals(mdl1)
subplot(3,1,2)
plotResiduals(mdl2)
subplot(3,1,3)
plotResiduals(mdl3)

The models have similar residuals. It is not clear which fits the data better.

Interestingly, the more complex models have larger maximum deviations of the residuals:

Rrange1 = [min(mdl1.Residuals.Raw),max(mdl1.Residuals.Raw)];
Rrange2 = [min(mdl2.Residuals.Raw),max(mdl2.Residuals.Raw)];
Rrange3 = [min(mdl3.Residuals.Raw),max(mdl3.Residuals.Raw)];
Rranges = [Rrange1;Rrange2;Rrange3]

Rranges = 3×2

  -10.7725   14.7314
  -11.4407   16.7562
  -12.2723   16.7927