Deep Exploratory Analysis and Random Forest Classification

Author(s): Greg Postalian-Yrausquin

Originally published on Towards AI.

Decision tree types of classification algorithms have the advantage that they produce results that are relatively easier to explain in terms of the impact of the predictors when compared to other supervised training algorithms, like neural networks, which is more like a closed box.

To illustrate this, I use RandomForest and descriptive analytics see if it is possible to predict if a person earns more or less than 50k dollars

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import klib
import dabl
from scipy.stats import chi2_contingency
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.metrics import confusion_matrix
from sklearn import tree

From the description of the dataset, this is what we know about the predictors:

age: continuous.

workclass: Private, Self-emp-not-inc, Self-emp-inc, Federal-gov, Local-gov, State-gov, Without-pay, Never-worked.

fnlwgt: continuous.

education: Bachelors, Some-college, 11th, HS-grad, Prof-school, Assoc-acdm, Assoc-voc, 9th, 7th-8th, 12th, Masters, 1st-4th, 10th, Doctorate, 5th-6th, Preschool.

education-num: continuous.

marital-status: Married-civ-spouse, Divorced, Never-married, Separated, Widowed, Married-spouse-absent, Married-AF-spouse.

occupation: Tech-support, Craft-repair, Other-service, Sales, Exec-managerial, Prof-specialty, Handlers-cleaners, Machine-op-inspct, Adm-clerical, Farming-fishing, Transport-moving, Priv-house-serv, Protective-serv, Armed-Forces.

relationship: Wife, Own-child, Husband, Not-in-family, Other-relative, Unmarried.

race: White, Asian-Pac-Islander, Amer-Indian-Eskimo, Other, Black.

sex: Female, Male.

capital-gain: continuous.

capital-loss: continuous.

hours-per-week: continuous.

native-country: United-States, Cambodia, England, Puerto-Rico, Canada, Germany, Outlying-US(Guam-USVI-etc), India, Japan, Greece, South, China, Cuba, Iran, Honduras, Philippines, Italy, Poland, Jamaica, Vietnam, Mexico, Portugal, Ireland, France, Dominican-Republic, Laos, Ecuador, Taiwan, Haiti, Columbia, Hungary, Guatemala, Nicaragua, Scotland, Thailand, Yugoslavia, El-Salvador, Trinadad&Tobago, Peru, Hong, Holand-Netherlands.

maindataset = pd.read_csv("adultsalary.csv")
print(maindataset)
print(maindataset.columns)

Right from the start, I can see that some columns are going to be correlated (marital status — relationship, education and education numerical) and some others are basically empty.

The next step is a review of the numerical variables

klib.dist_plot(maindataset['age'])
klib.dist_plot(maindataset[' education-num'])
klib.dist_plot(maindataset[' hours-per-week'])

For the categorical variables:

#I am going to group the countries of origin to reduce the number of values of the variable
us = [' United-States']
latin = [' Cuba',' South',' Puerto-Rico',' Honduras',' Columbia',' Ecuador',' Dominican-Republic',' El-Salvador',' Guatemala',' Peru',' Nicaragua']
carib = [' Jamaica',' Haiti',' Outlying-US(Guam-USVI-etc)',' Trinadad&Tobago']
mexico = [' Mexico']
canada = [' Canada']
westeur = [' England',' Germany',' France',' Scotland',' Ireland',' Holand-Netherlands']
easteur = [' Poland',' Yugoslavia',' Hungary']
southeur = [' Italy',' Portugal',' Greece']
iran = [' Iran']
india = [' India']
seasia = [' Philippines',' Cambodia',' Thailand',' Laos',' Vietnam']
china = [' China']
easia = [' Taiwan',' Japan',' Hong']

maindataset['region'] = np.nan
maindataset['region'] = np.where(maindataset[' native-country'].isin(us),'US',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(latin),'Other LA',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(carib),'Caribbean',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(mexico),'Mexico',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(canada),'Canada',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(westeur),'West Euro',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(easteur),'East Euro',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(southeur),'South Euro',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(iran),'Iran',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(india),'India',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(seasia),'SE Asia',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(china),'China',maindataset['region'])
maindataset['region'] = np.where(maindataset[' native-country'].isin(easia),'E Asia',maindataset['region'])


maindataset[' salary'] = np.where(maindataset[' salary']==' <=50K.',' <=50K',maindataset[' salary'])
maindataset[' salary'] = np.where(maindataset[' salary']==' >50K.',' >50K',maindataset[' salary'])

klib.cat_plot(maindataset[[' workclass', ' occupation', ' relationship', ' race', ' sex', 'region', ' salary']], top=5, bottom=5)

Now I apply the filters and clean the data

maindatasetF = maindataset[['age',' education-num',' hours-per-week',' workclass', ' occupation', ' relationship', ' race', ' sex', 'region', ' salary']]
maindatasetF = maindatasetF[maindatasetF[' workclass']!=' ?']
maindatasetF = maindatasetF[maindatasetF[' occupation']!=' ?']
maindatasetF = maindatasetF.dropna()
maindatasetF

Exploring the correlations.

For numerical data:

Which are not strong, positive, correlations

For the categorical data two tests are important:

1. Chi-square: to test for independence between the categorical variables.
2. ANOVA: to determine if between the values of the categorical variables there is a difference in the outcome (salary range).

For the Chi-square test:

maindatasetF = maindatasetF.rename(columns={' salary':'salary', ' workclass':'workclass', ' occupation':'occupation', ' relationship':'relationship', ' race':'race', ' sex':'sex'})

catvar = ['workclass', 'occupation', 'relationship', 'race', 'sex', 'region', 'salary']

xtabs = []
pairs = []
for i in catvar:
 for j in catvar:
 if i!=j:
 pair = pd.DataFrame([[i,j]], columns=['i','j'])
 data_xtab = pd.crosstab(maindatasetF[i],maindatasetF[j],margins=False)
 xtabs.append(data_xtab)
 pairs.append(pair)
pairs = pd.concat(pairs, ignore_index=True, axis=0)

ps = []
for i in xtabs:
 stat, p, dof, expected = chi2_contingency(i)
 ps.append(p)

pairs['p values'] = ps
pairs

These very low values of p mean that there is dependency between the features.

ANOVA test

import statsmodels.api as sm
from statsmodels.formula.api import ols

model = ols('num_salary ~ region + workclass + occupation + relationship + race + sex', data=maindatasetF).fit()
aov_table = sm.stats.anova_lm(model, typ=2)
aov_table

These low values of P mean that we can reject the NULL hypothesis, specifying that there is a difference between the members of the categorical features, and so, there is predictive power in them.

Next, it’s a random forest classification, and get the variable importance.

num_columns = ['age',' education-num',' hours-per-week']
cat_columns = ['workclass', 'occupation', 'relationship', 'race', 'sex', 'region']
cat_preprocessor = OneHotEncoder(handle_unknown="ignore")
num_preprocessor = StandardScaler()
preprocessor = ColumnTransformer(
 [
 ("one-hot-encoder", cat_preprocessor, cat_columns),
 ("standard_scaler", num_preprocessor, num_columns),
 ])

train, test = train_test_split(maindatasetF, train_size=0.8)
train = train.dropna()

X_train = train[['age',' education-num',' hours-per-week','workclass', 'occupation', 'relationship', 'race', 'sex', 'region']]
Y_train = train[['num_salary']]

X = pd.DataFrame.sparse.from_spmatrix(preprocessor.fit_transform(X_train))
catnames = preprocessor.transformers_[0][1].get_feature_names_out(cat_columns).tolist()
numnames = preprocessor.transformers_[1][1].get_feature_names_out(num_columns).tolist()
featnames = catnames + numnames

rf = RandomForestClassifier(n_estimators=100)
rf.fit(X, Y_train)

Plot of the importance of the variables

imp = rf.feature_importances_
imp = pd.Series(imp, index=featnames)
std = pd.Series(np.std([tree.feature_importances_ for tree in rf.estimators_], axis=0), index=featnames)
fig, ax = plt.subplots()
imp.plot(kind='barh', yerr=std, ax=ax, figsize=(15,15))
ax.set_title("Feature importances using MDI")
ax.set_ylabel("Mean decrease in impurity")
fig.tight_layout()
plt.show()

Age, education, hours worked and the fact that the subject is a husband dominate in the outcome of the salary, which is in line with what we know of the society we live in.

Now, I am performing a test with a confusion matrix to see how accurately we can predict the salary from the Random Forest model.

X_test = test[['age',' education-num',' hours-per-week','workclass', 'occupation', 'relationship', 'race', 'sex', 'region']]
Y_test = test[['num_salary']]
Xt = pd.DataFrame.sparse.from_spmatrix(preprocessor.fit_transform(X_test))
Yt = rf.predict(Xt)
cm = confusion_matrix(Y_test['num_salary'], Yt)
sns.heatmap(cm, annot=True, cmap='PuBu', fmt='g')

Which can be improved, but it is a good start. At least we know that the model used for variable importance, which is what we want to understand, makes sense.

I would like to draw a decision tree as well, so see the variables act in a graph.

tr = tree.DecisionTreeClassifier(max_depth=4)
tr = tr.fit(X, Y_train)
#I am limiting the depth to three
fig, ax1 = plt.subplots(figsize=(25,15))
tree.plot_tree(tr, ax=ax1, feature_names=featnames, proportion=True, filled=True, fontsize=7)
plt.show()

Blue indicates higher tendency towards high salaries. Note that the numeric variables are standardized (0 is the mean). “Yes’s” are to the left, and “No’s” to the right.

So here is clear to see that 80% of: husbands, with higher education than the mean, that work more than the mean, and older than the mean, earn more than 50k. This group is only 9% of the total population.

On the other side 97% of wives that are under 1.3 standard deviation about the mean in education (so, they are more educated than the mean), work less than the mean earn less than 50%, and they represent 41% of the population.

Join thousands of data leaders on the AI newsletter. Join over 80,000 subscribers and keep up to date with the latest developments in AI. From research to projects and ideas. If you are building an AI startup, an AI-related product, or a service, we invite you to consider becoming a sponsor.

Published via Towards AI