PPOL670 | Interpretable Machine Learning | Week 12

Walkthrough on some approaches to model interpretation

Overview

In this notebook, we’ll leverage the data and best fit model that we learned from the walkthrough last week (Week 11). The aim is to extract insights from the machine learning model to better understand who is covered and who is not.

require(tidyverse)
require(caret) # for machine learning
require(recipes) # For preprocessing your data
require(vip) # For variable importance
require(pdp) # pdp and ice plots

# For parallelization (to run the models across many cores) -- speeds up computation!
# install.packages("doMC")
doMC::registerDoMC()

Data

The following data contains information regarding whether or not someone has health coverage. The outcome of interest is whether of not an individual has healthcare coverage or not. The available predictive features capture socio-economic and descriptive factors.

set.seed(1988)

# Again use a small random sample of the data so that it runs quicker
dat = suppressMessages(read_csv("health-coverage.csv")) %>%
  
  # Convert education into an ordered category.
  mutate(
    
    educ = case_when(
      educ == 'Less than HS' ~ 0,
      educ == 'HS Degree' ~ 1,
      educ == 'Undergraduate Degree' ~ 2,
      educ == 'Graduate Degree' ~ 3),
    
    race = case_when(
      race == 'White' ~ 'white',
      race == 'Black' ~ 'black',
      race == 'Asian' ~ 'asian',
      T ~ 'other'),
    
    race = factor(race,level=c("white",'black','asian','other')),
    
    mar = factor(mar,levels= c('Never Married',
                               'Divorced','Married',
                               'Separated','Widowed'))
    
    ) %>% 
  
  # Convert all remaining factor variables into character variables
  mutate_if(is.character,as.factor) %>% 
  
  # Make sure that the category you want to predict is the first category in
  # your factor variable
  mutate(coverage = factor(coverage,
                           levels = c("Coverage","No_Coverage"))) %>% 
  
  sample_n(5000)

head(dat) # Peek at the data just to make sure everything was read in correctly. 

Split the Sample: Training and test data

As before, let’s split the sample up into a training and test dataset. We’ll completely hold off on viewing the test data. We’ll leverage the test data to explore some of the insights we learn from the training data. Let’s use rsample to split the data.

set.seed(123)
# install.packages("rsample") # If you don't have the pack
split <- rsample::initial_split(dat,prop = .8,)
train_data <- rsample::training(split) 
test_data  <- rsample::testing(split) 


# Breakdown
nrow(train_data)/nrow(dat) 

[1] 0.8002

nrow(test_data)/nrow(dat)

[1] 0.1998

Data Pre-processing

In the Week 11 Walkthrough, we explore the data and outlined a preprocessing strategy. Let’s re-implement it here.

rcp <- 
  recipe(coverage ~  .,train_data) %>% 
  step_dummy(all_nominal(),-coverage) %>% #
  step_log(wage,offset = 1) %>% # Log the skewed wage variable
  step_range(all_numeric()) %>%  # Normalize scale
  prep()


# Apply the recipe to the training and test data
train_data2 <- bake(rcp,train_data)
test_data2 <- bake(rcp,test_data) 

Cross-validation

Let’s use the k-fold cross-validation method with 5 folds in the data to tune the model. We also need to specify that we’re dealing with a classification problem: recall we can do this by adding summaryFunction = twoClassSummary and classProbs = TRUE to the cross-validation function.

set.seed(1988) # set a seed for replication purposes 

# Partition the data into 5 equal folds
folds <- createFolds(train_data2$coverage, k = 5) 

# Cross validation object
control_conditions <- 
  trainControl(method='cv', # K-fold cross validation
               summaryFunction = twoClassSummary, # Two classes 
               classProbs = TRUE, # Need this b/c it's a classification problem
               index = folds # The indices for our folds
  )

Fit Model

In the Week 11 Walkthrough, we found that the random forest model had the best performance on predicting the outcome. Let’s re-run that model here.

mod_rf <-
  train(coverage ~ ., # Equation (outcome and everything else)
        data=train_data2, # Training data 
        method = "ranger", # random forest (ranger is much faster than rf)
        metric = "ROC", # area under the curve
        trControl = control_conditions
  )

Plot the model output. When we plot the performance of the above model, we see that we get a fairly high AUC (.83) on the held out data.

# AUC on the Training data 
prob_pred <- predict(mod_rf,type="prob")$Coverage
auc_rf <- Metrics::auc(train_data2$coverage=="Coverage",prob_pred)
auc_rf

[1] 0.8838971

# Accuracy on the Training data 
acc_rf <- Metrics::accuracy(train_data2$coverage,predict(mod_rf))
acc_rf

[1] 0.7853037

Model Interpretation

Variable Importance

Let’s now extract some substantive insights from the model using variable importance approach learned in class. Which variables does the model rely on most to make it’s predictions?

vi_plot<- 
  vip(mod_rf, # Machine learning model
      train = train_data2, # Training data 
      method="permute", # permuted importance
      nsim = 10, # number of times to impute
      geom = "boxplot", # Type of plot 
      target = "coverage", # outcome
      reference_class = "Coverage", # what class are you predicting
      metric = "accuracy",
      pred_wrapper = predict)

# Plot VIP
vi_plot

The variable importance plot shows:

• age, wage, and educ, is highly important to the models predictive accuracy
• Things that matter less, but still matter to whether you’re covered or not
  • If someone is married,
  • Citizenship status.

Looking at the variable importance alone, we could refine our original model to only include those features that matter. Let’s see if doing so boosts performance.

# Re-run the random forest, this time with most important variable only. 
mod_rf_selected <-
  train(coverage ~ age + wage + educ + mar_Married + cit_Non.citizen, 
        data=train_data2, # Training data 
        method = "ranger", # random forest (ranger is much faster than rf)
        metric = "ROC", # area under the curve
        trControl = control_conditions
  )

Does the model do better, worse, or about the same?

# AUC on the Training data 
prob_pred <- predict(mod_rf_selected,type="prob")$Coverage
auc_rf_sel <- Metrics::auc(train_data2$coverage=="Coverage",prob_pred)

cat("AUC of the selected model:",auc_rf_sel,"\n",
    "AUC of the full model:",auc_rf,"\n",
    "Difference:", auc_rf_sel-auc_rf,"\n") # Difference: a 4.2% increase in AUC

AUC of the selected model: 0.9201622 
 AUC of the full model: 0.8838971 
 Difference: 0.0362651 

# Accuracy on the Training data 
acc_rf_sel <- Metrics::accuracy(train_data2$coverage,predict(mod_rf_selected))

cat("Accuracy of the selected model:",acc_rf_sel,"\n",
    "Accuracy of the full model:",acc_rf,"\n",
    "Difference:", acc_rf_sel-acc_rf,"\n") # Difference: a 3.4% increase in accuracy

Accuracy of the selected model: 0.8322919 
 Accuracy of the full model: 0.7853037 
 Difference: 0.04698825 

Clearly the model performs better (in this case) when we only include the features that matter most.

Looking into which variables are important to the model helped us refine the model and the features worth using.

Marginal Effects

Partial Dependencies

It’s useful to know that age and wage are important to the model, but how do they relate to whether an individual is covered or not in the data? Let’s explore a partial dependency plot to unpack the average marginal effect for age and wage.

# PDP for age
age_pdp <- partial(mod_rf_selected, pred.var = "age", plot = TRUE,prob=T,
                   grid.resolution = 20, # choosing less points makes it run quicker
                   plot.engine = "ggplot2")

# PDP for wage
wage_pdp <- partial(mod_rf_selected, pred.var = "wage", plot = TRUE,prob=T,
                    grid.resolution = 20, # choosing less points makes it run quicker
                   plot.engine = "ggplot2")

# PDP for education level
educ_pdp <- partial(mod_rf_selected, pred.var = "educ", plot = TRUE,prob=T,
                    grid.resolution = 4,plot.engine = "ggplot2")

# PDP for Married (note that this is a dummy var, so only can take on one of two values)
mar_pdp <- partial(mod_rf_selected, pred.var = "mar_Married", plot = TRUE,prob=T,
                   plot.engine = "ggplot2")

# PDP for Not a citizen (note that this is a dummy var, so only can take on one of two values)
cit_pdp <- partial(mod_rf_selected, pred.var = "cit_Non.citizen", plot = TRUE,prob=T,
                   plot.engine = "ggplot2")

# GridExtra allows us to stack two ggplot objects into a single plot
gridExtra::grid.arrange(age_pdp,wage_pdp,educ_pdp,mar_pdp,cit_pdp)

Looking above at the partial dependencies for all the variables in the model, we can see some interesting things:

• For age, there appears to exist a threshold. Once you get past a certain age (50-sh), you’re incredibly more likely to be covered.

• As a side note, notice the bump at the lower end of the age distribution? Minimum age in the data is 16. This predicted increase in the probability of coverage likely reflects children who are covered by their parents.

• For wage, it looks like as you get wealthier, you’re more likely to be covered.

• The more educated one is, the higher the likelihood that they’ll be covered.

• For the dummy variables: being married corresponds with an increase in the probability of being covered, whereas not being a citizen corresponds with a decrease in the likelihood of being covered.

Individual Conditional Expectations

One additional question we might want to ask is whether there is are important interactions. One way we can explore this is if we look at the individual expectations for each observation in the data as we manipulate the values of the target variable. If the ICE curves diverge (have different trajectories), that can be evidence that there is a potential interaction.

A few things to note here:

• It looks like there is a potential interaction in wage. Notice how some of the curves have a very different trend toward the end of the distribution? The reason for this is likely due to some interaction.

• The idea that only old people are covered (our conclusion after looking at the PDP plot) isn’t entirely accurate. In fact, there are a number of respondents in the data that are predicted to have a higher probability of coverage at a young age; however, these observations are canceled out by other young respondents who are less likely to be covered.

• Likewise it looks like there is heterogeneity in the effect of education on the likelihood of being covered.

• Finally, for the dummy variables, note how the impact of being married or not a citizen impacts different respondents differently in the model.

Let’s now explore the potential interaction between wage and age, wage and educ, and educ and wage. Is it the combination of being wealthy and old that corresponds with coverage, or does wealth make younger people more willing to purchase coverage?

From the above plot, we can indeed see that there is an interaction between age and wage. Specifically,

• We still note the threshold effect past .65 on the age scale (that is, all older people are more likely to be covered).
• However, what this plot demonstrates is that wage mediates age. Even if you’re young, if you’re making a high salary you’re more likely to be covered. Chances are that a high wage is correlated with employer-provided health coverage.
• Finally, there appears to be a major region in which the model predicts a lower probability of coverage. Specifically for young people who don’t have a high paying wage.

W/r/t education and wage: the probability of coverage is highest for educated people who are on the high end of the wage distribution. However, the effect of education is much less (even for highly educated people) if they are lower on the wage distribution.

W/r/t education and age, it doesn’t look like there is much of an interaction.

Global Surrogate Model

Tree

The random forest model does a fantastic job at predicting the likelihood of coverage, but it’s difficult to interpret. In the above sections we explored variable importance and the marginal effects of dependencies of the variables on the predicted probabilities to gain some insight into how the model is using the data to make its predictions. Let’s now take this a step further by building a global surrogate model. The aim here is to build a more interpretable model on top of our black box model. That is, we can train an interpretable model to approximate the predictions of the random forest model. We can then interpret this more interpretable model.

train_data3 <- 
  train_data2 %>% 
  mutate(
    # Probability predictions from the random forest model
    coverage_probs = predict(mod_rf_selected,type = "prob")$Coverage
  )  


surrogate_tree <-
  
  # Decision tree model (directly us the model rather than use the implementation in caret)
  rpart::rpart(
    # Main selected model. The outcome is now the predicted probabilities from
    # the RF model
    coverage_probs ~ age + wage + educ + mar_Married + cit_Non.citizen,
    
    # Data is being passed by the pipe
    data = train_data3,
    
    # Note that we can control the depth of the tree
    # a deeper tree increase fit but reduces interpretability
    control = rpart::rpart.control(maxdepth = 4)
  )

Let’s now assess the model fit.

tibble(truth = train_data3$coverage_probs,
       estimate = predict(surrogate_tree)) %>% 
  rsq(truth,estimate)

An R-Squared of .78 isn’t bad. It indicates that the surrogate model does a good job approximating the black box model (i.e. the random forest model).

Let’s now look at the tree.

rattle::fancyRpartPlot(surrogate_tree,sub="",type=1)

The surrogate model provides a picture that chimes well with our other insights, but new ones clearly stand out.

• As we previously noted, there is an interaction between age and wage. Though there are two notable interactions given ones wage. Assuming one is less than 65 years of age (i.e. .64 on our age scale):

  • If one doesn’t make the upper quantile of the wage distribution, then coverage depends on whether they’re a citizen or not (though note the probabilities are low for both splits).
  • If one does make a high wage and is married he/she is more likely to be covered.

• If one is more than 65 years of age, coverage mainly depends on whether he/she is a citizen or not. If you’re above 65 and a citizen, you have a 95% probability of being covered.

NOTE: that if we built a deeper model, a richer story emerges. But that story doesn’t differ much from the one we’ve already been piecing together.

surrogate_tree2 <-
  rpart::rpart(coverage_probs ~ age + wage + educ + mar_Married + cit_Non.citizen,
    data = train_data3,
    control = rpart::rpart.control(maxdepth = 7)
  )
rattle::fancyRpartPlot(surrogate_tree2,sub="",type=1)

Linear Model

We can use the insights above to build out a linear model (actually, a logistic model b/c the outcome is binary) to help with interpretability.

linear_surrogate_model <- 
  train_data2 %>% 
  mutate(coverage_dummy = 1*(coverage == "Coverage")) %>% 
  glm(coverage_dummy ~ 
        age*wage + 
        wage*educ +
        mar_Married + 
        cit_Non.citizen,
      data=.,family = binomial("logit"))

# Print the regression results
stargazer::stargazer(linear_surrogate_model,type = "text")

length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changed

=============================================
                      Dependent variable:    
                  ---------------------------
                        coverage_dummy       
---------------------------------------------
age                        3.505***          
                            (0.230)          
                                             
wage                        0.418*           
                            (0.243)          
                                             
educ                       0.821***          
                            (0.247)          
                                             
mar_Married                0.690***          
                            (0.077)          
                                             
cit_Non.citizen            -1.551***         
                            (0.148)          
                                             
age:wage                   -1.925***         
                            (0.475)          
                                             
wage:educ                  2.345***          
                            (0.416)          
                                             
Constant                   -1.988***         
                            (0.128)          
                                             
---------------------------------------------
Observations                 4,001           
Log Likelihood            -2,230.519         
Akaike Inf. Crit.          4,477.038         
=============================================
Note:             *p<0.1; **p<0.05; ***p<0.01

We could now test these theoretical insights on new, previously unseen data… i.e. the test data!

How well does the theoretical model predict future data? Pretty well, which indicates that our insights likely generalize to the test data.

Metrics::auc(actual = 1*(test_data2$coverage=="Coverage"),
             predicted = predict(linear_surrogate_model,
                                 newdata = test_data2))

[1] 0.7827528

Theory Development

We’ve learned a lot from the training data and model about the factors that increase the likelihood of coverage. We’ve generated a number of interesting hypotheses:

• Retirement age US citizens are more likely to have healthcare coverage than non-retirement age citizens.
• Young married individuals making a high salary are more likely to be covered than their non-married counterparts.
• Coverage largely depends on an individuals income unless he/she is above 65 years of age.