<jupyter_start><jupyter_text>Credit Card Fraud Detection  --- [Introduction](intro) [Objective](objective) [Exploratory Data Analysis (EDA)](eda)* [Understand - Structure & Variable Distribution](understand)* [Clean - Feature Engineer, Scale, Balance, & Outlier Handling](clean)* [Analyze Relationships - Correlation](relate) [Modeling](model)* [Logistic Regression](logistic)* [Gradient Boosting Machine](gbm)* [Boosted Decision Tree](bdt)* [Extreme Gradient Boost](xgb)* [Random Forest](forest)* [Neural Network](nn) [Conclusions](outro) [Citations](citations)---  IntroductionAccording to the Federal Bureau of Investigation, "Credit card fraud is the unauthorized use of a credit or debit card, or similar payment tool (ACH, EFT, recurring charge, etc.), to fraudulently obtain money or property. Credit and debit card numbers can be stolen from unsecured websites or can be obtained in an identity theft scheme."[$^{1}$](1)In the FTC's Consumer Sentinel Network Data Book 2019 it was noted that "Credit card fraud tops the list of identity theft reports in 2019. The FTC received more than 271,000 reports from people who said their information was misused on an existing account or to open a new credit card account."[$^{2}$](2) The Nilson Report, the leading global card and mobile payments trade publication, releases annual fraud statistics and indicated in the November 2019 issue that global losses had amounted to \\$27.85 billion. The United States which represents 21.54% of transactions accounted for \\$9.47 or 33.99% of the gross worldwide losses.[$^{3}$](3)  With data breaches becoming more common in recent years we can imagine that these numbers will continue to rise presenting an opportunity for credit card companies to take a proactive stance on detecting fraudulent charges quickly.[$^{4}$](4) We will explore various models utilizing a credit card transaction data set from Kaggle and evaluating which has the best performance metrics for classifying an incoming transaction as fraudulent or normal. The data set is provided by [The Machine Learning Group (MLG)](https://mlg.ulb.ac.be/wordpress/) of [Université Libre de Bruxelles (ULB)](https://www.ulb.be) and contains transactions made by european cardholders in September 2013.  ObjectiveThe primary goal of this notebook is to develop a highly performant model for detecting fraudulent credit card transactions that could be utilized by our client the credit card company. Along the way we will outline a generalized end to end process for completing a machine learning project from data load to model validation and insights derived from the ultimate solution.   Exploratory Data Analysis (EDA)First we need to explore the data set to **understand** our variables,**clean** our dataset, and **analyze relationships** between variables. Ultimately EDA should allow us to answer the following questions from <NAME>'s article [A Gentle Introduction to Exploratory Data Analysis](https://towardsdatascience.com/a-gentle-introduction-to-exploratory-data-analysis-f11d843b8184)[$^{5}$](5):* *What question(s) are you trying to solve (or prove wrong)?** *What kind of data do you have and how do you treat different types?** *What’s missing from the data and how do you deal with it?** *Where are the outliers and why should you care about them?** *How can you add, change or remove features to get more out of your data?*While we won't explicitly answer these questions in this paper we will approach understanding, cleaning, and analyzing with them in mind. Before we begin it is important to remember that except for `Time` and `Amount` all explanatory variables are transformed due to privacy reasons. The explanatory variables were transformed through [prinicpal component analysis or PCA](https://en.wikipedia.org/wiki/Principal_component_analysis) and no further information can be provided regarding features `V1, V2, ..., V28`.`Time` contains the seconds elapsed between each transaction and the first transaction in the dataset. The feature `Amount` is the transaction amount, this feature can be used for example-dependant cost-senstive learning. Finally the feature `Class` is the response variable and it takes value 1 in case of fraud and 0 otherwise.  Understand - Structure & Variable DistributionWe start by reading in the credit card transaction data set CSV file as a pandas [DataFrame](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html) object and look a sample of some records using the [sample](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.sample.html) method. Then begin exploring using the [describe](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.describe.html) method to obtain descriptive statistics about each of the available variables.[$^{6}$](6)<jupyter_code>import pandas as pd

pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)
df = pd.read_csv('Data/creditcard.csv')
df.sample(5)
df.describe()<jupyter_output><empty_output><jupyter_text>The data set consists of `284807` observations of `31` variables, `30` explanatory and `1` response. Our data set does not contain any missing values making [imputation](https://en.wikipedia.org/wiki/Imputation_(statistics)) irrelevant. It should be noted that "the problem of missing data is relatively common in almost all research and can have a significant effect on the conclusions that can be drawn from the data." [$^{7}$](7) Various imputation methods can be taken when less than 30% of observations are missing at random, however, for scenarios where more than 30% of observations are missing or observations are not missing at random it is best practice to remove the variable(s) in question.<jupyter_code>df.isnull().sum()<jupyter_output><empty_output><jupyter_text>Looking at the response variable `Class` we see the dataset is very imbalanced with only **0.17%** of observations result in a fraudulent transactions. This is to be expected, if fraudulent transactions were very common credit card companies would lose money and likely cease to operate. The topic of balancing will be explored the next section but lets continue reviewing the structure and distribution of our available variables.<jupyter_code>df['Class'].value_counts(normalize=True)*100
ax = df['Class'].value_counts(sort=False).plot.bar(
    title="Imbalanced Classes\n0 = Normal\n1 = Fraud", 
    color=['grey','red'],
    rot=0,
    figsize=(13,5)
)
for i, v in enumerate(df['Class'].value_counts()):
    ax.annotate(str(round(v,4)),(i, v), ha='center', va='bottom');<jupyter_output><empty_output><jupyter_text>`Time` appears to be cyclical with periods every approximately `86000` seconds (24 hours). We will execute simple feature engineering to break out the various time components, e.g. hour, min, second. If the data contained at least 2 weeks worth of observations we could possibly extract the day of week component as there is often a weekly seasonality in transaction or sales data.[$^{30}$](30) We will subsequently scale the new time component features to align with the PCA variables.<jupyter_code>df['Time'].plot.hist(
    title="Time Frequency Distribution",
    bins=50, 
    color='grey',
    figsize=(13,5)
);<jupyter_output><empty_output><jupyter_text>`Amount` is a highly skewed variable with maximum observed value of `25,691.16` however most transactions averaging at only `88`. Again because of the wide distribution of values we will scale the values down and align with the PCA variables. We should also note that `Amount` does not have any unit associated with it. Given the fact that the data was sourced in Europe which has many currencies and the fact that we will be scaling the variable later we will leave the numeric value unitless.<jupyter_code>df['Amount'].plot.hist(
    title="Amount Frequency Distibution",
    bins=100,
    color='grey',
    figsize=(13,5),
    density=True
);<jupyter_output><empty_output><jupyter_text>All of our PCA variables,`V1, V2, ..., V28`, as mentioned before have been transformed.<jupyter_code>df.drop(['Time','Amount','Class'],axis=1).boxplot(figsize=(13,5));
df.drop(['Time','Amount','Class'],axis=1).hist(
    figsize=(13,15), 
    ylabelsize=False,
    grid=False
);<jupyter_output><empty_output><jupyter_text>Looking closer we can see some variety in the distributions even though all the PCA variables have mean values of `0`. Distributions of variables `V4`, `V12`, and `V28` are shown below. `V4` shows a slightly right skewed, `V12` appears to be left skewed, and `V28` shows a very narrow distribution. We will look at these distributions with respect to our target variable `Class` later as it might be that the normal transactions follow a normal distribution and that the fraudulent transactions are what are skewing these plots.<jupyter_code>df.V4.hist(
    bins=50,
    figsize=(13,5),
    grid=False
);
df.V12.hist(
    bins=50,
    figsize=(13,5),
    grid=False
);
df.V28.hist(
    bins=50,
    figsize=(13,5),
    grid=False
);<jupyter_output><empty_output><jupyter_text>Before we proceed to clean our data set we will make a copy to compare with later.<jupyter_code>original = df.copy()<jupyter_output><empty_output><jupyter_text>Clean - Feature Engineer, Scale, Balance, & Outlier Handling Now we will revisit the variables `Time` and `Amount` to extract and scale them. [Feature scaling](https://en.wikipedia.org/wiki/Feature_scaling) is a set of methods used to normalize the range of explanatory variables.[$^{8}$](8) The primary methods are: * [rescaling (min-max normalization)](https://en.wikipedia.org/wiki/Feature_scalingRescaling_(min-max_normalization))* [mean normalization](https://en.wikipedia.org/wiki/Feature_scalingMean_normalization)* [standardization (z-score normalization)](https://en.wikipedia.org/wiki/Feature_scalingStandardization_(Z-score_Normalization))* [scaling to unit length](https://en.wikipedia.org/wiki/Feature_scalingScaling_to_unit_length)<NAME> does a great job outlining several reasons why scaling is desirable in support of some machine learning algorithms. In [All about Feature Scaling](https://towardsdatascience.com/all-about-feature-scaling-bcc0ad75cb35) Roy notes that "Neural network gradient descent converge much faster with feature scaling than without it."He continues, "Scaling is critical while performing **Principal Component Analysis (PCA)**. PCA tries to get the features with maximum variance, and the variance is high for high magnitude features and skews the PCA towards high magnitude features." [$^{9}$](9)Given the use of PCA in our dataset it seems highly advisable to use feature scaling for the `Time` and `Amount` variables. `Time` doesn't have any substantial outliers, however `Amount` is heavily right skewed resulting in a large number of outliers as seen in the [Matplotlib](https://matplotlib.org) boxplots below.[$^{10}$](10)<jupyter_code>import matplotlib.pyplot as plt

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(13,5))
axes[0].get_xaxis().set_visible(False)
axes[0].set_title('Amount')
axes[0].boxplot(df['Amount'])
axes[1].get_xaxis().set_visible(False)
axes[1].set_title('Time')
axes[1].boxplot(df['Time']);<jupyter_output><empty_output><jupyter_text>Feature scaling particularly benefits machine learning models because it allows for improved gradient descent requiring less iterations to arrive at minimum cost $ J(\theta)$. There are a number of great resources around cost functions, $ J(\theta)$, including [<NAME>'s Machine Learning](https://www.coursera.org/learn/machine-learning) course on Coursera. "Put simply, a cost function is a measure of how wrong the model is in terms of its ability to estimate the relationship between X and y."[$^{38}$](38) Gradient descent is an optimization algorithm that aims at finding the local or global minima of the aforementioned cost function.  Before we scale however we need to take a moment to work through some feature engineering to extract the various time components (e.g. hour, minute, second) from the `Time` variable<jupyter_code>df['hour'] = pd.to_timedelta(df['Time'], unit='s').dt.components.hours
df['minute'] = pd.to_timedelta(df['Time'], unit='s').dt.components.minutes
df['second'] = pd.to_timedelta(df['Time'], unit='s').dt.components.seconds<jupyter_output><empty_output><jupyter_text>We will be using RobustScaler from the [Scikit-learn](https://scikit-learn.org/stable/index.html) package because it is as the name suggests robust to outliers. According to the [sklearn.preprocessing.RobustScaler](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.RobustScaler.html) documentation, "This Scaler removes the median and scales the data according to the quantile range (defaults to IQR: Interquartile Range). The IQR is the range between the 1st quartile (25th quantile) and the 3rd quartile (75th quantile).Centering and scaling happen independently on each feature by computing the relevant statistics on the samples in the training set. Median and interquartile range are then stored to be used on later data using the transform method."[$^{11}$](11)<jupyter_code>from sklearn.preprocessing import RobustScaler

scaler = RobustScaler()
hours_scaled = scaler.fit_transform(df['hour'].values.reshape(-1,1))
minutes_scaled = scaler.fit_transform(df['minute'].values.reshape(-1,1))
seconds_scaled = scaler.fit_transform(df['second'].values.reshape(-1,1))
amount_scaled = scaler.fit_transform(df['Amount'].values.reshape(-1,1))
df.drop(['Time','Amount','hour','minute','second'], axis=1, inplace=True)
df.insert(0, 'hour_scaled', hours_scaled)
df.insert(1, 'minute_scaled', minutes_scaled)
df.insert(2, 'second_scaled', seconds_scaled)
df.insert(3, 'amount_scaled', amount_scaled)<jupyter_output><empty_output><jupyter_text>Next we will look to balance the dataset to avoid overfitting in our downstream models. "Many real world classification problems are represented by highly unbalanced data sets, in which, the number of samples from one class is much smaller than from another. This is known as class imbalance problem and is often reported as an obstacle to constructing a model that can successfully discriminate the minority samples from the majority samples."[$^{12}$](12)To address the class imbalance problem we will utilize a combination of Random Under Sampling and SMOTE for over sampling. Random Under Sampling, a technique where we identify _n_, the minimum observations between our response variable classes, then take a random sample of size _n_ from the class with more observations to have a resulting 50/50 ratio. It should be noted that this process can result in loss of information but depending on the use case may have good results.[$^{13}$](13) SMOTE, which stands for Synthetic Minority Over-sampling Technique, over-samples the minority class by taking each minority class sample and introducing synthetic examples along the line segments joining any/all of the _k_ minority class nearest neighbors. Synthetic samples are generated by taking the difference between the sample under consideration and its nearest neighbor, multiplying the difference by a random number between 0 and 1 then adding it to the sample under consideration. This technique ultimately forces the decision region of the minority class to become more general.[$^{25}$](25)  We happen to have **492 Fraud** transactions and **284315 Normal** transactions so we will use Random Under Sampling to sample 49200 observations from the normal subset and SMOTE to generate 49200 Fraud observations. df['Class'].value_counts()<jupyter_code>from imblearn.over_sampling import SMOTE
seed = 1618 # random state for reproducability 

fraud = df[df['Class']==1]
normal = df[df['Class']==0].sample(len(fraud)*100, random_state=seed)
df = pd.concat([normal,fraud])
sm = SMOTE(random_state=seed)
x = df.drop(['Class'], axis = 1)
y = df["Class"]
x_res, y_res = sm.fit_resample(x, y)
df = x_res
df['Class'] = y_res
df = df.sample(frac=1, random_state=seed)
ax = df['Class'].value_counts(sort=False).plot.bar(
    title="Balanced Classes\n0 = Normal\n1 = Fraud", 
    color=['grey','red'],
    rot=0,
    figsize=(13,5)
)
for i, v in enumerate(df['Class'].value_counts()):
    ax.annotate(str(round(v,4)),(i, v), ha='center', va='bottom');<jupyter_output><empty_output><jupyter_text>Next up in the clean phase is outlier handling. Outliers can come in a variety of forms, the most common are **Global Anomalies** a.k.a. point anomalies, **Contextual Anomalies** a.k.a. conditional anomalies, and **Collective Anomalies**. These outliers are often the result of data entry errors, measurement errors, or are simply natural outliers.[$^{14}$](14)Jim Frost in [Guidelines for Removing and Handling Outliers in Data](https://statisticsbyjim.com/basics/remove-outliers/) notes, "Sometimes it’s best to keep outliers in your data. They can capture valuable information that is part of your study area."[$^{15}$](15) Others such as <NAME> and <NAME> in [The effects of handling outliers on the performance of bankruptcy prediction models](https://doi.org/10.1016/j.seps.2018.08.004) indicate that "There seems to be consensus in the literature on the necessity of handling outliers, at the same time, it is not clear how to define extreme values to be handled in order to maximize the predictive power of models."[$^{16}$](16)Given the seemingly conflicting ideas around outlier handling we will opt to look at very extreme values and use imputation to replace the removed values which should provide us with overall improvements to our model accuracy metrics. While there are a number of different outlier detection techniques we will be using NumPy[$^{28}$](28) and the IQR method, specifically targeting the most extreme values (3X IQR) applying a univariate approach. For imputation we will leverage [IterativeImputer](https://scikit-learn.org/stable/modules/generated/sklearn.impute.IterativeImputer.htmlsklearn.impute.IterativeImputer) from Scikit-learn to replace the removed values. IterativeImputer uses "\[a\] strategy for imputing missing values by modeling each feature with missing values as a function of other features in a round-robin fashion", which was inspired by the R MICE package (Multivariate Imputation by Chained Equations). This multivariate approach for imputation typically results in a significantly lower MSE when compared to simple univariate mean imputation. For a summary of several outlier detection techniques including Z-score method, Winsorization method, Isolation Forest, and more see <NAME>'s [Outlier!!! The Silent Killer](https://www.kaggle.com/nareshbhat/outlier-the-silent-killer) notebook on Kaggle.[$^{17}$](17)<jupyter_code>import numpy as np

outliers = {}
for col in df.drop(['hour_scaled','minute_scaled','second_scaled','Class'], axis=1).columns:
    q25 = np.percentile(df[col].loc[df['Class'] == 1].values, 25) 
    q75 = np.percentile(df[col].loc[df['Class'] == 1].values, 75)
    iqr = q75 - q25
    upper = q75 + (iqr * 3)
    lower = q25 - (iqr * 3)
    outliers[col] = [x for x in df[col].loc[df['Class'] == 1].values 
                if x < lower or x > upper]
    df[col].replace(df[(df[col] > upper) | (df[col] < lower)].index, np.nan, inplace=True)
from sklearn.experimental import enable_iterative_imputer  
from sklearn.impute import IterativeImputer

imputer = IterativeImputer(max_iter=10, random_state=seed)
imp_arr = imputer.fit_transform(df)
df = pd.DataFrame(imp_arr, columns=df.columns)<jupyter_output><empty_output><jupyter_text>Below we compare boxplots from our original dataframe and the now cleaned dataframe. Notice the y-axis scale is significantly reduced in our cleaned dataframe and the boxes are actually visible for most of our variables.<jupyter_code>fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(13,10));
axes[0].set_title('Original');
original.drop(['Time','Amount','Class'],axis=1).boxplot(figsize=(13,5), ax=axes[0]);
axes[1].set_title('Cleaned');
df.drop(['hour_scaled',
         'minute_scaled',
         'second_scaled',
         'amount_scaled',
         'Class'],axis=1).boxplot(figsize=(13,5), ax=axes[1]);<jupyter_output><empty_output><jupyter_text>Analyze Relationships - Correlation We continue by looking at the correlation across our explanatory variables. "Correlation analysis is a statistical method used to evaluate the strength of relationship between two quantitative variables. A high correlation means that two or more variables have a strong relationship with each other, while a weak correlation means that the variables are hardly related."[$^{18}$](18) We will be using the Pearson correlation coefficient (PCC) as our statistic for measuring the relationships between our variables. PCC, also referred to as Pearson's _r_, is a statistic that measures the linear correlation between two variables _X_ and _Y_ and has a value between -1 and +1. Variables with a positive correlation (PCC > 0) behave in a manner in which the increase in one of the variables results in an increase in the other. The opposite can be said of variables with a negative correlation (PCC < 0), an increase in one of the variables results in a decrease in the other. Variables with PCC close to 0 have no linear correlation.[$^{19}$](19)Below is a heatmap generated using the [Seaborn](https://seaborn.pydata.org/index.html) data visualization library and reflecting the Pearson correlation coefficient between all pairs of variables.[$^{20}$](20)  **Dark blue** indicating strong positive correlation (PCC $\approx$ +1) and  **dark red** indicating strong negative correlation (PCC $\approx$ -1).<jupyter_code>import seaborn as sns

corr = df.corr()
fig = plt.figure(figsize=(13,13))
sns.heatmap(
    corr, 
    cmap="RdBu",
    xticklabels=corr.columns,
    yticklabels=corr.columns
);<jupyter_output><empty_output><jupyter_text>Listing the PCC for all explanatory variables with respect to our dependent variable `Class` in ascending order<jupyter_code>print(corr['Class'].sort_values())<jupyter_output>V14             -0.758997
V12             -0.682848
V10             -0.636053
V16             -0.597100
V9              -0.574197
V17             -0.564852
V3              -0.563319
V7              -0.478326
V18             -0.471480
V1              -0.425659
V6              -0.423333
V5              -0.382766
hour_scaled     -0.184260
V24             -0.092008
V13             -0.056718
V15             -0.030294
minute_scaled   -0.003443
V22             -0.000122
V23              0.000481
second_scaled    0.002549
amount_scaled    0.040536
V26              0.052131
V25              0.070954
V8               0.084352
V27              0.093551
V28              0.106097
V21              0.141852
V20              0.187488
V19              0.277831
V2               0.503695
V11              0.699367
V4               0.714225
Class            1.000000
Name: Class, dtype: float64<jupyter_text>Using the statsmodels[$^{32}$](32) library to do a basic logistic regression we can see all the variables are significant with _p-value_ < _0.05_ and the model has a resulting **_AIC_ = 46.95** which we will use as a baseline for comparing against a simplified logisitic regression model later.<jupyter_code>import statsmodels.api as sm
sm_model = sm.Logit(df['Class'],df.drop(['Class'], axis=1)).fit(disp=0)
sm_model.summary()<jupyter_output>/opt/anaconda3/lib/python3.7/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm<jupyter_text>We will create a simple function for investigating goodness of fit using traditional AIC statistic. This statistic is a small component for comparing models we will discuss others in the modeling section.<jupyter_code>def AIC(true,pred,cols):
    resid = true - pred
    sse = sum(resid**2)
    k = len(cols)
    AIC = 2*k - 2*np.log(sse)
    
    return AIC
y_val_base = sm_model.predict(df.drop(['Class'],axis=1))
AIC_base = AIC(df['Class'],y_val_base,df.drop(['Class'],axis=1).columns)
print('Baseline AIC: ', AIC_base)<jupyter_output>Baseline AIC:  46.95134448365124<jupyter_text>The variables with the lowest correlations (+/- 0.1) can be removed as they will not provide any valuable information in our models because the distributions are too similar between the two classes.<jupyter_code>irr_vars = list(corr['Class'][abs(round(corr['Class'].sort_values(),2))<=0.10].index)
df.drop(irr_vars, axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>The variables with largest negative correlations are `V14, V12, V10, V16, V3, V9` and the ones with the largest positive correlations are `V4, V11, V2, V19`. We will look at each of these next to understand how the transactions compare across the two classes.<jupyter_code>df.boxplot(['V14','V12'], by='Class', figsize=(13,5), grid=False)
df.boxplot(['V10','V16'], by='Class', figsize=(13,5), grid=False)
df.boxplot(['V3','V9'], by='Class', figsize=(13,5), grid=False);
df.boxplot(['V4', 'V11'], by='Class', figsize=(13,5), grid=False)
df.boxplot(['V2', 'V19'], by='Class', figsize=(13,5), grid=False);<jupyter_output><empty_output><jupyter_text>Another way to compare the distributions of the remaining variables with respect to `Class` is by using histogram plots with kernel density estimates (KDE) overlay. We have color coded normal transactions in  **gray**  and fraud transactions in  **red**  as we have done earlier in this notebook and as indicated in the legend for each plot.<jupyter_code>exp_vars = df.drop(['Class'], axis=1).columns
fig, axes = plt.subplots(nrows=len(exp_vars), ncols=1, figsize=(13,len(exp_vars)*5));
for var in exp_vars:
    sns.histplot(
        data=df, 
        x=var, 
        bins=50, 
        kde=True, 
        hue='Class', 
        palette=['grey','red'],
        ax=axes[list(exp_vars).index(var)]
    );<jupyter_output><empty_output><jupyter_text>ModelingBefore we begin generating various predictive models from our Binary Classification problem it is important to discuss how we can assess model performance. [Amazon Machine Learning Developer Guide](https://docs.aws.amazon.com/machine-learning/latest/dg/binary-classification.html) indicates "typical metrics are accuracy (ACC), precision, recall, false positive rate, and F1-measure". [$^{21}$](21)These evaluation metrics are calculated from four basic combinations of actual data category and assigned category:* **True Positive (TP)** - correct positive assignment* **True Negative (TN)** - correct negative assignment* **False Positive (FP)**  \[Type I Error\] - incorrect positive assignment   * **False Negative (FN)**  \[Type II Error\] - incorrect negative assignment    A confusion matrix, also known as an error matrix, is a 2x2 table that reflects these assignments. Each row of the matrix represents the instances in an actual class while each column represents the instances in a predicted class (or vice versa)[$^{29}$](29)Accuracy (ACC) is the fraction of correct assignment (TP + TN)/(TP + TN + FP + FN), precision is the fraction of TP/(TP + FP), recall a.k.a. sensitivity is the fraction TP/(TP + FN), false positive rate is FP/(TN + FP), and finally F1-measure or F1-score is the "harmonic mean of precision and recall"[$^{22}$](22)[$^{23}$](23) $\normalsize F_1= \frac{2}{recall^{-1} + precision^{-1}} = 2  \frac{precision * recall}{precision + recall}=\frac{TP}{TP + \frac{1}{2}(FP + FN)}$Depending on the use case one might be more interested in a model that performs well for a specific metric or subset of metrics over others. For credit card fraud it is important that we have both high precision and high recall. F1-measure is often considered to be the most relevant measure in fraud detection due to its combination of precision and recall into a single metric.[$^{24}$](24) Now that we know we need to evaluate our models to favor high recall or high sensitivity we need to split our data into train and test sets.<jupyter_code>from sklearn.model_selection import train_test_split

train_df, test_df = train_test_split(df, test_size=0.20, shuffle=True, random_state=seed)
x_train = train_df.drop(['Class'], axis = 1)
y_train = train_df["Class"]
x_test = test_df.drop(['Class'], axis = 1)
y_test = test_df["Class"]<jupyter_output><empty_output><jupyter_text>Logisitic RegressionLogistic model uses the logistic or inverse logit function to model a binary or dichotomous variable. Similar to linear regression, logistic regression fits the data to a line however unlike linear regression the line is not straight but rather an S shaped curved called **Sigmoid**. The S shaped curve is bound on the Y axis between 0 and 1 which fits our goal of binary classification. The sigmoid function (shown below) provides us with a probability of an observation belonging to one of the two classes. [$^{26}$](26) [$^{27}$](27)$\normalsize sigmoid(x)= \frac{1}{1+e^{-x}} = \frac{e^{x}}{e^{x}+1}$Representative plots for a linear regression and logistic regression can be seen below. **NOTE:** _We will be using a random search approach for selecting optimal model hyperparameters for all models shown below. Additionally, we will be looking at the delta (difference) of F1 scores between train and test sets. If F1 delta is > 0 then we can conclude the model does not overfit and can generalize for unseen observations. If F1 delta is < 0 then we would say the model is overfitting and does not generalize for unseen observations. For more details on hyperparameter tuning please refer to "Hypterparameter Tuning for Machine Learning Models" by <NAME> and "Random Search for Hyper-Parameter Optimization" by <NAME> et al. [$^{33}$](33)[$^{34}$](34) For more details on overfitting please refer to "The Problem of Overfitting" by <NAME>._ [$^{31}$](31)<jupyter_code>from sklearn.linear_model import LogisticRegression
from sklearn.metrics import recall_score, precision_score, f1_score, classification_report
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import uniform as sp_randFloat
from scipy.stats import randint as sp_randInt    
import time

params_lr = {
    'penalty' : ['l1', 'l2'],
    'C' : sp_randFloat(0, 4),
    'solver' : ['lbfgs',],
    'random_state': [seed]
}

lr = RandomizedSearchCV(
    LogisticRegression(), 
    param_distributions=params_lr, 
    cv=5, 
    n_iter=5, 
    n_jobs=-1 
)

t1 = time.time()
best_lr = lr.fit(x_train, y_train)
t2 = time.time()
y_pred_lr = best_lr.predict(x_test)
t3 = time.time()
y_val_lr = best_lr.predict(x_train)
t_fit_lr = t2 - t1
t_pred_lr = t3 - t2

print('The best parameters across ALL searched params:\n', best_lr.best_params_)
f1_test_lr = f1_score(y_test, y_pred_lr)
f1_train_lr = f1_score(y_train, y_val_lr)
f1_delta_lr = round(f1_test_lr-f1_train_lr,5)
print('F1 Score Delta: {}'.format(f1_delta_lr))
print(classification_report(y_test, y_pred_lr, target_names=['Normal', 'Fraud'], digits=5))
AIC_lr = AIC(y_test,y_pred_lr,x_train.columns)
print('LR AIC: ', AIC_lr)<jupyter_output>LR AIC:  25.91892721956809<jupyter_text>Using the same statsmodels approach as we did for the basic logistic regression in the [analyze relationships - correlation](relate) section in the two cells below the resulting **_AIC_ = _25.1_**. This AIC value is very similar to that of the AIC from the Scikit-learn logisitc regression model above which had the value of **_25.9_**. This demonstrates that both the statsmodels library and the Scikit-learn library have similiar outcomes and that our efforts for feature reduction through the correlation analysis was the correct approach as these models fit our data better than the basic logisitic regression.<jupyter_code>sm_model = sm.Logit(y_train, x_train).fit(disp=0)
sm_model.summary()
y_val_sm = sm_model.predict(x_test)
AIC_sm = AIC(y_test, y_val_sm, x_train.columns)
print('Simplified AIC: ', AIC_sm)<jupyter_output>Simplified AIC:  25.11541884638514<jupyter_text>Gradient Boosting Machine Weak learners or base models are the building blocks for designing more complex models in ensemble theory. Since most errors in modeling are sourced from variance, noise, and bias ensemble methods that focus on reducing these errors typically fall into three "meta-algorithms", boosting, bagging, and stacking. The next couple of models will focus on boosting a technique where learners are combined sequentially.[$^{35}$](35)[$^{36}$](36) There are different approaches to optimizing boosting algorithms, hence the variety of boosting models available, however that is beyond scope for this notebook. The first boosting model we will evaluate is GradientBoostingClassier sometimes referred to as Gradient Boosting Machine which operates in the purest sense of the boosting method building weak learners on the previous weak learner.[$^{37}$](37)<jupyter_code>from sklearn.ensemble import GradientBoostingClassifier

params_gbm = {
    'learning_rate': sp_randFloat(),
    'subsample'    : sp_randFloat(),
    'n_estimators' : sp_randInt(100, 1000),
    'max_depth'    : sp_randInt(4, 10),
    'random_state': [seed] 
}

gbm = RandomizedSearchCV(
    GradientBoostingClassifier(), 
    param_distributions=params_gbm, 
    cv=5, 
    n_iter=5, 
    n_jobs=-1
)

t1 = time.time()
best_gbm = gbm.fit(x_train, y_train)
t2 = time.time()
y_pred_gbm = best_gbm.predict(x_test)
t3 = time.time()
y_val_gbm = best_gbm.predict(x_train)
t_fit_gbm = t2 - t1
t_pred_gbm = t3 - t2

print('The best parameters across ALL searched params:\n', best_gbm.best_params_)
f1_test_gbm = f1_score(y_test, y_pred_gbm)
f1_train_gbm = f1_score(y_train, y_val_gbm)
f1_delta_gbm = round(f1_test_gbm-f1_train_gbm,5)
print('F1 Score Delta: {}'.format(f1_delta_gbm))
print(classification_report(y_test, y_pred_gbm, target_names=['Normal', 'Fraud'], digits=5))
AIC_gbm = AIC(y_test,y_pred_gbm,x_train.columns)
print('GBM AIC: ', AIC_gbm)<jupyter_output>GBM AIC:  35.204209454403255<jupyter_text>Boosted Decision Tree The next boosting model we will evaluate AdaBoostRegressor which uses adaptive boosting is an ensemble technique that uses the weighted sum of the various weak learners. The weak learners in this specific example are DecisionTreeRegressor instances a standard decision tree model.[$^{37}$](37)<jupyter_code>from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import AdaBoostRegressor

params_bdt = {
    'learning_rate': sp_randFloat(),
    'n_estimators' : sp_randInt(100, 1000),
    'random_state': [seed]
}

bdt = RandomizedSearchCV(
    AdaBoostRegressor(DecisionTreeRegressor()),
    param_distributions=params_bdt, 
    cv=5, 
    n_iter=5, 
    n_jobs=-1
)

t1 = time.time()
best_bdt = bdt.fit(x_train, y_train)
t2 = time.time()
y_pred_bdt = best_bdt.predict(x_test)
t3 = time.time()
y_val_bdt = best_bdt.predict(x_train)
t_fit_bdt = t2 - t1
t_pred_bdt = t3 - t2

print('The best parameters across ALL searched params:\n', best_bdt.best_params_)
f1_test_bdt = f1_score(y_test, y_pred_bdt)
f1_train_bdt = f1_score(y_train, y_val_bdt)
f1_delta_bdt = round(f1_test_bdt-f1_train_bdt,5)
print('F1 Score Delta: {}'.format(f1_delta_bdt))
print(classification_report(y_test, y_pred_bdt, target_names=['Normal', 'Fraud'], digits=5))
AIC_bdt = AIC(y_test, y_pred_bdt, x_train.columns)
print('BDT AIC: ', AIC_bdt)<jupyter_output>BDT AIC:  35.84111691664033<jupyter_text>Extreme Gradient Boost Our final boosting technique, described as an "efficient and scalable implementation of gradient boosting framwork by Friedman et al." [$^{45}$](45), eXtreme Gradient Boost (XGB) uses a "more regularized model formalization to control over-fitting, which gives it better performance" [$^{46}$](46) compared with Gradient Boosting Machine.<jupyter_code>import xgboost as xgb

param_xgb = {
    'n_estimators': sp_randInt(150, 1000),
    'learning_rate': sp_randFloat(0.01, 0.6),
    'subsample': sp_randFloat(0.3, 0.9),
    'max_depth': [3, 4, 5, 6, 7, 8, 9],
    'colsample_bytree': sp_randFloat(0.5, 0.9),
    'min_child_weight': [1, 2, 3, 4],
    'random_state': [seed]
}

xgb = RandomizedSearchCV(
    xgb.XGBClassifier(objective = 'binary:logistic'), 
    param_distributions = param_xgb,
    cv = 5,  
    n_iter = 5, 
    n_jobs = -1
)


t1 = time.time()
best_xgb = xgb.fit(x_train, y_train)
t2 = time.time()
y_pred_xgb = best_xgb.predict(x_test)
t3 = time.time()
y_val_xgb = best_xgb.predict(x_train)
t_fit_xgb = t2 - t1
t_pred_xgb = t3 - t2

print('The best parameters across ALL searched params:\n', best_bdt.best_params_)
f1_test_xgb = f1_score(y_test, y_pred_xgb)
f1_train_xgb = f1_score(y_train, y_val_xgb)
f1_delta_xgb = round(f1_test_xgb-f1_train_xgb,5)
print('F1 Score Delta: {}'.format(f1_delta_xgb))
print(classification_report(y_test, y_pred_xgb, target_names=['Normal', 'Fraud'], digits=5))
AIC_xgb = AIC(y_test, y_pred_xgb, x_train.columns)
print('XGB AIC: ', AIC_xgb)<jupyter_output>XGB AIC:  34.33357331188757<jupyter_text>Random Forest Random Forest models are in a different "meta-algorithm" class than the previous boosting models called bagging. Bagging is a method where weak learners are created in parallel, as opposed to boosting's sequential approach, and ultimately aggregated to create a kind of average of their individual outputs.[$^{37}$](37) Random Forests utilize parallel decision trees as shown in the image below.<jupyter_code>from sklearn.ensemble import RandomForestClassifier

params_rf = { 
    'n_estimators':  sp_randInt(100, 1000),
    'max_features': ['auto', 'sqrt', 'log2'],
    'max_depth': sp_randInt(4, 10),
    'bootstrap': [True, False],
    'random_state': [seed]
}

rf = RandomizedSearchCV(
    RandomForestClassifier(),
    param_distributions=params_rf, 
    cv=5, 
    n_iter=5, 
    n_jobs=-1
)

t1 = time.time()
best_rf = rf.fit(x_train, y_train)
t2 = time.time()
y_pred_rf = best_rf.predict(x_test)
t3 = time.time()
y_val_rf = best_rf.predict(x_train)
t_fit_rf = t2 - t1
t_pred_rf = t3 - t2

print('The best parameters across ALL searched params:\n', best_rf.best_params_)
f1_test_rf = f1_score(y_test, y_pred_rf)
f1_train_rf = f1_score(y_train, y_val_rf)
f1_delta_rf = round(f1_test_rf-f1_train_rf,5)
print('F1 Score Delta: {}'.format(f1_delta_rf))
print(classification_report(y_test, y_pred_rf, target_names=['Normal', 'Fraud'], digits=5))
AIC_rf = AIC(y_test, y_pred_rf, x_train.columns)
print('RF AIC: ', AIC_rf)<jupyter_output>RF AIC:  28.119657494559135<jupyter_text>Neural NetworkScikit-learn's implementation of an artificial neural network is the MLPClassifier. Artifical neural networks, often referred to simply as neural networks, are "a wide class of flexible nonlinear regression and discriminant models, data reduction models, and nonlinear dynamic systems" inspired by biological neural networks. Typically represented  as layers of nodes connected to one another like neurons connected through synapses. Each node representing some simple mathematical process of its input. The layers are commonly designated as the input layer (first layer), output layer (last layer), and hidden layer(s) (inner layer(s)). Most neural networks are fully connected from one layer to another with particular weights applied to the connections. [$^{39}$](39) [$^{40}$](40)Neural networks are amongst some of the most complex models in use at scale across industries. This complexity of course comes at a cost according to Geman et al. however recent papers have been investigating the U-shaped test error curve and refute this claim.[$^{41}$](41) [$^{42}$](42) Whether or not the bias-variance tradeoff is present in a neural network, one thing is certain, computational resources are required to execute and this can significantly increase the time for delivering a prediction. Our use case of credit card fraud may prohibit the use of this model but for the sake of evaluating the accuracy and its comparison to the previous models we will proceed.<jupyter_code>from sklearn.neural_network import MLPClassifier
params_nn = { 
    'hidden_layer_sizes':  [(sp_randInt.rvs(10,200,1),
                             sp_randInt.rvs(10,200,1),),
                            (sp_randInt.rvs(10,200,1),)],
    'activation': ['tanh', 'relu', 'logistic'],
    'solver': ['sgd', 'adam', 'lbfgs'],
    'alpha': sp_randFloat(0.0001, 0.05),
    'learning_rate': ['constant','adaptive'],
    'random_state': [seed]
}

nn = RandomizedSearchCV(
    MLPClassifier(),
    param_distributions=params_nn, 
    cv=5, 
    n_iter=5, 
    n_jobs=-1
)


t1 = time.time()
best_nn = nn.fit(x_train,y_train)
t2 = time.time()
y_pred_nn = best_nn.predict(x_test)
t3 = time.time()
y_val_nn = best_nn.predict(x_train)
t_fit_nn = t2 - t1
t_pred_nn = t3 - t2

print('The best parameters across ALL searched params:\n', best_nn.best_params_)
f1_test_nn = f1_score(y_test, y_pred_nn)
f1_train_nn = f1_score(y_train, y_val_nn)
f1_delta_nn = round(f1_test_nn-f1_train_nn,5)
print('F1 Score Delta: {}'.format(f1_delta_nn))
print(classification_report(y_test, y_pred_nn, target_names=['Normal', 'Fraud'], digits=5))
AIC_nn = AIC(y_test, y_pred_nn, x_train.columns)
print('NN AIC: ', AIC_nn)<jupyter_output>NN AIC:  34.33357331188757<jupyter_text>ConclusionsAfter evaluating all of our models and exclusively looking at their F1 scores, the **Boosted Decision Tree**, which leveraged an adaptive boosting technique performed the best. Interestingly if we were to use traditional statistics for evaluating model goodness of fit the **Logistic Regression** would be selected based on it minimizing the AIC. The LR model was also the only one that did not show signs of overfitting as seen with an F1 delta > 0 although the XGB model did have the smallest measure of overfitting.<jupyter_code>def bold(var):
    s = '\033[1m' + str(round(var,5)) + '\033[0m'
    return s

lr_str = 'Logistic Regression\nF1 score: {:.5f}\tF1 delta: {}\tAIC: {}\tTrain time: {}\tPredict time: {}\n'
gbm_str = 'Gradient Boosting Machine\nF1 score: {:.5f}\tF1 delta: {:.5f}\tAIC: {:.5f}\tTrain time: {:.5f}\tPredict time: {:.5f}\n'
bdt_str = 'Boosted Decision Tree\nF1 score: {}\tF1 delta: {}\tAIC: {:.5f}\tTrain time: {:.5f}\tPredict time: {:.5f}\n'
xgb_str = 'Extreme Gradient Boost\nF1 score: {:.5f}\tF1 delta: {}\tAIC: {:.5f}\tTrain time: {:.5f}\tPredict time: {:.5f}\n'
rf_str = 'Random Forest\nF1 score: {:.5f}\tF1 delta: {:.5f}\tAIC: {:.5f}\tTrain time: {:.5f}\tPredict time: {:.5f}\n'
nn_str = 'Neural Network\nF1 score: {:.5f}\tF1 delta: {:.5f}\tAIC: {:.5f}\tTrain time: {:.5f}\tPredict time: {:.5f}\n'

print(lr_str.format(f1_test_lr, bold(f1_delta_lr), bold(AIC_lr),bold(t_fit_lr), bold(t_pred_lr)))
print(gbm_str.format(f1_test_gbm, f1_delta_gbm, AIC_gbm,t_fit_gbm, t_pred_gbm))
print(bdt_str.format(bold(f1_test_bdt), f1_delta_bdt, AIC_bdt,t_fit_bdt, t_pred_bdt))
print(xgb_str.format(f1_test_xgb, f1_delta_xgb, AIC_xgb,t_fit_xgb, t_pred_xgb))
print(rf_str.format(f1_test_rf, f1_delta_rf, AIC_rf,t_fit_rf, t_pred_rf))
print(nn_str.format(f1_test_nn, f1_delta_nn, AIC_nn,t_fit_nn, t_pred_nn))<jupyter_output>Logistic Regression
F1 score: 0.94000	F1 delta: [1m0.00306[0m	AIC: [1m25.91893[0m	Train time: [1m5.85409[0m	Predict time: [1m0.0024[0m

Gradient Boosting Machine
F1 score: 0.99944	F1 delta: -0.00056	AIC: 35.20421	Train time: 1305.74339	Predict time: 0.26333

Boosted Decision Tree
F1 score: [1m0.99959[0m	F1 delta: -0.00041	AIC: 35.84112	Train time: 2333.66775	Predict time: 1.59219

Extreme Gradient Boost
F1 score: 0.99914	F1 delta: -0.00082	AIC: 34.33357	Train time: 29.10641	Predict time: 0.03720

Random Forest
F1 score: 0.98041	F1 delta: -0.00188	AIC: 28.11966	Train time: 425.23541	Predict time: 0.48133

Neural Network
F1 score: 0.99914	F1 delta: -0.00086	AIC: 34.33357	Train time: 367.89017	Predict time: 0.02989<jupyter_text>Given that the goal of our project was to determine the most performative model it is worthwhile for us to assess the speed of the models for both fitting and predicting. We should also note here that this dataset contains just two days of transactions and the true population size is unknown however we can assume that the population is considerably larger. This is important to know because some hyperparameter changes would likely be required for implementation, specifically the LR **solver** and **n_jobs** parameters. As noted in the Scikit-learn documentation _"For small datasets, ‘liblinear’ is a good choice, whereas ‘sag’ and ‘saga’ are faster for large ones."_<jupyter_code>fit_times = [t_fit_lr, t_fit_gbm, t_fit_bdt, t_fit_xgb, t_fit_rf, t_fit_nn]
pred_times = [t_pred_lr, t_pred_gbm, t_pred_bdt, t_pred_xgb, t_pred_rf, t_pred_nn]
tick_label=['LR','GBM','BDT','XGB','RF','NN']

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(13,5));
axes[0].set_title('Time to Fit');
axes[0].bar(x=[l for l in range(6)], height=fit_times, width=0.8, tick_label=tick_label)
for i, v in enumerate(fit_times):
    axes[0].annotate(str(round(v,2)),(i, v), ha='center', va='bottom')
axes[1].set_title('Time to Predict');
axes[1].bar(x=[l for l in range(6)], height=pred_times, width=0.8, tick_label=tick_label)
for i, v in enumerate(pred_times):
    axes[1].annotate(str(round(v,4)),(i, v), ha='center', va='bottom')
plt.show();<jupyter_output><empty_output><jupyter_text>While not in scope for our discussions there is a whole field of study around distributed computing and the architecture behind deploying ML solutions at scale that could shift practitioners decisions on the optimal solution.[$^{43}$](43) Based on our criterion for model comparison and evaluation the **Logistic Regression** has been identified as the best solution for our credit card fraud detection problem. Lets dive a bit deeper into this solution and take a look into the effects quantification and feature importance.<jupyter_code>final_model = LogisticRegression(
    C=1.3321057202868407, 
    penalty='l2', 
    solver='lbfgs',
    random_state=seed
).fit(x_train,y_train)
coefs = pd.DataFrame(final_model.coef_,columns=x_train.columns, index=['coefs'])
odds_ratios = pd.DataFrame(np.exp(final_model.coef_),columns=x_train.columns, index=['odds_ratio']) 
pd.concat([coefs,odds_ratios])<jupyter_output><empty_output><jupyter_text>While the principal components are unknown due to the anonymization process behind PCA we can see from the coefficients that the log-odds of a transaction being fraudulent is **1.063191** for a transaction having 1 unit increase of `V4` holding all other variables constant. Another way of describing this relationship is with odds ratio in which case we would state 1 unit increase of `V4` increases the odds of a transaction being fraudulent by **2.90** times. While odds ratio descriptions are significantly easier to understand than log-odds, this is still a bit ambiguous when compared to linear regression interpretations. Marginal effects are an alternative metric that describe the impact of the predictor on the outcome variable with the Average Marginal Effect (AME) representing the impact on the probability of the outcome.[$^{44}$](44) Using the marginal effect approach we can interpret 1 unit increase of `V4` increasing the probability of fraud by **73.8%**.<jupyter_code>from sklearn.inspection import plot_partial_dependence, partial_dependence

fig,ax= plt.subplots(figsize=(13, 5))
plot_partial_dependence(final_model, x_train, ['V4'], fig=fig)
V4_AME = round(np.mean(partial_dependence(final_model, x_train, [4])[0]),4)

print('Average Marginal Effect of V4: {}'.format(V4_AME))<jupyter_output>/opt/anaconda3/lib/python3.7/site-packages/sklearn/inspection/_plot/partial_dependence.py:323: FutureWarning: The fig parameter is deprecated in version 0.22 and will be removed in version 0.24
  FutureWarning)<jupyter_text>We wrap up the notebook with a quick review of feature importance. Scikit-learn's LogisticRegression model does not directly output feature importance like the ExtraTreesClassifier, RandomForestClassifier, or RecursiveFeatureElimination methods, however, one can calculate the relative feature importance as shown below. This approach uses the relative weights of the coefficients produced by the model. The top 5 important features are `V14, V12, V4, V10, V17` which unsurprisingly are also the 5 variables with the largest separation between the fraud and normal distributions as seen in the [histogram/KDE plots](kde) from the [analyze relationships - correlation section](relate).<jupyter_code>feature_importance = abs(final_model.coef_[0])
feature_importance = 100.0 * (feature_importance / feature_importance.max())
sorted_idx = np.argsort(feature_importance)
pos = np.arange(sorted_idx.shape[0]) + .75

fig = plt.figure(figsize=(13,5))
ax = fig.add_subplot(1, 1, 1)
ax.barh(pos, feature_importance[sorted_idx], align='center')
ax.set_yticks(pos)
ax.set_yticklabels(np.array(x_train.columns)[sorted_idx])
ax.set_title('Relative Feature Importance')
plt.show();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10165025938482905, 'fraction_numerical': 0.03257628364854377, 'mean_word_length': 5.28778733787081, 'pattern_counts': {'":': 0, '<': 140, '<?xml version=': 0, '>': 139, 'https://': 24, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '34', 'dir': 'jupyter-structured-clean-dedup', 'id': '521047', 'license': 'MIT', 'n_tokens_mistral': 16056, 'n_tokens_neox': 14632, 'n_words': 5367, 'path': 'credit-card-fraud.ipynb', 'repo_name': 'cristobalmitchell/credit-card-fraud', 'size': '49025'}