Variable importance logistic and random forest

Question

I have created variable importance plots using varImp in R for both a logistic and random forest model. I want to compare how the logistic and random forest differ in the variables they find important. Of course, they do this in a different way: logistic takes the absolute value of the t-statistic and the random forest the mean decrease in Gini.

Now, when I plot the variable importance plots for the logistic and the random forest, I find that the logistic and the random forest model handle factorial variables in a different way, whilst the random forest model takes the total group, the logistic regression takes one of the possible factor outcomes.

For example, I have a group US states with 50 factors, the logistic regression would take New York as an important factor, whilst the random forest model takes the group US states as the important factor.

How can I solve this for the logistic regression? Could I just add all the values for the variable importance to get the total group importance?

Answer 1

It is most likely not a good idea.

1. If you have many coefficients that are not very useful, i.e low T statistics, but adding up 50 of them might give you something huge... which just doesn't make sense.

2. T-statistic doesn't take into account the explained variance. Worst scenario, one of one of your categories end up in a sweet spot, it has low number of observations and by chance a small standard error, a huge t-statistic. Adding this up to your term inflates the importance.

We can use an example below:

library(survival)
library(randomForest)
library(caret)

da = survival::diabetic[,-1]
# make age categories
da$age = cut(diabetic$age,10)
da$status = factor(da$status) 

glm_mdl = glm(status ~ .,data=da,family=binomial)

rf_mdl = randomForest(status ~ .,data=da)

If we look at the summary of glm, seems like age has an effect, but if you sum up the tstat for all age, you end up with something huge:

    Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     1.063128   1.101749   0.965   0.3346    
laserargon     -0.048476   1.151578  -0.042   0.9664    
age(6.7,12.4]   0.964098   0.501488   1.922   0.0545 .  
age(12.4,18.1]  0.500876   0.525536   0.953   0.3406    
age(18.1,23.8]  2.191287   1.144998   1.914   0.0556 .  
age(23.8,29.5]  0.945382   1.333947   0.709   0.4785    
age(29.5,35.2]  0.849438   1.361294   0.624   0.5326    
age(35.2,40.9]  1.497774   1.425724   1.051   0.2935    
age(40.9,46.6]  0.545537   1.312921   0.416   0.6778    
age(46.6,52.3]  1.565862   1.385946   1.130   0.2586    
age(52.3,58.1]  0.945929   1.500791   0.630   0.5285    
eyeright        0.484579   0.293928   1.649   0.0992 .  
trt            -1.098955   0.295500  -3.719   0.0002 ***
risk            0.097595   0.103325   0.945   0.3449    
time           -0.094334   0.009613  -9.814   <2e-16 ***

We check the change in deviance (how good it is at reducing prediction error), it's actually quite little:

anova(glm_mdl)

          Df Deviance Resid. Df Resid. Dev
NULL                    393     528.15
laser  1    0.317       392     527.84
age    9    3.716       383     524.12
eye    1    3.110       382     521.01
trt    1   26.404       381     494.61
risk   1    5.107       380     489.50
time   1  179.399       379     310.10

If you like the variable importance to reflect how useful the variable is at predicting correctly, I think a fairer comparison might be change in deviance, so we can try something like:

v_glm = anova(glm_mdl)[-1,2,drop=FALSE]
v_glm = v_glm[order(v_glm[,1]),drop=FALSE,]
v_glm[,1] = 100*v_glm[,1]/max(v_glm[,1])

v_rf = as.matrix(varImp(rf_mdl))
v_rf =  v_rf[order(v_rf),]

And we get the estimate if we sum up the importance as you raised:

v_glm_sum = as.matrix(varImp(glm_mdl))
age_row = grepl("age",rownames(v_glm_sum))
v_glm_sum = rbind(age=sum(v_glm_sum[age_row,]),v_glm_sum[!age_row,drop=FALSE,])
v_glm_sum =  v_glm_sum[order(v_glm_sum),]

Now plot and we can see the sum of the importance of categories will be inflated, so most likely the deviance is something closer, for comparison:

par(mfrow=c(1,3))
barplot(t(v_rf),horiz=TRUE,main="rf",las=2)
barplot(t(v_glm),horiz=TRUE,main="glm_deviance",las=2)
barplot(t(v_glm_sum),horiz=TRUE,main="glm_sum_scores",las=2)