Machine Learning Workflow: Classification Trees

Econ 425T

Display system information for reproducibility.

Python

import IPython
print(IPython.sys_info())

{'commit_hash': 'add5877a4',
 'commit_source': 'installation',
 'default_encoding': 'utf-8',
 'ipython_path': '/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/IPython',
 'ipython_version': '8.8.0',
 'os_name': 'posix',
 'platform': 'macOS-10.16-x86_64-i386-64bit',
 'sys_executable': '/Library/Frameworks/Python.framework/Versions/3.10/bin/python3',
 'sys_platform': 'darwin',
 'sys_version': '3.10.9 (v3.10.9:1dd9be6584, Dec  6 2022, 14:37:36) [Clang '
                '13.0.0 (clang-1300.0.29.30)]'}

1 Overview

We illustrate the typical machine learning workflow for regression trees using the Heart data set from R ISLR2 package.

1. Initial splitting to test and non-test sets.

2. Pre-processing of data: not much is needed for regression trees.

3. Tune the cost complexity pruning hyper-parameter(s) using 10-fold cross-validation (CV) on the non-test data.

4. Choose the best model by CV and refit it on the whole non-test data.

5. Final prediction on the test data.

2 Heart data

The goal is to predict the binary outcome AHD (Yes or No) of patients.

Python

# Load the pandas library
import pandas as pd
# Load numpy for array manipulation
import numpy as np
# Load seaborn plotting library
import seaborn as sns
import matplotlib.pyplot as plt

# Set font sizes in plots
sns.set(font_scale = 1.2)
# Display all columns
pd.set_option('display.max_columns', None)

Heart = pd.read_csv("../data/Heart.csv")
Heart

     Unnamed: 0  Age  Sex     ChestPain  RestBP  Chol  Fbs  RestECG  MaxHR  \
0             1   63    1       typical     145   233    1        2    150   
1             2   67    1  asymptomatic     160   286    0        2    108   
2             3   67    1  asymptomatic     120   229    0        2    129   
3             4   37    1    nonanginal     130   250    0        0    187   
4             5   41    0    nontypical     130   204    0        2    172   
..          ...  ...  ...           ...     ...   ...  ...      ...    ...   
298         299   45    1       typical     110   264    0        0    132   
299         300   68    1  asymptomatic     144   193    1        0    141   
300         301   57    1  asymptomatic     130   131    0        0    115   
301         302   57    0    nontypical     130   236    0        2    174   
302         303   38    1    nonanginal     138   175    0        0    173   

     ExAng  Oldpeak  Slope   Ca        Thal  AHD  
0        0      2.3      3  0.0       fixed   No  
1        1      1.5      2  3.0      normal  Yes  
2        1      2.6      2  2.0  reversable  Yes  
3        0      3.5      3  0.0      normal   No  
4        0      1.4      1  0.0      normal   No  
..     ...      ...    ...  ...         ...  ...  
298      0      1.2      2  0.0  reversable  Yes  
299      0      3.4      2  2.0  reversable  Yes  
300      1      1.2      2  1.0  reversable  Yes  
301      0      0.0      2  1.0      normal  Yes  
302      0      0.0      1  NaN      normal   No  

[303 rows x 15 columns]

# Numerical summaries
Heart.describe(include = 'all')

        Unnamed: 0         Age         Sex     ChestPain      RestBP  \
count   303.000000  303.000000  303.000000           303  303.000000   
unique         NaN         NaN         NaN             4         NaN   
top            NaN         NaN         NaN  asymptomatic         NaN   
freq           NaN         NaN         NaN           144         NaN   
mean    152.000000   54.438944    0.679868           NaN  131.689769   
std      87.612784    9.038662    0.467299           NaN   17.599748   
min       1.000000   29.000000    0.000000           NaN   94.000000   
25%      76.500000   48.000000    0.000000           NaN  120.000000   
50%     152.000000   56.000000    1.000000           NaN  130.000000   
75%     227.500000   61.000000    1.000000           NaN  140.000000   
max     303.000000   77.000000    1.000000           NaN  200.000000   

              Chol         Fbs     RestECG       MaxHR       ExAng  \
count   303.000000  303.000000  303.000000  303.000000  303.000000   
unique         NaN         NaN         NaN         NaN         NaN   
top            NaN         NaN         NaN         NaN         NaN   
freq           NaN         NaN         NaN         NaN         NaN   
mean    246.693069    0.148515    0.990099  149.607261    0.326733   
std      51.776918    0.356198    0.994971   22.875003    0.469794   
min     126.000000    0.000000    0.000000   71.000000    0.000000   
25%     211.000000    0.000000    0.000000  133.500000    0.000000   
50%     241.000000    0.000000    1.000000  153.000000    0.000000   
75%     275.000000    0.000000    2.000000  166.000000    1.000000   
max     564.000000    1.000000    2.000000  202.000000    1.000000   

           Oldpeak       Slope          Ca    Thal  AHD  
count   303.000000  303.000000  299.000000     301  303  
unique         NaN         NaN         NaN       3    2  
top            NaN         NaN         NaN  normal   No  
freq           NaN         NaN         NaN     166  164  
mean      1.039604    1.600660    0.672241     NaN  NaN  
std       1.161075    0.616226    0.937438     NaN  NaN  
min       0.000000    1.000000    0.000000     NaN  NaN  
25%       0.000000    1.000000    0.000000     NaN  NaN  
50%       0.800000    2.000000    0.000000     NaN  NaN  
75%       1.600000    2.000000    1.000000     NaN  NaN  
max       6.200000    3.000000    3.000000     NaN  NaN  

Graphical summary:

# Graphical summaries
plt.figure()
sns.pairplot(data = Heart);
plt.show()

3 Initial split into test and non-test sets

We randomly split the data in half of test data and another half of non-test data. Stratify on AHD.

Python

from sklearn.model_selection import train_test_split

Heart_other, Heart_test = train_test_split(
  Heart, 
  train_size = 0.75,
  random_state = 425, # seed
  stratify = Heart.AHD
  )
Heart_test.shape

(76, 15)

Heart_other.shape

(227, 15)

Separate \(X\) and \(y\). We will use 13 features.

num_features = ['Age', 'Sex', 'RestBP', 'Chol', 'Fbs', 'RestECG', 'MaxHR', 'ExAng', 'Oldpeak', 'Slope', 'Ca']
cat_features = ['ChestPain', 'Thal']
features = np.concatenate([num_features, cat_features])
# Non-test X and y
X_other = Heart_other[features]
y_other = Heart_other.AHD
# Test X and y
X_test = Heart_test[features]
y_test = Heart_test.AHD

4 Preprocessing (Python) or recipe (R)

Python

There are missing values in Ca (quantitative) and Thal (qualitative) variables. We are going to use simple mean imputation for Ca and most_frequent imputation for Thal. This is suboptimal. Better strategy is to use multiple imputation.

# How many NaNs
Heart.isna().sum()

Unnamed: 0    0
Age           0
Sex           0
ChestPain     0
RestBP        0
Chol          0
Fbs           0
RestECG       0
MaxHR         0
ExAng         0
Oldpeak       0
Slope         0
Ca            4
Thal          2
AHD           0
dtype: int64

In principle, decision trees should be able to handle categorical predictors. However scikit-learn and xgboost implementations don’t allow categorical predictors and require one-hot encoding.

from sklearn.preprocessing import OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline

# Transformer for categorical variables
categorical_tf = Pipeline(steps = [
  ("cat_impute", SimpleImputer(strategy = 'most_frequent')),
  ("encoder", OneHotEncoder())
])

# Transformer for continuous variables
numeric_tf = Pipeline(steps = [
  ("num_impute", SimpleImputer(strategy = 'mean')),
])

# Column transformer
col_tf = ColumnTransformer(transformers = [
  ('num', numeric_tf, num_features),
  ('cat', categorical_tf, cat_features)
])

5 Model

Python

from sklearn.tree import DecisionTreeClassifier, plot_tree

classtree_mod = DecisionTreeClassifier(
  criterion = 'gini',
  random_state = 425
  )

6 Pipeline (Python) or workflow (R)

Here we bundle the preprocessing step (Python) or recipe (R) and model.

Python

from sklearn.pipeline import Pipeline

pipe = Pipeline(steps = [
  ("col_tf", col_tf),
  ("model", classtree_mod)
  ])
pipe

Pipeline(steps=[('col_tf',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('num_impute',
                                                                   SimpleImputer())]),
                                                  ['Age', 'Sex', 'RestBP',
                                                   'Chol', 'Fbs', 'RestECG',
                                                   'MaxHR', 'ExAng', 'Oldpeak',
                                                   'Slope', 'Ca']),
                                                 ('cat',
                                                  Pipeline(steps=[('cat_impute',
                                                                   SimpleImputer(strategy='most_frequent')),
                                                                  ('encoder',
                                                                   OneHotEncoder())]),
                                                  ['ChestPain', 'Thal'])])),
                ('model', DecisionTreeClassifier(random_state=425))])

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

7 Tuning grid

ccp_alpha is the Minimal Cost-Complexity Pruning parameter. Greater values of ccp_alpha increase the number of nodes pruned.

Python

# Tune hyper-parameter(s)
ccp_alpha_grid = np.linspace(start = 0.0, stop = 0.05, num = 100)
tuned_parameters = {
  "model__ccp_alpha": ccp_alpha_grid
  }
tuned_parameters  

{'model__ccp_alpha': array([0.        , 0.00050505, 0.0010101 , 0.00151515, 0.0020202 ,
       0.00252525, 0.0030303 , 0.00353535, 0.0040404 , 0.00454545,
       0.00505051, 0.00555556, 0.00606061, 0.00656566, 0.00707071,
       0.00757576, 0.00808081, 0.00858586, 0.00909091, 0.00959596,
       0.01010101, 0.01060606, 0.01111111, 0.01161616, 0.01212121,
       0.01262626, 0.01313131, 0.01363636, 0.01414141, 0.01464646,
       0.01515152, 0.01565657, 0.01616162, 0.01666667, 0.01717172,
       0.01767677, 0.01818182, 0.01868687, 0.01919192, 0.01969697,
       0.02020202, 0.02070707, 0.02121212, 0.02171717, 0.02222222,
       0.02272727, 0.02323232, 0.02373737, 0.02424242, 0.02474747,
       0.02525253, 0.02575758, 0.02626263, 0.02676768, 0.02727273,
       0.02777778, 0.02828283, 0.02878788, 0.02929293, 0.02979798,
       0.03030303, 0.03080808, 0.03131313, 0.03181818, 0.03232323,
       0.03282828, 0.03333333, 0.03383838, 0.03434343, 0.03484848,
       0.03535354, 0.03585859, 0.03636364, 0.03686869, 0.03737374,
       0.03787879, 0.03838384, 0.03888889, 0.03939394, 0.03989899,
       0.04040404, 0.04090909, 0.04141414, 0.04191919, 0.04242424,
       0.04292929, 0.04343434, 0.04393939, 0.04444444, 0.04494949,
       0.04545455, 0.0459596 , 0.04646465, 0.0469697 , 0.04747475,
       0.0479798 , 0.04848485, 0.0489899 , 0.04949495, 0.05      ])}

8 Cross-validation (CV)

Python

Set up CV partitions and CV criterion.

from sklearn.model_selection import GridSearchCV

# Set up CV
n_folds = 5
search = GridSearchCV(
  pipe,
  tuned_parameters,
  cv = n_folds, 
  scoring = "roc_auc",
  # Refit the best model on the whole data set
  refit = True
  )

Fit CV. This is typically the most time-consuming step.

# Fit CV
search.fit(X_other, y_other)

GridSearchCV(cv=5,
             estimator=Pipeline(steps=[('col_tf',
                                        ColumnTransformer(transformers=[('num',
                                                                         Pipeline(steps=[('num_impute',
                                                                                          SimpleImputer())]),
                                                                         ['Age',
                                                                          'Sex',
                                                                          'RestBP',
                                                                          'Chol',
                                                                          'Fbs',
                                                                          'RestECG',
                                                                          'MaxHR',
                                                                          'ExAng',
                                                                          'Oldpeak',
                                                                          'Slope',
                                                                          'Ca']),
                                                                        ('cat',
                                                                         Pipeline(steps=[('cat_impute',
                                                                                          SimpleImputer(strategy='most_frequent')),
                                                                                         ('encoder',
                                                                                          OneHotEncoder())]),
                                                                         ['ChestPain',
                                                                          '...
       0.03282828, 0.03333333, 0.03383838, 0.03434343, 0.03484848,
       0.03535354, 0.03585859, 0.03636364, 0.03686869, 0.03737374,
       0.03787879, 0.03838384, 0.03888889, 0.03939394, 0.03989899,
       0.04040404, 0.04090909, 0.04141414, 0.04191919, 0.04242424,
       0.04292929, 0.04343434, 0.04393939, 0.04444444, 0.04494949,
       0.04545455, 0.0459596 , 0.04646465, 0.0469697 , 0.04747475,
       0.0479798 , 0.04848485, 0.0489899 , 0.04949495, 0.05      ])},
             scoring='roc_auc')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Visualize CV results.

Code

cv_res = pd.DataFrame({
  "ccp_alpha": np.array(search.cv_results_["param_model__ccp_alpha"]),
  "auc": search.cv_results_["mean_test_score"]
  })

plt.figure()
sns.relplot(
  # kind = "line",
  data = cv_res,
  x = "ccp_alpha",
  y = "auc"
  ).set(
    xlabel = "CCP Alpha",
    ylabel = "CV AUC"
);
plt.show()

Best CV AUC:

search.best_score_

0.803652380952381

The training accuracy is

from sklearn.metrics import accuracy_score, roc_auc_score

accuracy_score(
  y_other,
  search.best_estimator_.predict(X_other)
  )

0.8149779735682819

9 Finalize our model

Now we are done tuning. Finally, let’s fit this final model to the whole training data and use our test data to estimate the model performance we expect to see with new data.

Python

Since we called GridSearchCV with refit = True, the best model fit on the whole non-test data is readily available.

search.best_estimator_

Pipeline(steps=[('col_tf',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('num_impute',
                                                                   SimpleImputer())]),
                                                  ['Age', 'Sex', 'RestBP',
                                                   'Chol', 'Fbs', 'RestECG',
                                                   'MaxHR', 'ExAng', 'Oldpeak',
                                                   'Slope', 'Ca']),
                                                 ('cat',
                                                  Pipeline(steps=[('cat_impute',
                                                                   SimpleImputer(strategy='most_frequent')),
                                                                  ('encoder',
                                                                   OneHotEncoder())]),
                                                  ['ChestPain', 'Thal'])])),
                ('model',
                 DecisionTreeClassifier(ccp_alpha=0.0202020202020202,
                                        random_state=425))])

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Visualize the best classification tree.

features = np.concatenate([
    features[:-2], 
    ['ChestPain:asymptomatic', 'ChestPain:nonanginal', 'ChestPain:nontypical', 'ChestPain:typical'],
    ['Thal:fixed', 'Thal:normal', 'Thal:reversable']
    ])

plt.figure()
plot_tree(
  search.best_estimator_['model'],
  feature_names = features
  );
plt.show()

Feature importances:

vi_df = pd.DataFrame({
  "feature": features,
  "vi": search.best_estimator_['model'].feature_importances_
  })

plt.figure()
sns.barplot(
  data = vi_df,
  x = "feature",
  y = "vi"
  ).set(
    xlabel = "Feature",
    ylabel = "VI"
);
plt.xticks(rotation = 90);
plt.show()

The final AUC on the test set is

roc_auc_score(
  y_test,
  search.best_estimator_.predict_proba(X_test)[:, 1]
  )

0.8574912891986062

The final classification accuracy on the test set is

accuracy_score(
  y_test, 
  search.best_estimator_.predict(X_test)
  )

0.8421052631578947