2
$\begingroup$

I am facing the following piece of legacy code, and there is no chance to speak with its author unfortunately.

Model1_LM <- lm(result ~ ., data = data)
fit1_LM <- stepAIC(Model1_LM, direction = 'backward')

Here data is a dataframe which contains features and the target to be predicted, the latter is the "result" column of the data. I have not used R, but I do know say Python and have some experience with learning, linear and non-linear. I am trying to understand what does stepAIC do here.

In my understanding, and judging from the documentation, stepAIC chooses the best model from the class according to the Akaike information criterion (AIC). As far as I understood, the class of models in my case are linear models that try to predict the result column based on all other ones. Now, is this class parametrized in R: in case that were say neural nets, I would think that the class would be parametrized by number of hidden nodes, but in case of a linear model in R I am no sure how the selection is done. Also, I am surprised that for this selection no separation in training and validation is needed: a single dataframe is provided. Perhaps, the split is done under the hood using some default method, but I have not found information about this in the docs.

Bottom line: I would be happy if someone told me what exactly is happening under the hood of stepAIC in the code above. Thanks.

$\endgroup$
8
$\begingroup$

Here's a simple example that illuminates how the stepAIC() function from the MASS package works.

# consider the built-in mtcars data set; select variables of interest 
data <- mtcars[,c("mpg","disp","hp","wt")]

# examine the first 6 records of the data 
head(data)

The first 6 records in this example data set are as follows:

> head(data)
                   mpg disp  hp    wt
Mazda RX4         21.0  160 110 2.620
Mazda RX4 Wag     21.0  160 110 2.875
Datsun 710        22.8  108  93 2.320
Hornet 4 Drive    21.4  258 110 3.215
Hornet Sportabout 18.7  360 175 3.440
Valiant           18.1  225 105 3.460

The example data set includes a total of 32 records, one for each type of car, on the variables:

mpg  Miles/(US) gallon  
disp  Displacement (cu.in.)  
hp  Gross horsepower  
wt  Weight (1000 lbs)  

We'll use the variable mpg as the response variable and the variables disp, hp and wt as the predictor variables.

The stepAIC() function performs backward model selection by starting from a "maximal" model, which is then trimmed down.

The "maximal" model in our example is given by the commands:

Model1_LM <- lm(mpg ~ ., data = data)

summary(Model1_LM)

The summary associated with the maximal model is as follows:

>     summary(Model1_LM)

Call:
lm(formula = mpg ~ ., data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-3.891 -1.640 -0.172  1.061  5.861 

Coefficients:
             Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 37.105505   2.110815  17.579 < 0.0000000000000002 ***
disp        -0.000937   0.010350  -0.091              0.92851    
hp          -0.031157   0.011436  -2.724              0.01097 *  
wt          -3.800891   1.066191  -3.565              0.00133 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.639 on 28 degrees of freedom
Multiple R-squared:  0.8268,    Adjusted R-squared:  0.8083 
F-statistic: 44.57 on 3 and 28 DF,  p-value: 0.0000000000865

The "maximal" model is a linear regression model which assumes independent model errors and includes only main effects for the predictor variables.

The "maximal" model is used as input for the stepAIC() function:

require(MASS)

fit1_LM <- stepAIC(Model1_LM, direction = 'backward')

The backward model selection produced by this function is:

Start:  AIC=65.83                    
mpg ~ disp + hp + wt

       Df Sum of Sq    RSS    AIC
- disp  1     0.057 195.05 63.840    <- "maximal" model without disp
<none>              194.99 65.831    <- "maximal" model
- hp    1    51.692 246.68 71.356    <- "maximal" model without hp
- wt    1    88.503 283.49 75.806    <- "maximal" model without wt

Step:  AIC=63.84
mpg ~ hp + wt

       Df Sum of Sq    RSS    AIC
<none>              195.05 63.840  <- "improved" model 
- hp    1    83.274 278.32 73.217  <- "improved" model without hp 
- wt    1   252.627 447.67 88.427  <- "improved" model without wt 

The first two rows of this output show the AIC and model formula for the "maximal" model.

The next portion of the output shows what happens to the AIC value of the reduced models obtained by eliminating, in turns, each of the predictors from the "maximal" model. If the elimination of a predictor improves the "maximal" model, the AIC should become smaller.

It turns out that when eliminate the predictor hp from the "maximal" model, we achieve a smaller AIC (i.e., AIC = 63.840), than the one corresponding to the "maximal" model (i.e., AIC = 65.831). This suggests that the "maximal" model can be improved by simply discarding hp. The "improved" model includes only hp and wt as predictor variables and has an AIC of 63.840.

We now repeat the same process as the one described in the paragraph above - we start out with the "improved" model, eliminate each of the remaining predictors hp and wt in turns from the model and monitor how that impacts the AIC. If the AIC becomes smaller, that means that the eliminated predictor does not contribute much to the model.

In our example, we cannot reduce the "improved" model's AIC by eliminating either hp or wt from the model. In fact, when we eliminate each of these predictor variables from the model, the AIC increases, suggesting the elimination has a detrimental effect on the model.

So the "final" model is the one posted by R when we type the following command:

fit1_LM 

In other words, the "final" model is:

> fit1_LM 

Call:
lm(formula = mpg ~ hp + wt, data = data)

Coefficients:
(Intercept)           hp           wt  
   37.22727     -0.03177     -3.87783  

This model includes only hp and wt as predictor variables. We can now examine the summary of the fit of this "final" model to our data:

summary(fit1_LM)

The summary looks like this:

> summary(fit1_LM)

Call:
lm(formula = mpg ~ hp + wt, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-3.941 -1.600 -0.182  1.050  5.854 

Coefficients:
        Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 37.22727    1.59879  23.285 < 0.0000000000000002 ***
hp          -0.03177    0.00903  -3.519              0.00145 ** 
wt          -3.87783    0.63273  -6.129           0.00000112 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.593 on 29 degrees of freedom
Multiple R-squared:  0.8268,    Adjusted R-squared:  0.8148 
F-statistic: 69.21 on 2 and 29 DF,  p-value: 0.000000000009109

The summary shows that both hp and wt have significant effects on mpg. Among cars with the same wt, hp is negatively associated with mpg. Among cars with the same hp, wt is negatively associated with mpg.

In this example, stepAIC() was applied to the entire data set. Note that AIC is recommended as a criterion for model selection in those situations where the model will be used for predictive purposes.

$\endgroup$
  • $\begingroup$ Actually, it can go backward, forward, or both depending on the argument you give it: rdocumentation.org/packages/MASS/versions/7.3-51.4/topics/… It's not obvious how 'both' would work. Would have to visit the source code to find out. Haven't tested, but I would guess 'both' might give different results from backwards or forwards. $\endgroup$ – wordsforthewise Sep 26 at 20:39
2
$\begingroup$

It is doing model selection based on Akaike's Information Criterion, which is calculated as AIC = 2k - 2lnL, where k is the number of parameters estimated by the model and lnL is the log likelihood of the data given the model. Basically, model selection with AIC attempts to select the model that best explains the data (highest likelihood), while still not fitting too many parameters.

As a metric, AIC only makes sense relative to other values; its absolute value has no meaning. So in the procedure in your code stepAIC() is starting at the most complex model (because direction = "backward"), and sequentially removing terms in an effort to lower the AIC. You can think of this as building a group of models from the possible combinations of predictors in your dataset, and determining which has the lowest AIC (FYI, you can also choose direction = "forward" which starts simple and adds terms, or direction = "both" which trials both adding and removing terms).

When the stepAIC() finds a model where it cannot lower the AIC any further by removing terms (if direction = "backward"), then it selects this model as the best model.

There are arguments as to the relative merits of this approach, and of using different selection metrics (AICc, BIC, etc.). If you search Cross Validated for 'model selection' there will be an abundance of results.

EDIT: Oh yeah, AFAIK, there is no separation of training or validation etc. This basic approach is agnostic to Cross Validation.

$\endgroup$
  • $\begingroup$ an ironic EDIT. So if I understand correctly, here we use a linear model (just usual linear regression, no regularization), and try to remove features one by one (backwards in stepAIC) - in the order of columns in the dataframe. Once AIC cannot be improved by removal of the feature, the model selection stops. Does it do several loops over the features, or just a single one? Say in the first loop the 5th feature was removed, but in the second one 3rd feature gets removed etc. $\endgroup$ – Ilya May 22 '18 at 21:23
  • $\begingroup$ @Ilya I am not certain, but I believe that it trials removing every term at each step and chooses the step which leads to the greatest reduction in AIC. $\endgroup$ – NatWH May 22 '18 at 21:24
  • $\begingroup$ thanks, that gives me some understanding! unless someone else will have a better info regarding the stepAIC in R particularly, I'll accept your answer $\endgroup$ – Ilya May 22 '18 at 21:26
  • 1
    $\begingroup$ Thanks @IsabellaGhement, your answer is much better. Have upvoted. $\endgroup$ – NatWH May 22 '18 at 21:40
  • 1
    $\begingroup$ Thanks, @NatWH - I thought your answer was much better, as it gave the big picture overview. Between our two answers, Ilya should be covered. $\endgroup$ – Isabella Ghement May 22 '18 at 21:46

Not the answer you're looking for? Browse other questions tagged or ask your own question.