# Set up
import numpy as np
from __future__ import division # division
import pandas as pd
import seaborn as sns # for visualiation
from scipy.stats import ttest_ind # t-tests
import statsmodels.formula.api as smf # linear modeling
import statsmodels.api as sm
import matplotlib.pyplot as plt # plotting
import matplotlib
from sklearn import metrics
matplotlib.style.use('ggplot')
%matplotlib inline

# Load the data, replace strings as numeric
df = pd.read_csv('./data/payment-default.csv')
df.default = df.default.replace(['Yes', 'No'], [1, 0]).astype(int)
df.student = df.student.replace(['Yes', 'No'], [1, 0]).astype(int)

"""
Explanation: Classification with Logistic Regression
In today's exercise, we'll introduce logistic regression as a mechanism for performing classification tasks. This exercise walks through a conceptual example discussd in chapter 4 of this book, and the data was extracted from the ISLR R package (specifically, the Default dataset). The question of interest for this exercise is, 

Can we predict which individuals will default on their credit card payments?

Set Up (0 minutes)
End of explanation
"""


# What is the **default rate** in the dataset (# of defaults / total)
rate = df[df.default == 1].count()[0]/df.shape[0] * 100 # 3% default rate

# What is the distribution of balances (for those who default, and those who don't)?

# You may want to create subsets of the data for those who defaulted, and those who did not
defaulted = df[df.default == 1]
paid = df[df.default != 1]

# Draw histograms of the distribution (perhaps overlapping histogram on the same chart)
plt.figure(figsize=(20,10))
plt.hist(paid.balance, alpha=0.5, label='Paid')
plt.hist(defaulted.balance, alpha=0.5, label='Defaulted')
plt.legend(loc='upper right')
plt.show()

# What is the relationship between balalance and income (show a scatter-plot, different colors for default/no default)?
# Hint: http://stackoverflow.com/questions/21654635/scatter-plots-in-pandas-pyplot-how-to-plot-by-category
fig, ax = plt.subplots(figsize=(10, 9))
ax.scatter(paid.income, paid.balance, marker="+", label='Paid', s=50)
ax.scatter(defaulted.income, defaulted.balance, label='Defaulted', marker="+", s=50, c='red')
ax.legend(numpoints=1, loc='upper right')
ax.set_xlabel('Income')
ax.set_ylabel('Balance')
plt.show()

"""
Explanation: Data Exploration (15 minutes)
In this first section, you'll explore your dataset to get a better handle on the distributions. The data for this exercise has 10,000 observations (people), who did, or did not, default on their credit card payment. The data has the following information for each person:

default: This is our outcome of interest, and is binary (yes/no)
student: This variable indicates if each person is a students (yes/no)
balance: Current balance on the credit card (continuous)
income: Annual income of the individual (continuous)

In this section, you'll write the code necessary to answer the following questions:

What is the default rate in the dataset (# of defaults / total)
What is the distribution of credit card balances (for those who default, and those who do not)?
What is the relationship between balance and income (show a scatter-plot, different colors for default/no default)?
End of explanation
"""


# Fit a linear model of the dependent variable (default) on balance
lm = smf.ols(formula='default ~ balance', data=df).fit()

# Generate predictions using your model (this can be interpreted as probability of default)
df['linear_preds'] = lm.predict()

# Visualize the results: scatterplot of balance v.s. predcited default rate. Bonus: add actual default rate.
fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(paid.balance, paid.default, c='blue', label="Paid")
ax.scatter(defaulted.balance, defaulted.default, c='red', label="Defaulted")
ax.plot(df.balance, df.linear_preds, c='black', label="Linear Fit")
ax.legend(numpoints=1, loc='upper left')
ax.set_xlabel('Balance')
ax.set_ylabel('Probability of Default')
plt.show()

# Interpret the **coefficients** and assess the **model fit**
lm.summary()

"""
Explanation: List pertinent observations from the above analysis:


Very few people in the dataset defaulted (3%)
Income and balance seem unrelated
The distribution of balances is much higher for those who defaulted


A Linear Approach (15 minutes)
Before introducing a new estimation approach, we'll observe the limitations of using a linear estimation of a binary outcome. In this section, you'll do the following:

Fit a linear model of the dependent variable (default) on one indepdendent variable (balance)
Generate predictions using your model (this can be interpreted as probability of default)
Visualize the results (comparing balance to the predicted default rate)
Interpret the coefficients and assess the model fit
End of explanation
"""


# Use the `smf.glm` function to fit `default` to `balance` (hint: use the `binomial` family)
logistic_model = smf.glm(formula = 'default ~ balance', data=df, family=sm.families.Binomial()).fit()

# Generate a set of predicted probabilities using your model
df['logistic_preds'] = logistic_model.predict()

# Visualize the predicted probabilities across balances (bonus: add the observed values as well)
df = df.sort_values(by='balance') # sort values for graphic
fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(paid.balance, paid.default, c='blue', label="Paid")
ax.scatter(defaulted.balance, defaulted.default, c='red', label="Defaulted")
ax.plot(df.balance, df.logistic_preds, c='black', label="Logistic Fit")
ax.legend(numpoints=1, loc='upper left')
ax.set_xlabel('Balance')
ax.set_ylabel('Probability of Default')
plt.show()

"""
Explanation: What is your interpretation of the coefficients and your assessment of the model fit?

A \$1 increase in balance is associated with a .0001 increase in the probability of default (which is statistically significant). As you can see, the model fit is poor (this is also reflected in the RSquared value, indicating that little variation in the outcome is explained by the (linear) variation in the predictor). The model also produces results that don't make sense (i.e, default probabilities that are negative).

The Logistic Function (10 minutes)
For this section, you'll simply familiarize yourself with the Logistic Function. In order to perform classification, we want to explicitly model the probability of falling into each category. Given the limitations of linear regression, we want to find a functional form with the following property:

All values output by the function should fall between 0 and 1.

Because we are modeling conditional probability, we will annotate our formula using the following shorthand:
\begin{equation}
p(X)= Pr\left(Y=1|X\right) \
\end{equation}
\begin{equation}
p(balance)=Pr(default=1|balance)
\end{equation}
While there are many functions that meet this requirement, a commonly used formula (and the one that we will use for logistic regression) is the logistic function:
\begin{equation}
p(X) = \frac { { e }^{ { { B }{ 0 } } + { B }{ 1 }X } }{ 1+{ e }^{ { { B }{ 0 } } + { B }{ 1 }X}} 
\end{equation}
Note, with some manipulaiton, we can rearrange the above equation into the following equality:
\begin{equation}
\frac { p(X) }{ 1 - p(X) } = { e }^{ { { B }{ 0 } }+{ B }{ 1 }X } 
\end{equation}
As you may have noticed, the left hand side of the equation is the odds of observing our dependent variable given our independent variable(s). This is because the odds is the ratio of the probability of an event, $p(X)$, relative to the probability of not observing the event, $1 - p(X)$.
By simply taking the log of both sides of the above equation, we arrive at the following formula:
\begin{equation}
\log \left( { \frac { p(X) }{ 1-p(X) } }\right) ={ { B }{ 0 } }+{ B }{ 1 }X
\end{equation}
On the left-hand side of the equation, we are left with the log odds, or logit, that has a linear relationship with the predictor variables. The betas are estimated using maximum likelihood methods, which are beyond the scope of this course. However, the intuition of this approach is that you are "finding the set of parameters for which the probability of the observed data is greatest" (source). 
Logistic models will generate a set of probabilities of an observation being a 0 or a 1. Given the functional form, the default threshold for classification is .5 (i.e., if a probability is less than .5, the model predicts it is not a positive case). 
Running a Logistic Regression (15 minutes)
We'll continue to use the statsmodels.formula.api to write R style formulas in our analytical approaches. Logistic Regression is considered a Generalized Linear Model, which is described well on Wikipedia:

The generalized linear model (GLM) is a flexible generalization of ordinary linear regression that allows for response variables that have error distribution models other than a normal distribution.

This class of models extends the linear model in a variety of ways, while providing an API that resembles the linear approach. 
In this section, you should: 

Use the smf.glm function to fit default to balance (hint: see this example)
Generate a set of predicted probabilities using your model
Visualize the predicted probabilities across balances
End of explanation
"""


# Retrieve beta values from the model
logistic_model.summary()

"""
Explanation: While the simplicity of the predict() method makes it easy to overlook the interworkings of the model. You should note that the estimates are generated by plugging values of $X$ into this equation (with your estimated Beta values):
\begin{equation}
p(X) = \frac { { e }^{ { { B }{ 0 } } + { B }{ 1 }X } }{ 1+{ e }^{ { { B }{ 0 } } + { B }{ 1 }X}} 
\end{equation}
Interpreting Betas (15 minutes)
As demonstrated above, the logistic formula is linear in logit space. This means that changes in your variables are associated with a change in the log odds of your dependent variable. There are a few implications of this:

P-Values retain their same interpretation (probability of observing the relationship by chance), and can be used to confirm/reject the null hypothesis
Direction of betas (+/-) is associated with a corresponding increase/decrease in the log odds of the outcome. This can be interpreted as a positive / negative relationship. 

Units
Units of the betas are more difficult to interpret. Each unit increase in $X$ is associated with an increase in the log odds of the dependent variable $Y$. As you can see in the curve above, increases in log odds are multiplicative in nature. This means that a unit increase in $X$ corresponds to a different change in $Y$ depending on your current value of $X$ (see how the slope of the curve changes). However, if we exponentiate the beta value, we can also describe proportional increase in the odds (not just the increase in log odds). Consider the following formulas for calculating the log odds for ${X}{i}$ and ${X}{i - 1}$ (a unit increase in $X$):
\begin{equation}
LogOddsX_{i} = \log \left( { \frac { p({X}{i}) }{ 1-p({X}{i}) } }\right) ={ { B }{ 0 } }+{ B }{ 1 }{X}_{i}
\end{equation}
\begin{equation}
LogOddsX_{i-1} = \log \left( { \frac { p({X}{i-1}) }{ 1-p({X}{i-1}) } }\right) ={ { B }{ 0 } }+{ B }{ 1 }{X}_{i-1}
\end{equation}
Because of the linear relationship in logit space, we can then easily calculate the difference in log odds between ${X}{i}$ and ${X}{i + 1}$ as $B_{1}$:
\begin{equation}
\log \left( { \frac { p({X}{i}) }{ 1-p({X}{i}) } }\right) - \log \left( { \frac { p({X}{i-1}) }{ 1-p({X}{i-1}) } }\right) = {B}_{1}
\end{equation}
If we exponentiate each side of the equation, we can see that the proportional increase in odds betweeen $P(X_{i-1})$ and $P(X_{i})$ is simply ${e}^{{B}_{1}}$:
\begin{equation}
\frac { \left( { \frac { p({ X }{ i }) }{ 1-p({ X }{ i }) }  } \right)  }{ \left( { \frac { p({ X }{ i-1 }) }{ 1-p({ X }{ i-1 }) }  } \right)  } ={ e }^{ B_{ 1} }
\end{equation}
We can then interpret the exponentiated beta value as the observed multiplicative unit increase in odds in the dependent variable.
In the section below, you should get the summary of your logistic model to retrieve betas values from your model.
End of explanation
"""


# Convert balance to units of $100 and re-run the regression
tmp = df.copy()
tmp['balance'] = tmp['balance']/100
logistic_model2 = smf.glm(formula = 'default ~ balance', data=tmp, family=sm.families.Binomial()).fit()

logistic_model2.summary()
np.exp(.5499)

"""
Explanation: Using the beta values produced by the model, describe the relationship between balance and default (make sure to note the p-value, direction, and units of the coefficient):


P Value of .000 indicates that the relationship is statistically significant
Direction of the relationship indicates that increases in balance are associated with increases in the log odds of the probability of default
Units are more interpretable if we exponentiate our Beta value. The odds of defaulting increase by a factor of 1.006 (.6%) for dollar unit increase in balance.


A More Interpretable Model (10 minutes)
The proportional increase of .6% in the model above is a bit small to have a good intuitive sense of. However, we can easily change our units of analysis to make this more interpretable. In the section below, you should do the following:

Convert balance to units of \$100 and re-run the regression (hint, you'll need to create a new column or copy your dataframe)
Extract the beta value for balance and interpret how an increase of your balance (by \$100) is associated with the odds of defaulting
End of explanation
"""


# First, using a threshold of .5, use your model to predict a binary outcome (each case as 0 or 1)
df['default_preds'] = np.where(df.logistic_preds > .5, 1, 0)

# What is the accuracy of your model (how often does the prediction match the data)?
df[df.default_preds == df.default].count() # 97% accurate!

# What are the sensitivity, specificity, Type I error rate, and Type II error rates?
cm = metrics.confusion_matrix(df.default, df.default_preds)
true_pos = cm[1,1]
false_pos = cm[0,1]
true_neg = cm[0,0]
false_neg = cm[1,0]

# Sensitivity, specificity, and error rates
sensitivity = true_pos / (true_pos + false_neg)
specificity = true_neg / (true_neg + false_pos)
type_one_error = 1 - specificity
type_two_error = 1 - sensitivity
print(type_one_error, type_two_error)

"""
Explanation: What is the interpretation of your beta value from your regression?

Each unit (\$100) increase is associated with increasing the odds of defaulting by 73%. 

Evaluation Metrics (10 minutes)
In the section above, you've developed a strategy for interpreting the betas produced by your model. However, you also need to assess how well your model fits your data. In order to do this, we need to return to the purpose of our model. We want to predict if people will default on their loans. Before diving back into the code, we'll introduce a confusion matrix to represent the relationship between the data and our estimates.
Sensitivity and Specificity
The concepts of sensitivity and specificity express how well your model accurately predicts positive and negative classes. In order to compute these, it is easy to first build a confusion matrix to show the total number of cases that were accurately predicted (and those that were not).
<table>
  <tr>
    <td/>
    <td/>
    <td colspan="2"><em>Predicted Class</em></td>
  </tr>
  <tr>
    <td/>
    <td/>
    <td><strong>Negative (-)</strong></td>
    <td><strong>Positive (+)</strong></td>
    <td><strong>Total</strong></td>
  </tr>
  <tr>
    <td rowspan="2"><em>True Class</em></td>
    <td><strong>Negative (-)</strong></td>
    <td>True Negative (TN)</td>
    <td>False Positives (FP)</td>
    <td>Total non-cases (N)</td>
  </tr>
  <tr>
    <td><strong>Positive (+)</strong></td>
    <td>False Negatives (FN)</td>
    <td>True Positives (TP)</td>
    <td>Total cases (P)</td>
  </tr>
</table>

In the matrix above, the rows represent the data, while the columns represent the predictions. For example, the top left cell (True Negatives) represents the number of instances in which someone did not default, and the model predicted that they did not default. Here are the definitions of each cell:

True Positives: The instances in which your model predicts a case (and the case is true in the data)
False Positives: The instances in which your model predicts a case (and the case is not true in the data)
False Negatives: The instances in which your model does not predict a case (and the case is true in the data)
True Negatives: The instances in which your model does not predict a case (and the case is not true in the data)

Using those values, we can compute the sensitivity and specificity of our model:
\begin{equation}
Sensitivity = \frac { True Positives }{ True Positives+False Negatives } 
\end{equation}
\begin{equation}
Specificity = \frac { TrueNegatives }{ TrueNegatives+FalsePositives } 
\end{equation}
These two metrics lend important insights regarding the quality of our model:

Sensitivity, also referred to as the true positive rate, tells us, of all of the cases in the data, how many did we accurately predict? This indicates the model's ability to detect cases. In other words, how sensitively does the model pick up on cases?
Specificity, also referred to as the true negative rate, tells us, of all of the non-cases in the data, how many did we accurately predict? This indicates the model's ability to assign non-cases.

These metrics are directly used to calculate Type I and Type II error rate, which are analagous to Type I and Type II errors in statistical tests (incorrectly reject a true null hypothesis, incorrectly accept a false null hypothesis). 

Type I Error rate is the proportion of instances which are incorrectly classified as positive cases (relative to the total number of negative cases). It is calculated as $1-specificity$, or simply the false positives relative to the total non-cases in the data, $FP/N$.
Type II Error rate is the proportion of instances which are incorrectly classified as negative cases (relative to the total number of positive cases). It is calculated as $1-sensitivity$, or simply the false negatives relative to the total cases in the data, $FN/P$.

Metric Calculation (15 minutes)
In this section, you'll evaluate the fit of your model by calculating the following:

First, using a threshold of .5, use your model to predict a binary outcome (each case as 0 or 1)
What is the accuracy of your model (how often does the prediction match the data)?
What are the sensitivity, specificity, Type I error rate, and Type II error rates?
End of explanation
"""


# Generate data for the ROC curve using the `metrics.roc_curve` function
fpr, tpr, thresholds  = metrics.roc_curve(df.default, df.logistic_preds)

# Draw your ROC curve
plt.figure(figsize=(20,10))
plt.plot(fpr, tpr)
plt.plot(fpr, fpr, 'b--')
plt.axis([-.02,1.02, -.02,1.02])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')

# Calculate the area under your ROC curve using the metrics.roc_auc_curve function
metrics.roc_auc_score(df.default, df.logistic_preds)

"""
Explanation: What is your interpretation of these values?

The model is really good at predicting who will not default (specificity). It has a Type I error rate of .4%. However, it is really bad at predicting who will default. The Type II error rate is 70% (meaning that it does not discover 70% of cases!). Note, the accuracy seems quite impressive -- it can accurately interpret 95% of cases. However, recall that only 3% of cases defaulted. This means that a naive model that predicts no one defaults would have an accuracy of 97%. We can choose a specific threshold (default is .5) and assess the sensitivity and specificity of our model.

ROC Curve (10 minutes)
It is common to compare the true positive rate (sensitivity) to the false positive rate (1 - specificity) at each threshold for classification in an ROC Curve. This interactive visualization may help you better understand the relationship between thresholds and the ROC. In this section, you'll do the following:

Generate data for the ROC curve using the metrics.roc_curve function
Draw your ROC curve
Calculate the area under your ROC curve using the metrics.roc_auc_curve function
End of explanation
"""


