#!/usr/bin/env python
# coding: utf-8

# # Logistic Regression

# In[1]:


import pandas as pd
import numpy as np

pd.options.display.max_columns=None
pd.options.display.max_rows=None


# ## 알고 가야 할 것 들
# 
# ### 1. Evidence of Weight
# ### 2. Information Value
# ### 3. VIF
# ### 4. C 통계량
# ### 5. AIC
# ### 6. ROC, AUC

# ## 데이터 셋을 로딩합니다.

# In[2]:


path = 'https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data'
dataset = pd.read_csv(path, delimiter=' ', header=None)


# In[4]:


dataset.shape


# In[5]:


dataset.head()


# In[6]:


COL = [
    'Status_of_existing_checking_account', 
    'Duration_in_month',
    'Credit_history',
    'Purpose', 
    'Credit_amount', 
    'Savings_account_bonds', 
    'Present_employment_since', 
    'Installment_rate_in_percentage_of_disposable_income',
    'Personal_status_and_sex',
    'Other_debtors_guarantors',
    'Present_residence_since', 
    'Property', 
    'Age_in_years', 
    'Other_installment_plans', 
    'Housing',
    'Number_of_existing_credits_at_this_bank',
    'Job',
    'Number_of_people_being_liable_to_provide_maintenance_for',
    'Telephone',
    'foreign_worker',
    'Target'
]


# In[7]:


dataset.columns = COL


# In[9]:


# 총 20개의 설명 변수와 1개의 종속 변수
dataset.head()


# Label
#  - 1 : Good (갚을 수 있다 / 신용도가 좋다) $\to$ 0으로 변경
#  - 2 : Bad (갚을 수 없다 / 신용도가 나쁘다) $\to$ 1로 변경

# - Ham(0) / Spam(1)
# - 정상사람(0) / 암환자(1)
# - 이벤트가 발생하지 않음(0) / 이벤트가 발생함(1)
# - 귀무가설(0) / 대립가설(1)

# In[10]:


# Label을 1/2 -> 0/1로 변경
dataset['Target'] = dataset['Target'] - 1


# ## 1. Information Value (정보 가치)
# 모델에서 변수의 사용유무를 판단하는 feature selection에서 유용한 방법입니다. 주로 모델을 학습하기전 첫 단계에서 변수들을 제거하는 데 사용합니다. 최종 모델에서는 대략 10개 내외의 변수를 사용하도록 합니다. IV와 WOE 신용채무능력이 가능한(good) 고객과 불가능한(bad) 고객을 예측하는 로지스틱 회귀 모델링과 밀접한 관계가 있습니다. 신용 정보 관련분야에서는 good customer는 부채를 갚을 수 있는 고객, bad customer는 부채를 갚을 수 없는 고객을 뜻합니다. 일반적으로 이야기할 때는 good customer는 non-events를 의미하고 bad customer는 events를 의미합니다.
# 
# > 신용 관련 분야
# 
# $
# \begin{align}
# &\text{WOE} = ln{\frac{\text{distribution of good}}{\text{distribution of bad}}} \\
# &\text{IV} = \sum{(\text{WOE} \times (\text{distribution of good} - \text{distribution of bad}))}\\
# \end{align}
# $
# 
# >일반적
# 
# $
# \begin{align}
# &\text{WOE} = ln{\frac{\text{distribution of non-events}}{\text{distribution of events}}} \\
# &\text{IV} = \sum{(\text{WOE} \times (\text{distribution of non-events} - \text{distribution of events}))}\\
# \end{align}
# $
# 
# 
# - 구간을 나눌 때, `pd.qcut`을 사용하였습니다. bin의 갯수만큼 값의 구간을 나누고 값을 해당 구간에 mapping 시켜줍니다.
# - Information value를 측정할 때, EOW가 inf, -inf가 나올 수 있습니다. 분포값이 0이 될 수 있기 때문입니다. 이 때는 EOW를 0으로 변경합니다.
# 
# 
# |<center>Information Value</center>|<center>예측력</center>|
# |:-----------:|:------------------------------:|
# |<center>0 to 0.02</center>|<center>무의미</center>|
# |<center>0.02 to 0.1</center>|<center>낮은 예측</center>|
# |<center>0.1 to 0.3</center>|<center>중간 예측</center>|
# |<center>0.3 to 0.5</center>|<center>강한 예측</center>|
# |<center>0.5 to 1</center>|<center>너무 강한 예측(의심되는 수치)</center>|

# 예제)
# 
# <img src="./images/information_value.png" style="height:400px width:400px">

# In[83]:


max_bin = 10  # 전체 데이터의 예측력에 해를 가하지 않는 한에서 구간을 테스트하였습니다.
def calc_iv(df, col, label, max_bin = max_bin):
    """IV helper function"""
    bin_df = df[[col, label]].copy()
    # Categorical column
    if bin_df[col].dtype == 'object':
        bin_df = bin_df.groupby(col)[label].agg(['count', 'sum'])
    # Numerical column
    else:
        bin_df.loc[:, 'bins'] = pd.qcut(bin_df[col].rank(method='first'), max_bin)
#         bin_df.loc[:, 'bins'] = pd.cut(bin_df[col], max_bin)
        bin_df = bin_df.groupby('bins')[label].agg(['count', 'sum'])
    
    bin_df.columns = ['total', 'abuse']
    bin_df['normal'] = bin_df['total'] - bin_df['abuse']
    bin_df['normal_dist'] = bin_df['normal'] / sum(bin_df['normal'])
    bin_df['abuse_dist'] = bin_df['abuse'] / sum(bin_df['abuse'])
    bin_df['woe'] = np.log(bin_df['normal_dist'] / bin_df['abuse_dist'])
    bin_df['iv'] = bin_df['woe'] * (bin_df['normal_dist'] - bin_df['abuse_dist'])
    
    bin_df.replace([np.inf, -np.inf], 0, inplace=True)
    bin_df = bin_df[bin_df['total'] > 0]
    iv_val = sum(filter(lambda x: x != float('inf'), bin_df['iv']))
    
    return bin_df, col, iv_val


# In[84]:


ch_df, ch, ch_i_val = calc_iv(dataset,'Credit_history', 'Target')
ch_df


# In[85]:


# 중간 정도의 예측 능력
print('information value', ch_i_val)


# In[86]:


dim_df, dim, dim_i_val = calc_iv(dataset,'Duration_in_month', 'Target')
dim_df


# In[87]:


# 높은 예측 능력
dim_i_val


# In[88]:


# 함수를 만들어서 전체 iv를 살펴보자
col_iv = {}
for col in [idx for idx in dataset.columns.tolist()]:
    if col == 'Target':
        continue
    _, col, iv = calc_iv(dataset, col, 'Target')
    col_iv[col] = iv


# In[47]:


col_iv


# In[89]:


import operator
candidates = sorted(col_iv.items(), key=operator.itemgetter(1), reverse=True)
display(candidates)


# Status_of_existing_checking_account가 가장 강력한 예측 능력(의심 될 정도로)를 가졌고, 다음이 Credit_history이다. Telephone은 가장 낮은 예측력을 가지고 있다.

# 사용할 컬럼은 본인이 선택해야한다. 나은 예측을 위해, 모델을 생성하고, 테스트하고 다시 돌아와서 컬럼을 생성하고를 반복한다.

# 우선 상위 15개만을 선택하자

# In[93]:


# 분석하려는 feature의 갯수와 예측력의 trade-off를 조정하였습니다.
iv_cols = [key for key, iv in candidates if iv >= 0.044]
display(len(iv_cols))
display(sorted(iv_cols))


# In[94]:


# 상수로 설정합니다.
IV_COL = iv_cols[:]


# Information Value에 의하여 상위 15개의 feature를 선정하였습니다. 전체 feature를 사용하여 로지스틱 회귀를 구했을 때의 accuracy와 신용불량자에 대한 precision의 감소폭이 적으면서 (전체 데이터셋을 가장 잘 대변해 줄 수 있으면서) feature의 수를 줄이는 과정을 몇 번 테스트 하면서 얻은 값입니다.
# 
# 먼저, IV에서 선택된 컬럼들로 이루어진 모델에서부터 시작합니다.

# In[111]:


# Category column들은 one hot vectore로 변환합니다.
cate_features = {}
num_fetures = []
for col in IV_COL:
    if dataset[col].dtype == 'object':
        cate_features[col] = pd.get_dummies(dataset[col], prefix=col)
    else:
        num_fetures.append(col)


# In[112]:


cate_features.keys()


# In[113]:


cate_features['Status_of_existing_checking_account'].head()


# 더미 변수를 생성하게 되면, 자유도는 더미 변수가 가질 수 있는 unique 값이 $n$ 이라고 했을때 $n-1$이 되므로, 각 변수에서 더미변수 축을 한개씩 제거해준다. 제거하지 않으면 `NaN`값이 발생합니다.
# 
# 후진 제거법을 통해 무의미한 변수를 하나씩 제거해 나가기 위해 별도의 리스트를 하나 생성합니다. 이 리스트는 실행이 반복될 때마다 가장 무의미한 변수(p-value가 큰 변수)와 다중공선성을 갖는 변수를 제거 한다.

# In[114]:


removed_features = []


# In[115]:


for col, dummies in cate_features.items():
    dropped_col = dummies.columns[-1]
    removed_features.append(dropped_col)
    cate_features[col] = dummies.drop(dropped_col, axis=1)


# In[116]:


removed_features


# In[117]:


cate_features['Status_of_existing_checking_account'].head()


# In[126]:


final_dataset = dataset[num_fetures]


# In[127]:


for col, df in cate_features.items():
    final_dataset = pd.concat([final_dataset, df], axis=1)


# In[128]:


final_dataset.shape


# In[129]:


final_dataset.head()


# ## 2. Training set, test set을 나눕니다. 

# In[110]:


from sklearn.model_selection import train_test_split


# In[132]:


train_x, test_x, train_y, test_y = \
    train_test_split( \
        final_dataset, dataset['Target'], test_size=0.2, random_state=42)


# In[133]:


print('train_set', train_x.shape)
print('test_set', test_x.shape)


# ## 3. Statsmodels를 이용하여 로지스틱 회귀를 돌립니다.

# In[134]:


import statsmodels.api as sm


# In[153]:


logistic_model = sm.Logit(
    train_y,
    sm.add_constant(train_x)
).fit()


# ### 성능 측정

# In[154]:


from sklearn.metrics import auc
from sklearn.metrics import accuracy_score
from sklearn.metrics import classification_report


# In[155]:


train_pred = pd.DataFrame({
    'probs': logistic_model.predict(sm.add_constant(train_x)),
    'class': train_y
})

train_pred['y_pred'] = 0
train_pred.loc[train_pred['probs'] > 0.5, 'y_pred'] = 1

# Test prediction
test_pred = pd.DataFrame({
    'probs': logistic_model.predict(sm.add_constant(test_x)),
    'class': test_y
})

test_pred['y_pred'] = 0
test_pred.loc[test_pred['probs'] > 0.5, 'y_pred'] = 1


# In[156]:


print('\nTraining Confusion matrix:')
display(pd.crosstab(train_pred['y_pred'], train_pred['class'], rownames=['Predict'], colnames=['Actual'], margins=True))
print('\nTraining Accuracy: ', round(accuracy_score(train_pred['class'], train_pred['y_pred']), 4))
print('\nTraining classification report:\n', classification_report(train_pred['class'], train_pred['y_pred'], digits=4))


# In[157]:


print('Test Confusion matrix:')
display(pd.crosstab(test_pred['y_pred'], test_pred['class'], rownames=['Predict'], colnames=['Actual'], margins=True))
print('\nTest Accuracy: ', round(accuracy_score(test_pred['class'], test_pred['y_pred']), 4))
print('\nTest classification report:\n', classification_report(test_pred['class'], test_pred['y_pred'], digits=4))


# ## 2. 후진제거법
#  - 가장 무의미한 변수 (p-value)
#  - 다중공선성을 갖는 변수 (Variance Inflation Factor)

# ### (A) P-value
# P-value가 가장 높은 변수를 제거합니다. 여기서 P-value는 각 독립변수의 유의성 검정(t-test)에서 얻어진 값입니다.

# In[179]:


unnecesarries = []


# In[180]:


logistic_model.summary2()


# 가장 무의미한 변수는 P값이 0.9580인 Personal_status_and_sex_A92이다.

# 체크하고 변수를 1개씩 제거하고를 반복해야한다.

# In[181]:


# Help 함수
def get_unrelated_cols(model, pvalue):
    cols = model.pvalues[model.pvalues >= pvalue].keys().tolist()
    print(len(cols))
    print(cols)
    
    return cols


# In[182]:


unrelated_cols = get_unrelated_cols(logistic_model, 0.8187)


# Personal_status_and_sex_A92을 삭제한다.

# In[174]:


unnecesarries.append('Personal_status_and_sex_A92')


# ### (B) 다중공선성
# 테스트와 상관없다고 판단되는 변수들을 제거한후에 다중공선성을 체크합니다. **다중공선성(multicollinearity)** 이란 독립변수들간의 선형관계가 존재하는 것을 나타냅니다. 독립변수 전체의 set에서 하나의 독립변수를 골라 종속변수로 보고, 그 독립변수를 제외한 set을 이용하여 선형회귀분석을 하여 구할 수 있습니다.
# 
# 사실 변수간의 선형관계 알아보기 위해서는 산점도 or 상관계수를 사용해도 되지만 **VIF(분산 팽창 지수, Variance Inflation Factor)** 를 이용하여 구한다.

# In[183]:


def get_max_vif(df, removal_cols):
    vifs = []
    cnames = df.drop(removal_cols, axis=1).columns.tolist()
    for i in range(len(cnames)):
        xvar = cnames[:]
        yvar = xvar.pop(i)
        model = sm.OLS(
            df.drop(removal_cols, axis=1)[yvar], 
            sm.add_constant(df.drop(removal_cols, axis=1)[xvar]))
        res = model.fit()
        vif = 1 / (1 - res.rsquared)
        vifs.append((yvar, round(vif, 3)))
    vifs = sorted(vifs, key=operator.itemgetter(1), reverse=True)
    return vifs


# In[184]:


vifs = get_max_vif(train_x, unnecesarries)


# In[185]:


vifs


# In[177]:


unnecesarries.append(vifs[0][0])


# In[178]:


unnecesarries


# 무의미한 변수(가장 높은 p-value)는 Personal_status_and_sex_A92였고 pvalue가 0.9580, VIF가 3.302이다.
# 
# 다중공선성이 가장 높았던 변수는 Housing_A152였고 pvalue가 0.4770, VIF가 5.846이다.
# 
# **무의미한 변수를 먼저 삭제해야하며, 지우기 전에, 두 변수의 pvalue와 VIF를 각각 비교해보는것이 좋다**.
# 
# 확실히 Personal_status_and_sex_A92는 삭제해야 한다.

# 모델을 삭제하기 앞서 모델 채택 기준을 살펴봐야한다.

# ## 3. 모델 채택 기준
# ### (A) C통계량 (concordance statistics)
#  - 로지스틱 회귀 모델에서 이항 결과의 적합도(goodness of fit)에 관한 품질을 측정하는 척도
#  - ROC curve와 동일하다.

# $$C=Pr[\pi(B|x_{i}) > \pi(B|x_{j}) | Y_{i} = 1, Y_{j} = 0]$$
# $$C 통계량(C 인덱스) = 0.5 + (\frac{일치쌍 비율 - 불일치쌍 비율}{2})$$

# In[189]:


def get_c_stat(iv_pred):
    noraml_test_df = iv_pred[iv_pred['class'] == 0][['class', 'probs']]
    spammer_test_df = iv_pred[iv_pred['class'] == 1][['class', 'probs']]

    noraml_test_df['key'] = 0
    spammer_test_df['key'] = 0

    cross_join_df = noraml_test_df.merge(spammer_test_df, how='outer', on='key').drop('key', axis=1)

    cross_join_df['concordance'] = cross_join_df['probs_x'] < cross_join_df['probs_y']
    cross_join_df['in_concordance'] = cross_join_df['probs_x'] > cross_join_df['probs_y']
    cross_join_df['tie'] = cross_join_df['probs_x'] == cross_join_df['probs_y']

    results = cross_join_df.agg({'concordance': np.sum, 'in_concordance': np.sum, 'tie': np.sum}) / len(cross_join_df)
    c_stat = 0.5 + (results['concordance'] - results['in_concordance']) / 2
    
    return c_stat


# ### (B) AIC, Likelihood
# 
# - $AIC=-2 \times ln(L) + 2 \times k$
# 
# Akaike information criterion는 주어진 데이터 집합에 관해 통계 모델의 상대적 품질을 측정한다. 이 척도는 편향과 분산의 트레이드 오프이다. 더 작은 AIC를 선호한다.
# 
# - $L = \Pi_{i=1}^{n}p(x_i)^{y_i}(1-p(x_i))^{1-y_i}$
# 
# log likelihood는 큰 값을 선호한다.

# In[190]:


def get_aic_value(model):
    return -2 * model.llf + 2 * (len(model.params) - 1)


# > 무의미한 변수를 지우기 전

# In[211]:


c_stat, aic = get_c_stat(train_pred), get_aic_value(logistic_model)
print('c_stat:', c_stat)
print('aic:', aic)
print('loglikehood:', logistic_model.llf)


# In[193]:


c_stat, aic = get_c_stat(test_pred), get_aic_value(logistic_model)
print('c_stat:', c_stat)
print('aic:', aic)
print('loglikehood:', logistic_model.llf)


# > 지우고 난 후

# In[202]:


unnecesarries.extend(['Personal_status_and_sex_A92', 'Housing_A152'])


# In[208]:


drop_logistic_model = sm.Logit(
    train_y,
    sm.add_constant(train_x.drop(unnecesarries, axis=1))
).fit()


# In[213]:


train_pred = pd.DataFrame({
    'probs': drop_logistic_model.predict(sm.add_constant(train_x.drop(unnecesarries, axis=1))),
    'class': train_y
})

train_pred['y_pred'] = 0
train_pred.loc[train_pred['probs'] > 0.5, 'y_pred'] = 1

# Test prediction
test_pred = pd.DataFrame({
    'probs': drop_logistic_model.predict(sm.add_constant(test_x.drop(unnecesarries, axis=1))),
    'class': test_y
})

test_pred['y_pred'] = 0
test_pred.loc[test_pred['probs'] > 0.5, 'y_pred'] = 1


# In[214]:


c_stat, aic = get_c_stat(train_pred), get_aic_value(drop_logistic_model)
print('c_stat:', c_stat)
print('aic:', aic)
print('loglikehood:', drop_logistic_model.llf)


# In[215]:


c_stat, aic = get_c_stat(test_pred), get_aic_value(drop_logistic_model)
print('c_stat:', c_stat)
print('aic:', aic)
print('loglikehood:', drop_logistic_model.llf)


# In[206]:


drop_logistic_model.summary2()


# ## 4. 무의미한 변수와 다중공선성이 있는 변수가 없을때까지 반복한다. (기준은 C통계량, AIC, Likelihood)

# 원래는 이 과정을 각각 변수를 확인하면서, 변수를 하나씩 제거해야 한다. 하지만 시간 관계상..

# #### Step 1. Pvalue

# In[218]:


# Pvalue를 확인하고, 상한선을 구한다.
unrelated_cols = get_unrelated_cols(drop_logistic_model, 0.5)


# In[219]:


unnecesarries.extend(unrelated_cols)


# #### Step 2. VIF

# In[221]:


vifs = get_max_vif(train_x, unnecesarries)


# In[222]:


vifs[:10]


# In[226]:


unnecesarries.extend([elem[0] for elem in vifs[:10]])


# In[228]:


drop_logistic_model = sm.Logit(
    train_y,
    sm.add_constant(train_x.drop(unnecesarries, axis=1))
).fit()


# In[229]:


train_pred = pd.DataFrame({
    'probs': drop_logistic_model.predict(sm.add_constant(train_x.drop(unnecesarries, axis=1))),
    'class': train_y
})

train_pred['y_pred'] = 0
train_pred.loc[train_pred['probs'] > 0.5, 'y_pred'] = 1

# Test prediction
test_pred = pd.DataFrame({
    'probs': drop_logistic_model.predict(sm.add_constant(test_x.drop(unnecesarries, axis=1))),
    'class': test_y
})

test_pred['y_pred'] = 0
test_pred.loc[test_pred['probs'] > 0.5, 'y_pred'] = 1


# In[230]:


c_stat, aic = get_c_stat(train_pred), get_aic_value(drop_logistic_model)
print('c_stat:', c_stat)
print('aic:', aic)
print('loglikehood:', drop_logistic_model.llf)


# In[231]:


c_stat, aic = get_c_stat(test_pred), get_aic_value(drop_logistic_model)
print('c_stat:', c_stat)
print('aic:', aic)
print('loglikehood:', drop_logistic_model.llf)


# In[232]:


drop_logistic_model.summary2()


# In[235]:


print('Train Confusion matrix:')
display(pd.crosstab(train_pred['y_pred'], train_pred['class'], rownames=['Predict'], colnames=['Actual'], margins=True))


# In[233]:


print('Test Confusion matrix:')
display(pd.crosstab(test_pred['y_pred'], test_pred['class'], rownames=['Predict'], colnames=['Actual'], margins=True))


# ## 5. ROC curve, AUC
#  - ROC = FPR vs TPR(Recall) 의 비율
#  - AUC = ROC 곡선 아래의 면적

# In[237]:


import matplotlib.pyplot as plt
from sklearn import metrics
from sklearn.metrics import auc

plt.rcParams["figure.figsize"] = (7,7)


# In[238]:


fpr, tpr, thresholds = metrics.roc_curve(train_pred['class'], train_pred['probs'], pos_label=1)


# In[240]:


roc_auc = auc(fpr, tpr)


# In[241]:


plt.figure()
lw = 2
plt.plot(fpr, tpr, color='darkorange', lw=lw, label='ROC curve (area = %0.4f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel('False Postive Ratio')
plt.ylabel('True Postive Ratio(Recall)')
plt.title('ROC curve')
plt.legend(loc='lower right')
plt.show()


# ## 6. 결정경계

# In[253]:


for i in np.arange(0, 1, 0.1):
    train_pred['y_pred'] = 0
    train_pred.loc[train_pred['probs'] > i, 'y_pred'] = 1
    acc = round(accuracy_score(train_pred['class'], train_pred['y_pred']), 4)
    print("Threshold", round(i,1), "Train accuracy:",acc)


# In[254]:


for i in np.arange(0, 1, 0.1):
    test_pred['y_pred'] = 0
    test_pred.loc[test_pred['probs'] > i, 'y_pred'] = 1
    acc = round(accuracy_score(test_pred['class'], test_pred['y_pred']), 4)
    print("Threshold", round(i,1), "Test accuracy:",acc)


# In[261]:


# Test prediction
test_pred = pd.DataFrame({
    'probs': drop_logistic_model.predict(sm.add_constant(test_x.drop(unnecesarries, axis=1))),
    'class': test_y
})

test_pred['y_pred'] = 0
test_pred.loc[test_pred['probs'] > 0.4, 'y_pred'] = 1


# In[262]:


print('Test Confusion matrix:')
display(pd.crosstab(test_pred['y_pred'], test_pred['class'], rownames=['Predict'], colnames=['Actual'], margins=True))