<jupyter_start><jupyter_text>CS-S109A Introduction to Data Science  Lecture 11: NNs and Visualizating Prediction Models**Harvard University****Summer 2020****Instructors:** <NAME>**Authors:** <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>---<jupyter_code>## RUN THIS CELL TO GET THE RIGHT FORMATTING 
import requests
from IPython.core.display import HTML
styles = requests.get("https://raw.githubusercontent.com/Harvard-IACS/2018-CS109A/master/content/styles/cs109.css").text
HTML(styles)<jupyter_output><empty_output><jupyter_text>Table of Contents  Learning Goals   Brief Tree Review       Architecture of Artificial Neural Networks (ANNs)   Variable Importances   Interpreting Prediction Models<jupyter_code>import pandas as pd
import sys
import numpy as np
import scipy as sp
import sklearn as sk
import matplotlib.pyplot as plt


from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
from sklearn import tree
from sklearn import ensemble


# Here are the decision trees
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import AdaBoostClassifier

import tensorflow as tf

print(tf.__version__)  # You should see a 2.0.0 here!

# sns.set(style="ticks")
# %matplotlib inline<jupyter_output>2.2.0<jupyter_text>Learning GoalsThis Jupyter notebook accompanies Lecture 11. By the end of this lecture, you should be able to:- have a better grasp of neural network archetecture- interpret a few different types of variable importances- interpret a prediction model by exploring the relationships of predictors with the response through prediction plots.  Part 0: Data WranglingFor this notebook we will be using the heart data set we've used all semester for performing classification:<jupyter_code>heart_df = pd.read_csv('../data/Heart.csv')
print(heart_df.shape)
heart_df.head()
heart_df.describe()
# Split into X and y
X = heart_df[['Age','Sex','ChestPain','RestBP','Chol','Fbs','RestECG','MaxHR','ExAng','Oldpeak','Slope','Ca','Thal']]
y = 1*(heart_df['AHD']=='Yes')
# fix categorical data types for maching learning methods...don't worry about the warning message

X['ChestPain']=X['ChestPain'].astype('category')
X['ChestPain']=X['ChestPain'].cat.codes

X['Thal']=X['Thal'].astype('category')
X['Thal']=X['Thal'].cat.codes

X.dtypes
# imputing zeroes for the missing values in `CA`

X['Ca']=X['Ca'].fillna(0)
X.describe()
# split into train and test

from sklearn.model_selection import train_test_split
itrain, itest = train_test_split(range(X.shape[0]), train_size=0.80)

X_train = X.iloc[itrain, :]
X_test = X.iloc[itest, :]
y_train = y.iloc[itrain]
y_test = y.iloc[itest]<jupyter_output><empty_output><jupyter_text>Part 1: tree-based modelsBelow `max_depth=3` and `max_depth=10` decision trees are fit.<jupyter_code>#fit the simple (depth = 3) decision tree classifier
dt3= tree.DecisionTreeClassifier(max_depth = 3)
dt3.fit(X_train,y_train)

#fit the an overfit (depth = 10) decision tree classifier
dt10 = tree.DecisionTreeClassifier(max_depth = 10)
dt10.fit(X_train,y_train)<jupyter_output><empty_output><jupyter_text>**Q1.1** Calculate the AUC on both train and test, and interpret the results<jupyter_code>######
#n Your code here
######

print("AUC on train for dt10:",sk.metrics.roc_auc_score(y_train,dt10.predict_proba(X_train)[:,1]))
print("AUC on test for dt10:",sk.metrics.roc_auc_score(y_test,dt10.predict_proba(X_test)[:,1]))

print("AUC on train for dt3:",sk.metrics.roc_auc_score(y_train,dt3.predict_proba(X_train)[:,1]))
print("AUC on test for dt3:",sk.metrics.roc_auc_score(y_test,dt3.predict_proba(X_test)[:,1]))
roc_dt3_train = sk.metrics.roc_curve(y_train,dt3.predict_proba(X_train)[:,1], drop_intermediate=False)
plt.plot(roc_dt3_train[0],roc_dt3_train[1]);
roc_dt3_test = sk.metrics.roc_curve(y_test,dt3.predict_proba(X_test)[:,1], drop_intermediate=False)
plt.plot(roc_dt3_test[0],roc_dt3_test[1]);

roc_dt10_train = sk.metrics.roc_curve(y_train,dt10.predict_proba(X_train)[:,1], drop_intermediate=False)
plt.plot(roc_dt10_train[0],roc_dt10_train[1]);
roc_dt10_test = sk.metrics.roc_curve(y_test,dt10.predict_proba(X_test)[:,1], drop_intermediate=False)
plt.plot(roc_dt10_test[0],roc_dt10_test[1]);<jupyter_output><empty_output><jupyter_text>**Solution:**The AUC is higher in train for the more complex tree with depth 10, but is lower in the test set.  This is a clear illustration of overfitting: the model performs will with the same set of data that it was fit on, but performs poorly with out-of-sample data. We continue fitting tree based models: first with a random forest, and then a boosted tree model.  Note: these are untuned.<jupyter_code>np.random.seed(109)
randomforest = RandomForestClassifier(n_estimators=100, max_features='sqrt', max_depth=10)
randomforest.fit(X_train,y_train);

adaboost = AdaBoostClassifier(
    base_estimator=DecisionTreeClassifier(max_depth=4),
    n_estimators=500,
    learning_rate=.75)
adaboost.fit(X_train,y_train);
# evaluating
print("AUC on train for randomforest:",sk.metrics.roc_auc_score(y_train,randomforest.predict_proba(X_train)[:,1]))
print("AUC on test for randomforest:",sk.metrics.roc_auc_score(y_test,randomforest.predict_proba(X_test)[:,1]))

print("AUC on train for adaboost:",sk.metrics.roc_auc_score(y_train,adaboost.predict_proba(X_train)[:,1]))
print("AUC on test for adaboost:",sk.metrics.roc_auc_score(y_test,adaboost.predict_proba(X_test)[:,1]))<jupyter_output>AUC on train for randomforest: 1.0
AUC on test for randomforest: 0.9308278867102396
AUC on train for adaboost: 1.0
AUC on test for adaboost: 0.9052287581699346<jupyter_text>**Q1.2** What would happen to the above AUC on train and test (random forest and adaboost) if the number of estimators (base trees) were increased for each? **Solution:**The AUC cannot get better on the train set for either model (1.0 is the maximum), but on the test there is a different story in expectation: the random forest would likely have a very similar AUC on test (may increase slightly) but the adaboost model would eventually get worse and worse as the model overfit more to the train set.   Part 2: NN modelBelow we build our first NN model for these data<jupyter_code>model_NN = tf.keras.models.Sequential([
    tf.keras.layers.Dense(100, input_shape=(pd.DataFrame(X_train).shape[1],), activation='relu'),
    tf.keras.layers.Dense(25, activation='tanh'),
    tf.keras.layers.Dense(1, activation='linear'),
])
model_NN.summary()<jupyter_output>Model: "sequential"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
dense (Dense)                (None, 100)               1400      
_________________________________________________________________
dense_1 (Dense)              (None, 25)                2525      
_________________________________________________________________
dense_2 (Dense)              (None, 1)                 26        
=================================================================
Total params: 3,951
Trainable params: 3,951
Non-trainable params: 0
_________________________________________________________________<jupyter_text>**Q2.1** How many hidden layers does this model have?  What should be the loss function for this model?  What is incorrect in the model architecture above?  Be sure to fix it... **Solution:**This model has 2 hidden layers (and one output layer).  Since this is a classification problem, the loss should be `binary_crossentropy` and the activation in the last (output) layer should be `sigmoid`, not `linear`.This model has 3951 parameters:1. The 100 nodes in layer-1 has 13 inputs (predictors) plus a bias each, so 100*(13+1)=1400 parameters in this layer.2. The 25 nodes in layer-2 has 100 inputs (nodes from layer-1) plus a bias each, so 25*(100+1)=2525 parameters in this layer.3. The single node in the output layer has 25 inputs (nodes from layer-2) plus a bias term, so 26 parameters in this layer.<jupyter_code>from numpy.random import seed
seed(109)
tf.random.set_seed(121)

# now fit the model, and evaluate:
model_NN = tf.keras.models.Sequential([
    tf.keras.layers.Dense(100, input_shape=(pd.DataFrame(X_train).shape[1],), activation='relu'),
    tf.keras.layers.Dense(25, activation='tanh'),
    tf.keras.layers.Dense(1, activation='sigmoid'),
])


model_NN.compile(optimizer='ADAM', loss='binary_crossentropy', metrics=['acc'])
history = model_NN.fit(X_train, y_train, epochs=100, batch_size=64, verbose=0)

print("AUC on train for NN_model:",sk.metrics.roc_auc_score(y_train,model_NN.predict(X_train)))
print("AUC on test for NN_model:",sk.metrics.roc_auc_score(y_test,model_NN.predict(X_test)))<jupyter_output>AUC on train for NN_model: 0.8735233516483517
AUC on test for NN_model: 0.863834422657952<jupyter_text>**Q2.2** Create a new NN model called `model_NN2` that improves upon the fixed model above.  Why do you suppose it is doing a better job?<jupyter_code>######
# your code here
######
# now fit the model, and evaluate:
model_NN2 = tf.keras.models.Sequential([
    tf.keras.layers.Dense(100, input_shape=(pd.DataFrame(X_train).shape[1],), activation='relu'),
    tf.keras.layers.Dense(25, activation='tanh'),
    tf.keras.layers.Dense(100, activation='relu'),
    tf.keras.layers.Dense(50, activation='tanh'),
    tf.keras.layers.Dense(1, activation='sigmoid'),
])


model_NN2.compile(optimizer='ADAM', loss='binary_crossentropy', metrics=['acc'])
history = model_NN2.fit(X_train, y_train, epochs=200, batch_size=128, verbose=0)

print("AUC on train for NN_model:",sk.metrics.roc_auc_score(y_train,model_NN2.predict(X_train)))
print("AUC on test for NN_model:",sk.metrics.roc_auc_score(y_test,model_NN2.predict(X_test)))<jupyter_output>AUC on train for NN_model: 0.9425824175824176
AUC on test for NN_model: 0.8006535947712419<jupyter_text>**Solution:**This model added two intermediate layer that repeated the hidden layers from before, and increased the number of epochs to 200.  This model may be performing better on the test set since it has more flexibility to approximate the true signal better; the first NN model may have been a little biased.  Part 3: Variable ImportanceBelow the variable importances are created for the 4 tree-based models:<jupyter_code>#Default Variable Importance

plt.figure(figsize=(24,6))
#plt.set_xticks()
#plt.set_xticklabels(X.columns)
num=10 

plt.subplot(1, 4, 1)
dt3_importances = dt3.feature_importances_
order = np.flip(np.argsort(dt3_importances))[0:num]
plt.barh(range(num),dt3_importances[order],tick_label=X.columns[order]);
plt.title("Relative Variable Importance for dt3")

plt.subplot(1, 4, 2)
dt10_importances = dt10.feature_importances_
order = np.flip(np.argsort(dt10_importances))[0:num]
plt.barh(range(num),dt10_importances[order],tick_label=X.columns[order]);
plt.title("Relative Variable Importance for dt10")

plt.subplot(1, 4, 3)
rf_importances = randomforest.feature_importances_
order = np.flip(np.argsort(rf_importances))[0:num]
plt.barh(range(num),rf_importances[order],tick_label=X.columns[order]);
plt.title("Relative Variable Importance for rf")

plt.subplot(1, 4, 4)
adaboost_importances = adaboost.feature_importances_
order = np.flip(np.argsort(adaboost_importances))[0:num]
plt.barh(range(num),adaboost_importances[order],tick_label=X.columns[order]);
plt.title("Relative Variable Importance for adaboost");<jupyter_output><empty_output><jupyter_text>**Q3.1** Interpret the plots above: why do they make sense?  How would the random forest variable imporance change if `max_features` was altered? **Solution:**The plots make sense: the decision tree of depth 3 has all of its splits defined by just a handful of variables as there are only 7 splits in this tree.  The tree of depth 10 has a lot of importance on just one variable (`ChestPain`, with relative importance of about 0.26), which defines the most influential splits on MSE when building the tree.  The random forest and adaboost both have a lot more uniformly spread-out importance, which makes sense since the random forest there are some base trees that do not consider `ChestPain` in the first splits, and in adaboost, the learning rate reduces the influence that any one single split contributes to the reduciton of MSE. Below we use the [`eli5`](https://eli5.readthedocs.io/en/latest/autodocs/sklearn.htmleli5.sklearn.permutation_importance.PermutationImportance) package to perform permutation importance for the random forest model.<jupyter_code>#pip install eli5
#permutation importance
import eli5
from eli5.sklearn import PermutationImportance
from eli5.permutation_importance import get_score_importances


perm = PermutationImportance(randomforest).fit(X_test, y_test)
#eli5.show_weights(perm,feature_names=X.columns)
print(X.columns)
eli5.show_weights(perm, feature_names = X_train.columns.tolist())<jupyter_output>/Users/kevinrader/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/deprecation.py:144: FutureWarning: The sklearn.metrics.scorer module is  deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.metrics. Anything that cannot be imported from sklearn.metrics is now part of the private API.
  warnings.warn(message, FutureWarning)
/Users/kevinrader/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/deprecation.py:144: FutureWarning: The sklearn.feature_selection.base module is  deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.feature_selection. Anything that cannot be imported from sklearn.feature_selection is now part of the private API.
  warnings.warn(message, FutureWarning)<jupyter_text>**Q3.2** How do the permutation importances compare to the default feature importance?  What is the difference in interpretation? `ChestPain` is still the most influential, but `Thal` was the second most important based on relative importance, but is near the bottom in the permutation importance.  This indicates (and can be explained by the dfact) that this variable is highly correlated with other predictors used in the model.<jupyter_code>#Note: eli5 does not behave well with Keras, by default.

perm = PermutationImportance(model_NN, random_state=1).fit(X_train,y_train)

#Note: eli5 does not behave well with Keras, by default.  But a user-


def score(model,X, y):
    y_pred = model.predict(X)
    return sk.metrics.accuracy_score(y, y_pred>0.5)

perm = PermutationImportance(model_NN, scoring=score, random_state=1).fit(X_train,y_train)

eli5.show_weights(perm, feature_names = X_train.columns.tolist())<jupyter_output><empty_output><jupyter_text>Part 4: Plotting PredictionsBelow we start to interpret relationships from various models based on the predictions from those models<jupyter_code>yhat_rf_train = randomforest.predict_proba(X_train)[:,1]
plt.scatter(X_train[['Age']],yhat_rf_train);
yhat_rf_test = randomforest.predict_proba(X_test)[:,1]
plt.scatter(X_test[['Age']],yhat_rf_test,marker='x');
plt.title("Predicted Probabilities vs. Age from the RF in train and test");<jupyter_output><empty_output><jupyter_text>**Q4.1** What does the above plot showing?  How can it be interpreted? **Solution:**This plot shows that there is not much directional relationship being predicted between Age and AHD in the random forest model.  It does suggest there are 2 subgroups being predicted in the train set (one gorup near a probability of zero and one group near 1) which suggests a categorical predictor is causing these differences in predicted values (likely `ChestPain` is the culprit). **Q4.1** Reproduce the above plot for your neural netowrk model.  How does it compare?  What does it say about Age's relationship with Cardiac Arrest?<jupyter_code>######
# Your code here
# Note: we also show it for a stronger predictor: MaxHR
######
plt.figure(figsize=(14,6))

plt.subplot(1, 2, 1)
yhat_nn_train = model_NN2.predict(X_train)
plt.scatter(X_train[['Age']],yhat_nn_train);
yhat_nn_test = model_NN2.predict(X_test)
plt.scatter(X_test[['Age']],yhat_nn_test,marker='x');
plt.title("Predicted Probabilities vs. Age from the RF in train and test");

plt.subplot(1, 2, 2)
yhat_nn_train = model_NN2.predict(X_train)
plt.scatter(X_train[['MaxHR']],yhat_nn_train);
yhat_nn_test = model_NN2.predict(X_test)
plt.scatter(X_test[['MaxHR']],yhat_nn_test,marker='x');
plt.title("Predicted Probabilities vs. MaxHR from the RF in train and test");<jupyter_output><empty_output><jupyter_text>**Q4.3** Fit a logistic regression to the predicted response from your NN model based on Age (in train).  Interpret the result<jupyter_code>from sklearn.linear_model import LogisticRegression

######
# your code here
######
logit = LogisticRegression(C=10000)
logit.fit(X_test[['Age']],np.ravel(1*(yhat_nn_test>0.5)))
logit.coef_<jupyter_output><empty_output><jupyter_text>**Solution:**The estimated $\hat{\beta}_1$ is estimated to be 0.10255, which suggests a positive relationship, on average.  A 1-year increase is associated with a $e^{0.10255}=1.108$ multiplicative increase (a 10.8% increase) in the odds of having a cardiac event on average based on the Neural Network model.   Below, several different plots are created:1. The predicted probabilities vs. age for any reasonable value of age at the mean values for the other predictors2. The predicted probabilties for each individual vs. Age (sometimes called profile plots) and the averaged individual probabilities vs. Age.3. The median of these individual predcited probability curves, along with the middle 95% ranges at any particular value of Age.<jupyter_code>means1 = X_train.mean(axis = 0)
#means1 =pd.Series(means)
means_df = (means1.to_frame()).transpose()
#df_repeated = pd.concat(means*3)
#print(df_repeated)
Ages = np.arange(np.min(X['Age']),np.max(X['Age']))
means_df  = pd.concat([means_df]*Ages.size,ignore_index=True)
means_df['Age'] = Ages
means_df.head()
#plots at means
yhat_nn = model_NN2.predict_proba(means_df)
plt.scatter(X_train['Age'],y_train)
plt.plot(means_df['Age'],yhat_nn,color="red")
plt.title("Predicted Probabilities vs. Age from NN in train");
#Plots for all pobservations.  And then averaged

means1 = X_train.mean(axis = 0)
#means1 =pd.Series(means)
means_df = (means1.to_frame()).transpose()
#df_repeated = pd.concat(means*3)
#print(df_repeated)
Ages = np.arange(np.min(X['Age']),np.max(X['Age']))
means_df  = pd.concat([means_df]*Ages.size,ignore_index=True)
means_df['Age'] = Ages
yhat_nns = []
for i in range(0,X_train.shape[0]):
    obs = X_train.iloc[i,:].to_frame().transpose()
    obs_df  = pd.concat([obs]*Ages.size,ignore_index=True)
    obs_df['Age'] = Ages
    yhat_nn = model_NN2.predict_proba(obs_df)
    yhat_nns.append(yhat_nn.transpose())
    plt.plot(obs_df['Age'],yhat_nn,color='blue',alpha=0.05)

plt.plot(obs_df['Age'],np.mean(yhat_nns,axis=0)[0],color='red',linewidth=2);
    
plt.ylim(0,1)
plt.title("Predicted Probabilities vs. Age from NN in train for all observations");
plt.plot(obs_df['Age'],np.median(yhat_nns,axis=0)[0],color='red');
plt.plot(obs_df['Age'],np.quantile(yhat_nns,q=.05,axis=0)[0],color='blue');
plt.plot(obs_df['Age'],np.quantile(yhat_nns,q=.95,axis=0)[0],color='blue');<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup