Helping detect a rare disease
Below I construct a decision-tree model and a random-forest model to help detect a condition with an incidence of less than 0.173% among patients (173 in 100,000). The models are built using the anonymized information of more than 280,000 patients, consisting of 30 anonymized variables per patient. Each patient is classified as either having the condition 'Class'=1, or not having it 'Class'=0.
Author: Leonardo Espin
Date: August 8, 2019
- First I will load all the libraries I will use for the exploratory data analysis
In [1]:
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
pd.set_option('display.max_columns', 35)
#the lines below get rid of the plotting library warnings
import warnings
warnings.filterwarnings("ignore")
- Next I will read the data into a dataframe, print the size of the dataset, the first few lines of data and the types of variables
In [2]:
raw=pd.read_csv('data_challenge_dataset.csv')
print(raw.shape)
print(raw.dtypes)
raw.head()
Out[2]:
Exploratory data analysis (EDA)¶
- First I print summary statistics of the dataset
In [3]:
raw.describe()
Out[3]:
- The data seems to have been centered, but it has not been scaled to unit variance
- Below I check for missing data, first along features, then across columns
In [4]:
print(raw.isna().sum())
print('Percentage of missing data in each column: {:.1f}%'.format(100*raw.V1.isna().sum()/len(raw)))
Therefore, all the columns have exactly the same amount of missing data. Then, if we look across columns:
In [5]:
naInRows=raw.isna().sum(axis=1)
print('Percentage of rows with missing data: {:.1f}%'.format(100*len(raw[raw.isna().sum(axis=1)>0])/len(raw)))
plt.figure(figsize=(16,4))
naInRows.plot()
plt.title('Missing data for each row', fontsize=14)
plt.xlabel('Row number', fontsize=14)
plt.ylabel('Number of missing', fontsize=14);
- All the columns have the same number of missing values (10% of total)
- Most rows have missing data, and the plot above shows that the missing data is more or less uniformly distributed along rows
- Importantly, there are not missing values in the target column
- I will use simple inputation before modeling the data, which I will explain in more detail in the modeling portion of the document
- Below I inspect the distribution of the target variable:
In [6]:
raw.Class.value_counts()
Out[6]:
- So, this is a heavily imbalanced dataset
- Now, I will inspect the distribution of individual variables, but before doing that I will partition the dataset into train/test portions and will keep the test portion untouched until I use it for testing the model I will create.
- Notice that I use stratified partitioning, in order to maintain the proportions of the classification variable in the test and train portions of the dataset
In [7]:
#Here I'm creating a list of features, and I separate the target variable
features=list(raw.columns)
target=features.pop()
In [8]:
from sklearn.model_selection import train_test_split
seed=491
train,test = train_test_split(raw,test_size=0.3,random_state=seed,stratify=raw[target])
In [9]:
#checking that the proportions have been maintained
print(100*train[target].value_counts()/len(train))
print(100*test[target].value_counts()/len(test))
- Below I plot the histograms of all the individual features to have a sense of the distribution of each variable
- Note that from this point on, I will only use the train portion of the dataset for analysis and modeling. The test portion will only be used for testing a model
In [10]:
plt.figure(figsize=(16,25))
for i in range(len(features)):
plt.subplot(8, 4, i+1)
ax=sns.distplot(train[features[i]].dropna(),kde=False,rug=False)#NaN causes problems
ax.set_yscale('log')
- All the variables have outliers (a few single values in the right or left tail of the count plots). Unfortunately, since the features are anonymous, it is difficult to guess the meaning of outlier values for each particular variable
- Below I plot the distribution of features, segregated by the value of the target variable:
In [11]:
def segregatedPlot(df,col,grouping='Class'):
'''conveniently making segregated plots'''
#plt.figure(figsize=(9,6))
sns.distplot(df.groupby(grouping)[col].get_group(0).dropna(),kde=False,rug=False, label=grouping.capitalize()+'=0');
ax=sns.distplot(df.groupby(grouping)[col].get_group(1).dropna(),kde=False,rug=False, label=grouping.capitalize()+'=1');
ax.set_yscale('log')
plt.legend();
In [12]:
plt.figure(figsize=(16,30))
for i in range(len(features)):
plt.subplot(10, 3, i+1)
segregatedPlot(train,features[i])
The plots above are very illuminating with respect to what variables have might have more predictive power. Some of the variables show different distributions depending on the value of the class that a given row belongs to. I can see six different patterns in the plots above:
- For some variables the target classification makes little difference in the distribution of the variable
- The variable distribution shows a lower peak count for class 1 (V3)
- The variable distribution shows a larger peak count for class 1 (V4)
- Lower peak count and lower tail values (V9, V10, V12, V14, V16, V17)
- Larger peak count and larger tail values (V11)
- Lower tail values (V18)
Because of the apparent importance of tail values on the class distribution of the target variable, I am not inclined to use classification methods that are sensitive to feature scaling (like logistic regression). So I will first try using decision trees and random forest, and see how well they perform on the the test data.
- Before start constructing a model, I inspect variable correlations below:
In [13]:
corr = train[features].corr()
plt.figure(figsize=(10,10))
#diverging colormap
cmap = sns.diverging_palette(220,10,as_cmap=True)
# Generate a mask for the upper triangle
mask = np.zeros_like(corr, dtype=np.bool)
#triu_indices_from() returns the indices for the upper-triangle of arr
mask[np.triu_indices_from(mask)] = True
sns.heatmap(corr, mask=mask, cmap=cmap, center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .5});
- There is little correlation between the features, with the exception of V29 (which from the plots above, seems to have very little predictive power). This must be a consequence of the variables being centered, as pointed out earlier in the document
- Note that V29 is the only variable that is in a scale that is significantly different than the other variables (only positive values, on the order of a thousand, it could be length of stay at a hospital, for example).
Modeling the data¶
Given the seemingly uniform distribution of the missing data, and the variable distribution plots shown above, I use simple imputation to fill the missing data
Baseline decision-tree model¶
In [14]:
from sklearn import tree
bench = tree.DecisionTreeClassifier(random_state=seed)
#simple imputation:
train_f=train.fillna(train.mean())
bench.fit(train_f[features],train_f[target])
Out[14]:
In [15]:
print('Maximum depth = {}'.format(bench.tree_.max_depth))
print('Default score = {:.3f}'.format(bench.score(train_f[features],train_f[target])))
print('Default feature importances (sorted by importance):')
tups=[(features[_],bench.feature_importances_[_]) for _ in range(len(features))]
tups.sort(key=lambda x:x[1], reverse=True)
for tup in tups:
print('{}, importance = {:.4f}'.format(tup[0],tup[1]))
- Unsurprisingly the benchmark decision tree overfits the data. An interesting observation is that the top three features for the tree are in my list of variables which I expected to have more predictive power.
- Below I use cross validation with stratified partitioning of the train dataset to find an optimal depth for a tree
In [16]:
from sklearn.model_selection import StratifiedKFold, GridSearchCV
benchCV = StratifiedKFold(n_splits=6,random_state=seed)
param_test = {
'max_depth':[2,5,20,30]}#[_ for _ in range(2,30)] produces 5
dptSearch = GridSearchCV(
estimator=tree.DecisionTreeClassifier(random_state=seed),
param_grid = param_test, scoring='roc_auc',n_jobs=4,
iid=False, cv=benchCV,return_train_score=True)
dptSearch.fit(train_f[features],train_f[target])
dptSearch.best_params_, dptSearch.best_score_
Out[16]:
- The best cross-validation AUC score obtained with this tree is 88%, which is good. Below I check the ROC-AUC score on the test set, after fitting a decision tree with a maximum depth of 5 nodes:
In [17]:
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix
bench = tree.DecisionTreeClassifier(random_state=seed,max_depth=5)
bench.fit(train_f[features],train_f[target])
#simple imputation of test dataset:
test_f=test.fillna(test.mean())
print('Decision-tree ROC-AUC score = {:.3f}'.format(roc_auc_score(test_f[target],bench.predict_proba(test_f[features])[:,1])))
print('Feature importances (sorted by importance):')
tups=[(features[_],bench.feature_importances_[_]) for _ in range(len(features))]
tups.sort(key=lambda x:x[1], reverse=True)
for tup in tups:
print('{}, importance = {:.4f}'.format(tup[0],tup[1]))
- How it should be expected, the ROC-AUC score on the test data is lower than on cross-validation, but not by much (84% vs. 88%)
- Note that the feature importances agree pretty well, for the most part, with the variables I highlighted in the EDA portion of this analysis
- Below I calculate the confusion matrix for the predictions of the decision-tree model
In [18]:
#the confusion matrix for the prediction is below
print('\tTn\tFp\n\tFn\tTp')
confusion_matrix(test_f[target],bench.predict(test_f[features]))
Out[18]:
- The confusion matrix above shows reasonable results as well (sensitivity: 73%, specificity: 99.9%)
Constructing a random-forest model¶
- Firs I use
SelectFromModel, below, to select features for the model based on importance weights. The threshold for selection calculated automatically bySelectFromModelis0.034483, and the features selected automatically are V10, V11, V12, V14, V16 and V17. - Note that I listed all these variables as likely to have more predictive power in the EDA portion of the analysis.
- Since there are a few variables that I listed that were not selected automatically (e.g. V3, V4, V9 and V18), I performed a brief sensitivity analysis and slightly reduced the automatic threshold to see if some of these variables would be picked up:
In [19]:
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectFromModel #for selecting features for the model
sel = SelectFromModel(RandomForestClassifier(n_estimators = 100,random_state=seed,n_jobs=4),
threshold=0.03)
sel.fit(train_f[features], train_f[target])
Out[19]:
In [20]:
selected=[features[i] for i in range(len(features)) if sel.get_support()[i]]
print('The selected variables for my random forest model are: ')
for _ in selected:
print(_)
- As noted above, I listed all these variables as likely to have more predictive power in the EDA portion of the analysis
- Below I use cross validation to find the optimal number of trees for the random forest model.
In [21]:
param_test2 = {
'n_estimators':[2,50,76,100]}#[_ for _ in range(2,100)] selects 76
forestSearch = GridSearchCV(
estimator=RandomForestClassifier(random_state=seed),
param_grid = param_test2, scoring='roc_auc',n_jobs=4,
iid=False, cv=benchCV,return_train_score=True)
forestSearch.fit(train_f[selected], train_f[target])
forestSearch.best_params_, forestSearch.best_score_
Out[21]:
- The best cross-validation AUC score obtained with this random-forest model is 92.9%.
- Below I calculate the confusion matrix for the predictions of the random-forest model, on the test set:
In [22]:
#confusion matrix for the best forest model
print('\tTn\tFp\n\tFn\tTp')
confusion_matrix(test_f[target],forestSearch.predict(test_f[selected]))
Out[22]:
- The random-forest model sensitivity and specificity are 79.7% and 99.9%, respectively
Comparing the models by constructing ROC curves¶
In [23]:
from sklearn.metrics import auc #will use below, when plotting the curves
from roc_curve import get_roc_curve #I put this code in a separate file for clarity
#making an array of thresholds
thresholds=np.linspace(0,1,15)
- Below I construct and plot an ROC curve for my decision-tree model
In [24]:
plt.figure(figsize=(10,7))
tprate=[]
fprate=[]
aucs=[]
counter=0
for tr_fold, test_fold in benchCV.split(train_f[selected],train_f[target]):
bench.fit(train_f[selected].iloc[tr_fold],train_f[target].iloc[tr_fold])
probs=bench.predict_proba(train_f[selected].iloc[test_fold])[:,1]#prob of class 1
tpr,fpr=get_roc_curve(train_f[target].iloc[test_fold].values,probs,thresholds)
tprate.append(tpr)
fprate.append(fpr)
auc_fold=auc(fpr, tpr)
aucs.append(auc_fold)
plt.plot(fpr,tpr,lw=1,alpha=0.3,
label='Fold {}, AUC={:.3f}'.format(counter,auc_fold))
#plt.plot(fpr[index],tpr[index],'r+');
counter+=1
#chance curve
plt.plot([0, 1], [0, 1],linestyle='--',lw=1, alpha=.8,
color='r',label='Chance')
mean_tpr = np.mean(tprate, axis=0)
mean_fpr = np.mean(fprate, axis=0)
mean_auc = np.mean(aucs)
std_auc = np.std(aucs)
plt.plot(mean_fpr, mean_tpr,lw=2, color='b',
label='Mean ROC (AUC={:.3f} $\pm$ {:.3f})'.format(mean_auc, std_auc))
#the cross shows the value corresponding to an standard threshold of 0.5
index=np.where(thresholds ==0.5)[0][0]#assuming that 0.5 is present in thresholds
plt.plot(mean_fpr[index],mean_tpr[np.where(thresholds ==0.5)[0][0]],'r+')
plt.xlabel('False Positive Rate', fontsize=14)
plt.ylabel('True Positive Rate', fontsize=14)
plt.title('ROC decision-tree model', fontsize=18)
plt.legend(loc="lower right", prop={'size': 12});
- The red cross in this plot shows the true and false positive rates for the standard classification threshold of 0.5. Similar crosses are also shown in the plots below
In [25]:
forest = RandomForestClassifier(random_state=seed,n_estimators=76)
In [26]:
plt.figure(figsize=(10,7))
tpr_forest=[]
fpr_forest=[]
aucs_forest=[]
counter=0
for tr_fold, test_fold in benchCV.split(train_f[selected],train_f[target]):
forest.fit(train_f[selected].iloc[tr_fold],train_f[target].iloc[tr_fold])
probs=forest.predict_proba(train_f[selected].iloc[test_fold])[:,1]#prob of class 1
tpr,fpr=get_roc_curve(train_f[target].iloc[test_fold].values,probs,thresholds)
tpr_forest.append(tpr)
fpr_forest.append(fpr)
auc_fold=auc(fpr, tpr)
aucs_forest.append(auc_fold)
plt.plot(fpr,tpr,lw=1,alpha=0.3,
label='Fold {}, AUC={:.3f}'.format(counter,auc_fold))
#plt.plot(fpr[index],tpr[index],'r+');
counter+=1
#chance curve
plt.plot([0, 1], [0, 1],linestyle='--',lw=1, alpha=.8,
color='r',label='Chance')
mean_tpr_f = np.mean(tpr_forest, axis=0)
mean_fpr_f = np.mean(fpr_forest, axis=0)
mean_auc_f = np.mean(aucs_forest)
std_auc_f = np.std(aucs_forest)
plt.plot(mean_fpr_f, mean_tpr_f,lw=2, color='b',
label='Mean ROC (AUC={:.3f} $\pm$ {:.3f})'.format(mean_auc_f, std_auc_f))
#the cross shows the value corresponding to an standard threshold of 0.5
index=np.where(thresholds ==0.5)[0][0]#assuming that 0.5 is present in thresholds
plt.plot(mean_fpr_f[index],mean_tpr_f[index],'r+')
plt.xlabel('False Positive Rate', fontsize=14)
plt.ylabel('True Positive Rate', fontsize=14)
plt.title('ROC random-forest model', fontsize=18)
plt.legend(loc="lower right", prop={'size': 12});
- Below I compare the mean ROC curves for both models:
In [27]:
#comparing the two models
plt.figure(figsize=(10,7))
#decision tree ROC
plt.plot(mean_fpr, mean_tpr,lw=2,
label='Dtree Mean ROC (AUC={:.2f} $\pm$ {:.2f})'.format(mean_auc, std_auc))
plt.plot(mean_fpr_f, mean_tpr_f,lw=2,
label='Rforest Mean ROC (AUC={:.2f} $\pm$ {:.2f})'.format(mean_auc_f, std_auc_f))
#chance curve
plt.plot([0, 1], [0, 1],linestyle='--',lw=1, alpha=.8,
color='r',label='Chance')
plt.plot(mean_fpr[index],mean_tpr[index],'r+')
plt.plot(mean_fpr_f[index],mean_tpr_f[index],'g+')
plt.xlabel('False Positive Rate', fontsize=14)
plt.ylabel('True Positive Rate', fontsize=14)
plt.title('ROC for the two models', fontsize=18)
plt.legend(loc="lower right", prop={'size': 12});
- Since the crosses show the standard classification threshold of 0.5, it is clear that we could reduce this threshold in order to improve the true positive rate, without significantly reducing the false positive rate.
- To investigate this in more detail, I show a table with the mean true and false positive rates for both the decision tree and random forest models, together with the corresponding thresholds below:
In [28]:
rates=pd.concat([pd.Series(thresholds,name='Thresholds'),
pd.Series(mean_tpr,name='tree_tpr'),
pd.Series(mean_fpr,name='tree_fpr'),
pd.Series(mean_tpr_f,name='forest_tpr'),
pd.Series(mean_fpr_f,name='forest_fpr')
],axis=1)
rates
Out[28]:
- As can be seen in the table above, we can indeed reduce the classification threshold without affecting the false positive rate considerably. This would be equivalent to increasing the number of true positives or reducing the false negatives.
- Since the number of actual positives is small, we would be increasing the number of false positives as well, which might not be desirable.
Conclusions¶
- The exploratory data analysis reveals several interesting trends that can be observed in the data. In particular, the distributions corresponding to each class of the target variable (0 or 1) are different for a subset of features (variables V3, V4, V9 to V12, V14, V16, V17 and V18).
- For a few varibles, outliers in the distribution seem to correlate with belonging to class 1
- Missing data seems to be distributed randomly and uniformly across the data, totaling around 10% for each column.
- Feature importances for my decision tree model, and feature selection with a random forest model correspond well with the variables highlighted in points 1 and 2
- The performance of a random forest model with 76 trees improves performance with respect to my decision-tree model, decreasing both the number of false positives and false negatives (sensitivity: 79.7%, specificity: 99.9%)
- As shown in the ROC-curve analysis section, the threshold for classification can be changed to optimize the number of false positives or false negatives. Selecting this threshold would probably require external input