<jupyter_start><jupyter_text>Car Accident Severity Classification ModelingThis notebook is using data collected and cleaned in the `car_accidents_eda.ipynb` notebook, and will be using that data to create a classification model for predicting the severity of a car accident.  TABLE OF CONTENTS1. [Prepare data for modeling](Split-the-data-into-Training,-Validation,-and-Test-Sets)2. [Modeling](Modeling)    - [LG1 Model (Logistic Regression)](LG1-Model-(Logistic-Regression))    - [LG2 Model (Logistic Regression with Undersampling)](LG2-Model-(Logistic-Regression-with-Undersampling))    - [LR3 Model (Logistic Regression with SMOTE)](LR3-Model-(Logistic-Regression-with-SMOTE))    - [RF1 Model (Random Forest with all features)](RF1-Model-(Random-Forest-with-all-features))    - [RF2 Model (Random Forest with narrowed features)](RF2-Model-(Random-Forest-with-narrowed-features))    - [XGB1 Model (XGBoost with all features)](XGB1-Model-(XGBoost-with-all-features))    - [LG20 Model (2020-only data)](LG20-Model-(2020-only-data))    First, import libraries and modules needed ><jupyter_code># First, install these packages if needed by un-hashing and running:
# !pip install -U imbalanced-learn
# !pip3 install xgboost
%pylab inline
%config InlineBackend.figure_formats = ['retina']
%matplotlib inline

import pickle
import pandas as pd

import seaborn as sns

from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
import xgboost as xgb

from sklearn.metrics import confusion_matrix
from sklearn.metrics import roc_auc_score, roc_curve, classification_report
from sklearn.metrics import precision_score, recall_score, precision_recall_curve
from sklearn.metrics import f1_score
from sklearn.metrics import accuracy_score

from imblearn.over_sampling import SMOTE
from imblearn.under_sampling import RandomUnderSampler
from imblearn.over_sampling import ADASYN

import warnings

if not sys.warnoptions:
    warnings.simplefilter("ignore")<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Now import the dataset created and prepared for modeling ><jupyter_code># now re-open the dataset up again to continue work and analysis on it

with open('pickle/model_data2.pickle','rb') as read_file:
    model_data2 = pickle.load(read_file)<jupyter_output><empty_output><jupyter_text>- Split the data into Training, Validation, and Test Sets<jupyter_code>X = model_data2.drop(columns='Severity')
y = model_data2['Severity']
X_train_val, X_test, y_train_val, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify = y)
X_train, X_val, y_train, y_val = train_test_split(X_train_val,
                                                  y_train_val,
                                                  test_size=.25,
                                                  random_state=42,
                                                  stratify = y_train_val
                                                 )<jupyter_output><empty_output><jupyter_text>- Examine each variable's relationship to the targetStarting with pairplots for the numberical variables to the target variable (Severity)<jupyter_code># first look at relationship with variables/features that are numerical

sns.pairplot(model_data2[['Severity','Temperature(F)',
                          'Visibility(mi)','Wind_Speed(mph)']],
             hue='Severity', plot_kws=dict(alpha=.5));<jupyter_output><empty_output><jupyter_text>Takeaways:- No strong indiction of relationship between severity and these variables Before getting into modeling, create a function that will display a visual of the confusion matrix for the models' performances.  We'll then use this function for each model >  Modeling  LG1 Model (Logistic Regression)Before doing more advanced work, let's start with a simple Logistic Regression model and test out various combinations of variables to see what can get best results.  This will give us an indication of direciton to take feature selections ><jupyter_code>X_train.columns
std_scale = StandardScaler()

#features we want to use in this model
v1_features = ['Weather_Condition_Rain/Thunderstorms', 'Sunrise_Sunset',
               'Traffic_Signal','Weather_Condition_Snow/Hail',
               'Weather_Condition_Foggy/Hazy','Shut_Down',
               'Crossing','Temperature(F)',
               'Weather_Condition_Clear/Fair',
               'Junction'
              ]

#create new X set of v1 features for training and test set
X_train_v1 = X_train[v1_features]
X_val_v1 = X_val[v1_features]
X_test_v1 = X_test[v1_features]

# scale all of the features for training and test set
X_train_scaled_v1 = std_scale.fit_transform(X_train_v1)
X_val_scaled_v1 = std_scale.fit_transform(X_val_v1)
X_test_scaled_v1 = std_scale.fit_transform(X_test_v1)

# create training model using scaled features
LR1 = LogisticRegression(solver='liblinear').fit(X_train_scaled_v1, y_train)
# create predicted values of LR1 model

pred_train_lr1 = LR1.predict(X_train_scaled_v1)
pred_val_lr1 = LR1.predict(X_val_scaled_v1)
pred_test_lr1 = LR1.predict(X_test_scaled_v1)

# create list of predicted probabilities for LR1 model

predict_prob_lr1 = LR1.predict_proba(X_test_scaled_v1)
# validate score/accuracy for the LR1 model

print(f'Training Accuracy Score for LR1: {LR1.score(X_train_scaled_v1, y_train)}')
print(f'Validation Accuracy Score for LR1: {LR1.score(X_val_scaled_v1, y_val)}')
print(f'Test Accuracy Score for LR1: {LR1.score(X_test_scaled_v1, y_test)}')
# f1 score for the LR1 model

print(f'Training F1 Score for LR1: {f1_score(y_train, pred_train_lr1)}')
print(f'Validation F1 Score for LR1: {f1_score(y_val, pred_val_lr1)}')
print(f'Test F1 Score for LR1: {f1_score(y_test, pred_test_lr1)}')
# create simple confusion matrix

confusion_matrix(y_test, pred_test_lr1)
# find coefficients for this LR1 model

cofficients_lr1 = dict(zip(X_test_v1.columns, [x for nest in LR1.coef_ for x in nest]))
pd.DataFrame.from_dict(cofficients_lr1, orient='index', columns=['LR1 Model Coefficients'])
print('INTERCEPT for LR1 Model:')
print(LR1.intercept_)
# create metrics for ROC curve = false positive rate, true positive rate, thresholds

fpr_lr1, tpr_lr1, thresholds_lr1 = roc_curve(y_test, predict_prob_lr1[:,1])

# Precision = TP / (TP + FP)
# Recall = TP/P = True positive rate
# false positive rate = FP / true negatives = FP / (FP + TN) 

plt.plot(fpr_lr1, tpr_lr1,lw=2)
plt.plot([0,1],[0,1],c='violet',ls='--')
plt.xlim([-0.05,1.05])
plt.ylim([-0.05,1.05])


plt.xlabel('False positive rate')
plt.ylabel('True positive rate')
plt.title('ROC curve for car accident severity in LR1 model');
# plt.savefig('images/lr1_rocauc.svg', bbox_inches="tight")
print("ROC AUC score for LR1 model = ", roc_auc_score(y_test, predict_prob_lr1[:,1]))
# Precision & Recall at default threshold of 0.5

print("For LR1 model:")
print(f'''    Precision: {round(precision_score(y_test, pred_test_lr1),3)}, 
    Recall: {round(recall_score(y_test, pred_test_lr1),3)}''')<jupyter_output>For LR1 model:
    Precision: 0.599, 
    Recall: 0.089<jupyter_text>Takeaways of LM1 model:- Traffic Signal seems to be a strong feature.  Interestingly, it's negative, indicating a car accident is more likely to be less severe if there is a traffic signal present.  Perhaps this means that a person is more likely to be caucious/paying attention- Majority of signs/signals didn't influence model with any significance.  Weather conditions played a bigger role- Visiblity was removed because did not have that much influence and also there is co-variance risk with it and weather conditions as found in our earlier EDA- Very minor change of .0001 to F1 score and .00001 to accuracy, but technically did see that adding 'Shut_Down' -- i.e. whether the accident was during the period of stay-at-home orders improved the model.  This starts to provide answers to our greater question around impact of stay-at-home orders on accidents.   - LG2 Model (Logistic Regression with Undersampling)Now let's try running the same model as LR1, but with using Undersampling to balance out our classes and see if that helps improve scores ><jupyter_code># create undersampling -- this will give us 408k samples in each class

rus = RandomUnderSampler(random_state=42)

X_train_under, y_train_under = rus.fit_sample(X_train_scaled_v1,y_train)
X_val_under, y_val_under = rus.fit_sample(X_val_scaled_v1,y_val)
# build LM2 model using scaled & undersampled data

LR2 = LogisticRegression(solver='liblinear').fit(X_train_under, y_train_under)
# create predicted values of LR2 model

pred_train_lr2 = LR2.predict(X_train_scaled_v1)
pred_val_lr2 = LR2.predict(X_val_scaled_v1)
pred_test_lr2 = LR2.predict(X_test_scaled_v1)

# create list of predicted probabilities for LR2 model

predict_prob_lr2 = LR2.predict_proba(X_test_scaled_v1)
# validate score/accuracy for the LR2 model

print(f'Training Accuracy Score for LR2: {LR2.score(X_train_scaled_v1, y_train)}')
print(f'Validation Accuracy Score for LR2: {LR2.score(X_val_scaled_v1, y_val)}')
print(f'Test Accuracy Score for LR2: {LR2.score(X_test_scaled_v1, y_test)}')
# f1 score for the LR2 model

print(f'Training F1 Score for LR2: {f1_score(y_train, pred_train_lr2)}')
print(f'Validation F1 Score for LR2: {f1_score(y_val, pred_val_lr2)}')
print(f'Test F1 Score for LR2: {f1_score(y_test, pred_test_lr2)}')
# create simple confusion matrix

confusion_matrix(y_test, pred_test_lr2)
# find coefficients for this LR2 model

cofficients_lr2 = dict(zip(X_test_v1.columns, [x for nest in LR2.coef_ for x in nest]))
pd.DataFrame.from_dict(cofficients_lr2, orient='index', columns=['LR2 Model Coefficients'])
print('INTERCEPT for LR2 Model:')
print(LR2.intercept_)<jupyter_output>INTERCEPT for LR2 Model:
[-0.16570569]<jupyter_text>Takeaways:- Accuracy of model descreased.  Coefficients are essentially the same -- it's just the intercept that shifted to weight more towards the minority class.  This is not a model we want to move forward with. - LR3 Model (Logistic Regression with SMOTE)Now let's try running the same model as LR1, but with using SMOTE to balance out our classes and see if that helps improve scores ><jupyter_code># Create SMOTE sample. This will provide more emphasis 
#    on the regions where the class imbalance is greatest.

smote = SMOTE(random_state=42)

X_train_smote, y_train_smote = smote.fit_sample(X_train_scaled_v1,y_train)
X_val_smote, y_val_smote = smote.fit_sample(X_val_scaled_v1,y_val)
# build LM3 model using scaled & smote'd data

LR3 = LogisticRegression(solver='liblinear').fit(X_train_smote, y_train_smote)
# create predicted values of LR3 model

pred_train_lr3 = LR3.predict(X_train_scaled_v1)
pred_val_lr3 = LR3.predict(X_val_scaled_v1)
pred_test_lr3 = LR3.predict(X_test_scaled_v1)

# create list of predicted probabilities for LR3 model

predict_prob_lr3 = LR3.predict_proba(X_test_scaled_v1)
# validate score/accuracy for the LR3 model

print(f'Training Accuracy Score for LR3: {LR3.score(X_train_scaled_v1, y_train)}')
print(f'Validation Accuracy Score for LR3: {LR3.score(X_val_scaled_v1, y_val)}')
print(f'Test Accuracy Score for LR3: {LR3.score(X_test_scaled_v1, y_test)}')
# create simple confusion matrix

confusion_matrix(y_test, pred_test_lr3)<jupyter_output><empty_output><jupyter_text>Takeaways:- Accuracy of model descreased again and results are very similiar to basic undersampling results (LR2 model).  While we do get more correct in minority class, there are more misclassifications (false positive or false negative) then in LR1. - RF1 Model (Random Forest with all features)Now we'll test whether using a Random Forest model can improve our model.  We'll begin with using all features.<jupyter_code>RF1 = RandomForestClassifier(n_estimators=100)

RF1.fit(X_train_val, y_train_val)
# create predicted values of RF1 model

pred_train_RF1 = RF1.predict(X_train_val)
pred_test_RF1 = RF1.predict(X_test)

# create list of predicted probabilities for rf1 model

predict_prob_rf1 = RF1.predict_proba(X_test)<jupyter_output><empty_output><jupyter_text>Now look at validation score:<jupyter_code># K-Fold Class Validation Scores

scores_rf1 = cross_val_score(RF1, X_train_val, y_train_val, cv=5, scoring='f1')
print(scores_rf1)
print(round(scores_rf1.mean(),4))
# find F1 scores for training and test sets

print(f'Training F1 Score: {f1_score(pred_train_RF1, y_train_val)}')
print(f'Test F1 Score: {f1_score(pred_test_RF1, y_test)}')
# find score/accuracy for the RF1 model

print(f'Training Accuracy Score for RF1: {RF1.score(X_train_val, y_train_val)}')
print(f'Test Accuracy Score for RF1: {RF1.score(X_test, y_test)}')<jupyter_output>Training Accuracy Score for RF1: 0.7331904778175322
Test Accuracy Score for RF1: 0.6664706068985379<jupyter_text>-**Confusion Matrix** ><jupyter_code># create simple confusion matrix

confusion_matrix(y_test, pred_test_RF1)<jupyter_output><empty_output><jupyter_text>-**Feature Importance** ><jupyter_code>feat_importance_rf1 = dict(zip(X_test.columns, RF1.feature_importances_))
feat_importance_rf1<jupyter_output><empty_output><jupyter_text>Takeaways:- Seeing nice improvement in F1 score.  - Training and test scores have enough separation that there is likely some overfitting going on.  If decide to continue with this model, will do a GridSearchCV to refine model.- Dusty Weather appears to have a strong feature importance -- a feature that had minimal to no effect on the Logistic Regression models (returned to them to test and verified this) - RF2 Model (Random Forest with narrowed features)For some comparison and to see if can help with overfitting, we'll try a version of Random Forest with narrowed features to see if that helps performance.<jupyter_code>#create new X set of v1 features for combined train & validation sets
X_train_val_v1 = X_train_val[v1_features]
RF2 = RandomForestClassifier(n_estimators=100)

RF2.fit(X_train_val_v1, y_train_val)
# create predicted values of RF2 model

pred_train_RF2 = RF2.predict(X_train_val_v1)
pred_test_RF2 = RF2.predict(X_test_v1)

# create list of predicted probabilities for RF2 model

predict_prob_rf2 = RF2.predict_proba(X_test_v1)
# K-Fold Class Validation Scores

scores_rf2 = cross_val_score(RF2, X_train_val_v1, y_train_val, cv=5, scoring='f1')
print(scores_rf2)
print(round(scores_rf2.mean(),4))
# find F1 scores for training and test sets

print(f'Training F1 Score: {f1_score(pred_train_RF2, y_train_val)}')
print(f'Test F1 Score: {f1_score(pred_test_RF2, y_test)}')<jupyter_output>Training F1 Score: 0.20428206936419682
Test F1 Score: 0.19264717447987562<jupyter_text>Takeaways:   -  Significant dips on scores with the reduction of features.  This model is not a direction to continue with. - XGB1 Model (XGBoost with all features)Now we'll try an XGBoost model ><jupyter_code>xgb1 = xgb.XGBRegressor(n_estimators=100,
                        max_depth=10,
                        objective='binary:logistic',
                        learning_rate=.1, 
                        subsample=1,
                        min_child_weight=1,
                        colsample_bytree=.8
                       )

eval_set=[(X_train,y_train),(X_val,y_val)]

xgb1.fit(X_train, y_train, 
         eval_set=eval_set,
         eval_metric='error',
         early_stopping_rounds=30,
         verbose=True
        )
y_pred = xgb1.predict(X_test)
predictions = [round(value) for value in y_pred]

# evaluate predictions
accuracy = accuracy_score(y_test, predictions)
print("Accuracy: %.2f%%" % (accuracy * 100.0))<jupyter_output>Accuracy: 67.81%<jupyter_text>Takeaways:- Extremely similiar to the LM1 model on accuracy -- no significant improvements with this model - LR20 Model (2020-only data)To see if being in shut down plays a bigger role in predicting accidents in 2020 alone, we'll follow our exact method for the LM1 model, but only use 2020 samples. ><jupyter_code># create new dataframe with only 2020 samples

model_2020 = model_data2[model_data2.Year == 2020]

model_2020.shape<jupyter_output><empty_output><jupyter_text>With our filtered data using only 2020 samples we see we have ~205k samples. Now we'll create new training and test sets and run exactly as we did the LM1 model. ><jupyter_code>X_2020 = model_2020.drop(columns='Severity')
y_2020 = model_2020['Severity']
X_train_val_2020, X_test_2020, y_train_val_2020, y_test_2020 = train_test_split(
    X_2020, y_2020, test_size=0.2, random_state=42, stratify = y_2020)
X_train_2020, X_val_2020, y_train_2020, y_val_2020 = train_test_split(X_train_val_2020,
                                                                      y_train_val_2020,
                                                                      test_size=.25,
                                                                      random_state=42,
                                                                      stratify = y_train_val_2020
                                                                     )
std_scale = StandardScaler()

#features we want to use for this model (and used for LM1)
v1_features = ['Weather_Condition_Rain/Thunderstorms', 'Sunrise_Sunset',
               'Traffic_Signal','Weather_Condition_Snow/Hail',
               'Weather_Condition_Foggy/Hazy','Shut_Down',
               'Crossing','Temperature(F)',
               'Weather_Condition_Clear/Fair',
               'Junction'
              ]

#create new X set of v1 features for training and test set
X_train_v1_2020 = X_train_2020[v1_features]
X_val_v1_2020 = X_val_2020[v1_features]
X_test_v1_2020 = X_test_2020[v1_features]

# scale all of the features for training and test set
X_train_scaled_v1_2020 = std_scale.fit_transform(X_train_v1_2020)
X_val_scaled_v1_2020 = std_scale.fit_transform(X_val_v1_2020)
X_test_scaled_v1_2020 = std_scale.fit_transform(X_test_v1_2020)

# create training model using scaled features
LR20 = LogisticRegression(solver='liblinear').fit(X_train_scaled_v1_2020, y_train_2020)
# create predicted values of LR1_2020 model

pred_train_lr20 = LR20.predict(X_train_scaled_v1_2020)
pred_val_lr20 = LR20.predict(X_val_scaled_v1_2020)
pred_test_lr20 = LR20.predict(X_test_scaled_v1_2020)

# create list of predicted probabilities for LR1_2020 model

predict_prob_lr20 = LR20.predict_proba(X_test_scaled_v1_2020)
# validate score/accuracy for the LR1_2020 model

print(f'Training Accuracy Score for LR20: {LR1_2020.score(X_train_scaled_v1_2020, y_train_2020)}')
print(f'Validation Accuracy Score for LR20: {LR1_2020.score(X_val_scaled_v1_2020, y_val_2020)}')
print(f'Test Accuracy Score for LR20: {LR1_2020.score(X_test_scaled_v1_2020, y_test_2020)}')
# f1 score for the LR1 model

print(f'Training F1 Score for LR20: {f1_score(y_train_2020, pred_train_lr20)}')
print(f'Validation F1 Score for LR20: {f1_score(y_val_2020, pred_val_lr20)}')
print(f'Test F1 Score for LR20: {f1_score(y_test_2020, pred_test_lr20)}')
# create simple confusion matrix

print(f'''
Confusion matrix for LR20 Model: 
{confusion_matrix(y_test_2020, pred_test_lr20)}
''')

# find coefficients for this LR1 model

cofficients_lr20 = dict(zip(X_test_v1_2020.columns, [x for nest in LR20.coef_ for x in nest]))
pd.DataFrame.from_dict(cofficients_lr20, orient='index', columns=['LR20 Model Coefficients'])
print('INTERCEPT for LR20 Model:')
print(LR20.intercept_)
# create metrics for ROC curve = false positive rate, true positive rate, thresholds

fpr_lr20, tpr_lr20, thresholds_lr20 = roc_curve(y_test_2020, predict_prob_lr20[:,1])

# Precision = TP / (TP + FP)
# Recall = TP/P = True positive rate
# false positive rate = FP / true negatives = FP / (FP + TN) 

plt.plot(fpr_lr20, tpr_lr20,lw=2)
plt.plot([0,1],[0,1],c='violet',ls='--')
plt.xlim([-0.05,1.05])
plt.ylim([-0.05,1.05])


plt.xlabel('False positive rate')
plt.ylabel('True positive rate')
plt.title('ROC curve for car accident severity in LR20 model');
# plt.savefig('images/lr1_rocauc.svg', bbox_inches="tight")
print("ROC AUC score for LR20 model = ", roc_auc_score(y_test_2020, predict_prob_lr20[:,1]))<jupyter_output>ROC AUC score for LR20 model =  0.6630394983825849
{'fraction_non_alphanumeric': 0.08639298211029521, 'fraction_numerical': 0.04213690798876349, 'mean_word_length': 4.712837837837838, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 65, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '394102', 'license': 'AAL', 'n_tokens_mistral': 7367, 'n_tokens_neox': 6690, 'n_words': 1860, 'path': 'modeling_and_analysis/car_accidents_model.ipynb', 'repo_name': 'cecann10/covid-impact-car-accidents', 'size': '20291'}