by Devon Bodey
To find sufficient data, I used Data.gov. Once I found a dataset that I found interesting and contained more than 100 data entries with one target variable, I then narrowed down the information by removing irrelevant information. The dataset I chose for this analysis was Accidental Drug Related Deaths 2012-June 2017. A copy of the original dataset can be found here To make the information more relevant and usable for this analysis, I removed columns such as "Case Number", "Death County", and "Location". After removing several variables, I was left with "Age", "Sex" and information indicating which drug lead to the overdose for each victim. Once I finalized cleaning the data, my spreadsheet looked like this. To be properly analyzed, the data needed to contain all numbers. Originally, the data consisted of "Y" and "NaN", in excel, I found and replaced all "Y" with "1" and "NaN" with "0". Furthermore, within the notebook, I coded "Male" as 1 and "Female" as "0" so that there were two opposite variables which could be compared and contrasted against the remaining information.
Training and Cross Validation was extremely difficult for me given the insufficient dataset that I originally chose to work with. Seeing as I chose a dataset which only included information regarding individuals who had overdosed, it was nearly impossible to form any algorithm to calculate the statistical change of a drug overdose based on age, sex, drug of choice, etc. This is something I did not realize initially. However, simply put, it is impossible to code an algorithm to calculate outcomes on this basis because there is no base information to 'train' off of. The data on its own inferred that all women and men were using drugs and that the drugs that were most common for overdoses were the most deadly or commonly used. However, we know that that information is not true. Or, if it is, this data alone is not enough to prove it.
In order to work on resolving this issue, and in hopes of at least having some data to analyze regardless of accuracy, a fake dataset was added into the code. This dataset mimicked the information contained in the original, however, this time the person did not overdose. This allowed a statistical analysis to be done on overdose rates.
To make the analyses accurate and useful in the real world, I would have had to input more information regarding census data for men and women of certain ages, and rates of use for the drugs being analyzed. I found useful census information here, and this website had strong statistical models on heroin use in the state of Connecticut.
Unfortunately, given the lack of accurate information, I attempted to use the fake dataset to mimic the percentage of men who use heroin by making it so that only 1/3 of the data resulted in overdoses. However, I quickly learned the troubles with that as well. That data would have insinuated that every male who uses heroin would overdose. Additionally, the overdose rate being so infrequent lead to recurring insufficient results.
The current data analysis provides that assuming the individual is a drug user, factoring in sex, age and drug of choice (Heroin, Fentanyl, Oxycodone, and Alcohol/EtOH), there is a 50% chance they will overdose. However, we know this information is inaccurate seeing as we coded the results to be 50% by adding the fake data (1,2). Yet, from this information, we can tell that out of those who have overdosed, the most common overdose is due to Heroin and the least likely overdose is from Alcohol/EtOH.
This model would provide a properly-formed prediction when fed properly-formed data. It is due to the data, or lack thereof that this model is not ready for real-world application. To be real world ready, the data would need to be bolstered by adding information regarding non-drug-users and information about drug-users who did not overdose. Additionally, for it to be truly accurate, the ratio of drug-users to non-drug-users would have to match the statistics in Connecticut, but specialized by drug. The ratio of Male:Female would also have to be accurate based on the Connecticut census, and additionally accurate for every age group.
For this model to be real-world ready, a lot of information would need to be added. Despite not having a practical statistical evaluation, I found this exercise to be extremely helpful in hilighting the possibilities of statistical machine evaluation through the use of algorithms.
import os
try:
inputFunc = raw_input
except NameError:
inputFunc = input
import pandas as pd
from pandas.tseries.holiday import USFederalHolidayCalendar
from pandas.tseries.offsets import CustomBusinessDay
import numpy as np
import seaborn as sns
from statsmodels.formula.api import ols
from sklearn import linear_model
from sklearn import metrics
from sklearn.linear_model import LogisticRegression
from patsy import dmatrices
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
import matplotlib.pyplot as plt
import random
# Custom functions
def evaluate(pred, labels_test):
acc = accuracy_score(pred, labels_test)
print ("Accuracey: %s"%acc)
tn, fp, fn, tp = confusion_matrix(labels_test, pred).ravel()
recall = tp / (tp + fp)
percision = tp / (tp + fn)
f1 = (2 / ((1/recall)+(1/percision)))
print ("")
print ("True Negatives: %s"%tn)
print ("False Positives: %s"%fp)
print ("False Negatives: %s"%fn)
print ("True Positives: %s"%tp)
print ("Recall: %s"%recall)
print ("Precision: %s"%percision)
print ("F1 Score: %s"%f1)
def plot_bound(Z_val,data,col1,col2,binary):
# Z-val equals "Yes" value. E.g., "Y" or "1".
# data equals df
# col1 and col2 defines which colums to use from data
# Plot binary decision boundary.
# For this, we will assign a color to each
# point in the mesh [x_min, m_max]x[y_min, y_max].
x_min = float(data.iloc[:,[col1]].min())-float(data.iloc[:,[col1]].min())*0.10
x_max = float(data.iloc[:,[col1]].max()+float(data.iloc[:,[col1]].min())*0.10)
y_min = 0.0;
y_max = float(training.iloc[:,[col2]].max())+float(training.iloc[:,[col2]].max())*0.10
h_x = (x_max-x_min)/100 # step size in the mesh
h_y = (y_max-y_min)/100 # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h_x), np.arange(y_min, y_max, h_y))
if binary == 1:
Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
Z = np.where(Z=="Y",1,0)
else:
Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.pcolormesh(xx, yy, Z)
plt.show()
Here we load the data we collected and get it all ready to feed to our statistical model(s). That is, we are trying to make a table with one target column and one or more features. Here I'm loading happiness.csv from: https://data.somervillema.gov/Happiness/Somerville-Happiness-Survey-responses-2011-2013-20/w898-3dfm Note: you can find information on the data elements at this link.
# Load and peek at your data. Change the file name as needed.
processed_data_df = pd.read_csv('overdose.csv')
processed_data_df.head()
# You can replace values in a column based on logic like so
# Note: I used the unique values found above to inform my logic.
# That is, I took the unique text lables and translated them into numbers.
# It's clear that different surveys had different buckets. So I probably
# sould limit myself to years using the same metrics, but for our purposes
# I'm just going to run with a quick and dirty translation.
processed_data_df.loc[processed_data_df['Sex'] == 'Male', 'Sex'] = 1
processed_data_df.loc[processed_data_df['Sex'] == 'Female', 'Sex'] = 0
processed_data_df.head()
# To make sure all of your columns are stored as numbers, use the pd.to_numeric method like so.
processed_data_df = processed_data_df.apply(pd.to_numeric, errors='coerce')
# errors='coerce' will set things that can't be converted to numbers to NaN
# so you'll want to drop these like so.
processed_data_df = processed_data_df.dropna()
processed_data_df.head()
# The second set will be for classifiers where the target is a class.
# Happiness
data_df = processed_data_df[[
'Sex',
'Age', 'Heroin','Oxycodone','EtOH','Fentanyl'
]].copy()
data_df.head()
data_df["OD"] = 1
data_df.head()
# not reall!
data_1 = data_df.copy()
data_2 = data_1.copy()
data_2["OD"] = 0
for i in range(1,2):
data_df = pd.concat([data_df, data_2],ignore_index=True)
data_df.head()
#data_df = data_df.dropna()
data_df.head()
data_df[data_df["OD"]==1]
data = data_df
holdout = data.sample(frac=0.2)
training = data.loc[~data.index.isin(holdout.index)]
print(len(training),len(holdout))
# Define the target (y) and feature(s) (X)
features_train = training.drop("OD", axis=1).as_matrix(columns=None)
labels_train = training["OD"].as_matrix(columns=None)
features_test = holdout.drop("OD", axis=1).as_matrix(columns=None)
labels_test = holdout["OD"].as_matrix(columns=None)
# What percentage of the time is target Y?
print("Percentage of 1s: %s\n"%(len(data[data["OD"]==1])/len(data)))
# Logistic Regression
model = LogisticRegression(fit_intercept = False, C = 1e9)
clf = model.fit(features_train, labels_train)
pred = clf.predict(features_test)
print("Logistic Regression")
evaluate(pred, labels_test)
print (clf.coef_)
from sklearn import tree
clf = tree.DecisionTreeClassifier(min_samples_split=40)
clf = clf.fit(features_train, labels_train)
pred = clf.predict(features_test)
print("\nDecision Tree")
evaluate(pred, labels_test)
from sklearn.ensemble import RandomForestClassifier
clf = RandomForestClassifier()
clf = clf.fit(features_train, labels_train)
pred = clf.predict(features_test)
print("Random Forest")
evaluate(pred, labels_test)
from sklearn.svm import SVC
clf = SVC(kernel="rbf",probability=True)
clf = clf.fit(features_train, labels_train)
pred = clf.predict(features_test)
print("SVM")
evaluate(pred, labels_test)