In [1]:
!date
import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
%matplotlib inline
sns.set_context('poster')
sns.set_style('darkgrid')
In [2]:
# set random seed, for reproducibility
np.random.seed(12345)
Load and clean PHMRC VA data:
In [3]:
df = pd.read_csv('IHME_PHMRC_VA_DATA_ADULT_Y2013M09D11_0.csv')
What do you think we should do about that warning?
In [4]:
# just do what the warning recommends
df = pd.read_csv('IHME_PHMRC_VA_DATA_ADULT_Y2013M09D11_0.csv', low_memory=False)
In [5]:
# take a look at what is in the DataFrame
df.head()
Out[5]:
In [6]:
df.columns
Out[6]:
In [7]:
# also load codebook (excel doc)
cb = pd.read_excel('IHME_PHMRC_VA_DATA_CODEBOOK_Y2013M09D11_0.xlsx')
cb
Out[7]:
Every column that starts with a{NUMBER}_stuff is part of the signs and symptoms reported about an adult death:
In [8]:
df.filter(regex='a\d_')
Out[8]:
There is also a "bag of words" section from the processed results of the free-text response section:
In [9]:
df.filter(like='word_')
Out[9]:
Make the feature vectors, by putting together everything from the adult module, and everything from the free response section.
In [10]:
X = np.hstack((np.array(df.filter(regex='a\d_')), np.array(df.filter(like='word_'))))
Did I get that right?
In [11]:
np.array(df.filter(regex='a\d_')).shape, np.array(df.filter(like='word_')).shape
Out[11]:
In [12]:
X.shape
Out[12]:
Actually, that won't work without more data cleaning, because sklearn expects the feature vectors to be numbers, not answers like "Yes" and "Don't know".
So let's just use the open-response text, which has already been processed.
In [13]:
X = np.array(df.filter(like='word_'))
The causes of death appear in the gs_text34 column (obviously...)
In [14]:
df.gs_text34
Out[14]:
For class purposes, we will simplify the prediction task: was the death due to stroke, diabetes, or something else?
In [15]:
df['Cause'] = df.gs_text34.map({'Stroke':'Stroke', 'Diabetes':'Diabetes'}).fillna('Other')
y = np.array(df.Cause)
Did that work?
In [16]:
np.unique(y)
Out[16]:
Now, let's train one of our (by now) standard classifiers to predict the underlying cause of death from the restricted cause list:
In [17]:
import sklearn.tree
clf = sklearn.tree.DecisionTreeClassifier()
clf.fit(X,y)
y_pred = clf.predict(X)
np.mean(y == y_pred)
Out[17]:
Does that look good? Too good? If it looks too good to be true, it probably is...
In [18]:
import sklearn.cross_validation
scores = sklearn.cross_validation.cross_val_score(clf, X, y)
print scores.mean(), scores.std()
In [19]:
sns.distplot(scores, rug=True)
Out[19]:
What have we done here?
In [20]:
help(sklearn.cross_validation.cross_val_score)
If you read through that, you will find out: we have done stratified K-Fold cross validation with $K = 3$. In other words, split the data into 3 equal parts (what if it's not divisible by three?) and then use two for training and the third for testing, in all three ways.
So how to do 10-fold CV as in the reading?
In [21]:
%%time
n = len(y)
cv = sklearn.cross_validation.KFold(n, n_folds=10, shuffle=True)
scores = sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)
print scores.mean(), scores.std()
In [22]:
sns.distplot(scores, rug=True)
Out[22]:
I think it is important to stratify the sampling, although our book does not make a big deal about that:
In [24]:
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True)
scores = sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)
print scores.mean(), scores.std()
In [25]:
sns.distplot(scores, rug=True)
Out[25]:
And I don't want you to ever rely on just 10 samples. How can you do ten repetitions of 10-fold cross validation?
In [28]:
scores = []
for rep in range(10):
print rep,
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True)
scores += list(sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv))
In [29]:
np.mean(scores), np.std(scores)
Out[29]:
In [30]:
sns.distplot(scores, rug=True)
Out[30]:
Although our books make a big deal about K-Fold cross-validation, and we have gone through all of the trouble to understand what it is doing, I am partial to what is sometimes called leave-group-out cross-validation, Monte Carlo cross-validation, or just plain-old cross-validation: split the data into training and test sets randomly a bunch of times.
Just to be perverse, in sklearn this is called ShuffleSplit cross-validation:
In [32]:
%%time
n = len(y)
cv = sklearn.cross_validation.ShuffleSplit(n, n_iter=10, test_size=.25)
scores = sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)
print scores.mean(), scores.std()
Of course, you can also do a stratified version of this:
In [34]:
cv = sklearn.cross_validation.StratifiedShuffleSplit(y, n_iter=10, test_size=.25)
scores = sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)
print scores.mean(), scores.std()
One last way to do cross-validation, which (in my humble opinion) was not emphasized sufficiently in the reading: split based on what you are actually interested in. In the case of VA, I am very interested in how the approach will generalize to other parts of the world. So I could hold-out based on site:
In [35]:
df.site.value_counts()
Out[35]:
In [36]:
cv = sklearn.cross_validation.LeaveOneLabelOut(df.site)
scores = sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)
print scores.mean(), scores.std()
There is also a version that holds out $P$ labels, and lets you choose $P$:
In [37]:
cv = sklearn.cross_validation.LeavePLabelOut(df.site, p=2)
scores = sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)
print scores.mean(), scores.std()
What did I not mention that whole time? Leave-one-out cross-validation. That is because it takes forever.
In [38]:
%%time
n = len(y)
cv = sklearn.cross_validation.LeaveOneOut(n)
scores = sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)
print scores.mean(), scores.std()
Let's return to the stratified, repeated, 10-fold cross-validation approach:
In [39]:
%%time
scores = []
for rep in range(10):
print rep,
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True)
scores += [sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)]
print
In [40]:
np.mean(scores), np.std(scores)
Out[40]:
Sure, we can have a couple of minutes of break-time running this thing, but how do you think this is actually doing? Last week, we used a random prediction method as a baseline, and that is reasonable, but this week let's use something even simpler. Predicting a single class:
In [41]:
y_single = 'Diabetes'
np.mean(y == y_single)
Out[41]:
In [42]:
y_single = 'Stroke'
np.mean(y == y_single)
Out[42]:
In [43]:
y_single = 'Other'
np.mean(y == y_single)
Out[43]:
Uh, oh... we don't need any fancy classifiers to get 80% accuracy, and we don't need to wait around for five minutes to see how well that fancy classifier does. We could just say "other" and call it a day. To be a really fair (yet still damning) comparison, let's make a stratified, repeated 10-fold cross-validation of the baseline classifier that always says "other".
To do so, it will help you to know how to use a sklearn.cross_validation object a little bit more:
In [44]:
cv = sklearn.cross_validation.KFold(10, n_folds=3, shuffle=False)
for train, test in cv:
print train, test
We can use that to make a fair comparison of the decision tree to the baseline classifier that always predicts "Other":
In [45]:
scores = []
for rep in range(10):
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True)
for train, test in cv:
scores.append(np.mean(y[test] == 'Other'))
np.mean(scores)
Out[45]:
To make this really, really fair, we should use the same random splits in each experiment...
In [46]:
%%time
scores = []
for rep in range(10):
print rep,
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True, random_state=123456+rep)
scores += [sklearn.cross_validation.cross_val_score(clf, X, y, cv=cv)]
print
print np.mean(scores)
In [47]:
scores = []
for rep in range(10):
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True, random_state=123456+rep)
for train, test in cv:
scores.append(np.mean(y[test] == 'Other'))
np.mean(scores)
Out[47]:
But that is just practicing "defensive research", to prepare for the pickiest of critics...
What we really need to be thinking about now is what was called "cost-sensitive learning" in the Data Mining book.
In [48]:
df.Cause.value_counts()
Out[48]:
Make weights for examples proportional to inverse of counts:
In [49]:
weights = 1000. / df.Cause.value_counts()
weights
Out[49]:
In [50]:
sample_weight = np.array(weights[y])
What is in sample_weights?
In [51]:
import sklearn.metrics
y_pred = ['Other' for i in range(len(y))]
sklearn.metrics.accuracy_score(y, y_pred, sample_weight=sample_weight)
Out[51]:
Is that what you expected?
This is a better metric for our classification challenge, and we should use it when building the classifier and when measuring its performance:
In [52]:
%%time
scores = []
for rep in range(10):
print rep,
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True, random_state=123456+rep)
for train, test in cv:
clf = sklearn.tree.DecisionTreeClassifier()
clf.fit(X[train], y[train], sample_weight=sample_weight[train])
y_pred = clf.predict(X[test])
scores.append(sklearn.metrics.accuracy_score(y[test], y_pred, sample_weight=sample_weight[test]))
print
print np.mean(scores)
In [53]:
scores = []
for rep in range(10):
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True, random_state=123456+rep)
for train, test in cv:
y_pred = ['Other' for i in range(len(test))]
scores.append(sklearn.metrics.accuracy_score(y[test], y_pred, sample_weight=sample_weight[test]))
np.mean(scores)
Out[53]:
The moral is: ML methods can do something, but you have to be careful about it!
This relates to a subtle point, that I think may be somewhat unique to population health metrics, where we are interested in the population-level mean of the predictions more than the individual-level predictions themselves. I will see how far we make it through this exercise and decide if we have time to really dig into it.
In short, the percent of test set examples with underlying cause Diabetes is equal to the percent of the training set with this underlying cause. So it is not really out-of-sample, even though the actual examples are completely disjoint.
In [54]:
def csmf_accuracy(y_true, y_pred):
csmf_diff = 0
csmf_true_min = 1.
for j in ['Stroke', 'Diabetes', 'Other']:
csmf_true_j = np.mean(np.array(y_true) == j)
csmf_pred_j = np.mean(np.array(y_pred) == j)
csmf_diff += np.abs(csmf_true_j - csmf_pred_j)
if csmf_true_j < csmf_true_min:
csmf_true_min = csmf_true_j
#print csmf_true_j, csmf_pred_j, csmf_diff, csmf_true_min
return 1 - csmf_diff / (2 * (1 - csmf_true_min))
# test this function
y_true = ['Stroke', 'Diabetes', 'Other']
y_pred = ['Diabetes', 'Other', 'Stroke']
assert np.allclose(csmf_accuracy(y_true, y_pred), 1)
y_true = ['Stroke', 'Stroke', 'Stroke']
y_pred = ['Diabetes', 'Other', 'Diabetes']
assert np.allclose(csmf_accuracy(y_true, y_pred), 0)
The baseline predictor is pretty good:
In [55]:
scores = []
for rep in range(10):
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True, random_state=123456+rep)
for train, test in cv:
y_pred = ['Other' for i in range(len(test))]
scores.append(csmf_accuracy(y[test], y_pred))
np.mean(scores)
Out[55]:
But because the validation environment is leaking information, a different silly predictor is even better. I call this "Random-From-Train":
In [56]:
scores = []
for rep in range(10):
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True, random_state=123456+rep)
for train, test in cv:
y_pred = np.random.choice(y[train], size=len(test), replace=True)
scores.append(csmf_accuracy(y[test], y_pred))
np.mean(scores)
Out[56]:
To fix it, you must resample the cause distribution of the test dataset. Super-hard extra bonus homework: figure out how to do this.
In [60]:
scores = []
for rep in range(10):
print rep,
cv = sklearn.cross_validation.StratifiedKFold(y, n_folds=10, shuffle=True, random_state=123456+rep)
for train, test in cv:
clf = sklearn.tree.DecisionTreeClassifier()
clf.fit(X[train], y[train], sample_weight=sample_weight[train])
# resample y[test] here to have a distribution chosen uniformly at random
rows = []
for j in np.unique(y[test]):
test_rows_with_cause_j = np.where(y[test] == j)[0] # np.where returns a tuple, just to be annoying
rows += list(np.random.choice(test_rows_with_cause_j, size=1000))
X_test = X[rows]
y_pred = clf.predict(X_test)
scores.append(csmf_accuracy(y[rows], y_pred))
print
np.mean(scores)
Out[60]:
Did I promise a little more for this week in the syllabus?
Probability prediction:
In [61]:
clf.predict_proba(X[test[:10]])
Out[61]:
Maybe more convincing with Naive Bayes:
In [62]:
import sklearn.naive_bayes
In [63]:
clf = sklearn.naive_bayes.MultinomialNB()
clf.fit(X[train], y[train], )
Out[63]:
In [64]:
clf.predict_proba(X[test[:10]])
Out[64]:
Use np.round to make this look nicer:
In [65]:
np.round(clf.predict_proba(X[test[:10]]), 2)
Out[65]:
It is a bit hard to understand those probabilities. What do they mean?
In [66]:
y[test[:10]]
Out[66]:
Let us hope that the middle column is the probability of 'Other'. The column names must be stored in the clf instance somewhere:
In [67]:
clf.classes_
Out[67]:
Make that into a pretty DataFrame
In [69]:
t = pd.DataFrame(clf.predict_proba(X[test[:10]]), columns=clf.classes_)
np.round(t, 2)
Out[69]:
And what about numeric prediction? I'll give us a little something to show the mechanics, although it is a weird example. Can you figure out what g1_07a is? Hint: remember that cb is the codebook.
In [70]:
def float_or_100(x):
try:
return float(x)
except:
return 100.
y = np.array(df.g1_07a.map(float_or_100))
Here is some MSE cross-validation:
In [78]:
scores = []
cv = sklearn.cross_validation.KFold(len(y), n_folds=10, shuffle=True, random_state=123456)
for train, test in cv:
clf = sklearn.tree.DecisionTreeRegressor()
clf.fit(X[train], y[train])
y_pred = clf.predict(X[test])
scores.append(sklearn.metrics.mean_squared_error(y[test], y_pred))
print
print np.mean(scores)
I prefer root mean squared error, myself. I find it more interpretable:
In [79]:
np.mean(np.sqrt(scores))
Out[79]:
There is also a mean absolute error in sklearn.
In [80]:
sklearn.metrics.mean_absolute_error
Out[80]:
Can you make a custom error function that doesn't charge for the missing values we clipped to 100?
In [85]:
def rmse_when_not_100(a, b):
diff = a-b
rows = np.where((a != 100) & (b != 100))[0]
return np.sqrt(np.mean(diff[rows]**2))
In [86]:
cv = sklearn.cross_validation.KFold(len(y), n_folds=10, shuffle=True, random_state=123456)
for train, test in cv:
clf = sklearn.tree.DecisionTreeRegressor()
clf.fit(X[train], y[train], sample_weight=sample_weight[train])
y_pred = clf.predict(X[test])
scores.append(rmse_when_not_100(y[test], y_pred))
print np.mean(scores)
Is this any good? How can you tell?
Answer: compare to a baseline method that makes a numeric prediction without looking at the feature vector
In [87]:
scores = []
for val in np.unique(y):
if val != 100:
scores.append(rmse_when_not_100(y, val))
max(scores)
Out[87]:
Yes, about twice as good as random.