[윤형준] - Section 2 Project - 위생검사 통과 여부 예측하기

이번 프로젝트는 머신러닝을 통해 위생검사 통과 여부를 예측하는 모델을 만들어 보겠습니다.

이번 프로젝트에서 중점적으로 다룰 내용은 다음과 같습니다.

**1. 특성별 다양한 인코더를 알맞게 사용하기**

**2. 모델 검증 방법 - Holdout Validation**

**3. 평가 지표 - F1 Score / ROC AUC**

**4. 사용 모델 -> DecisionTreeClassifier() / RandomForestClassifier()**

**5. Hyperparameter Tuning - RandomizedSearchCV**

**사용된 데이터:** 시카고에서 상업목적으로 비지니스가 이루어지고 있는 업체들을 대상으로 10년간 진행된 위생검사 통과 여부를 기록한 데이터를 기반으로 프로젝트를 진행하였습니다.

우리가 예측하고자 하는 **Target** 값은 **Inspection fail** 입니다. [0, 1] 로 이루어진 이진 데이터 이며, 0으로 표시된것은 위생검사를 통과한 업체이고, 1로 표시된것은 위생검사를 통과하지 못한 업체입니다. 우리가 만들어야 할 모델은 새로운 데이터를 모델에 넣어서 **Target** 값을 알맞게 예측하는 것입니다. (**A.K.A. 분류문제**)

Python Code - 데이터 불러오기

```py

import pandas as pd
import numpy as np

train_url = 'https://ds-lecture-data.s3.ap-northeast-2.amazonaws.com/food_inspection_sc23x/food_ins_train.csv'
test_url = 'https://ds-lecture-data.s3.ap-northeast-2.amazonaws.com/food_inspection_sc23x/food_ins_test.csv'

# train, test 데이터셋을 불러옵니다
train = pd.read_csv(train_url)
test = pd.read_csv(test_url)

```


**Exploratory Data Analysis (EDA)**

데이터는 언제나 깔끔하게 저장되어 있지 않습니다. 우리가 모델을 만들어 데이터를 학습 시키기 위해서는 데이터가 모두 숫자형데이터이여야 합니다. 그래서 숫자형이 아닌 데이터는 **feature engineering**을 통해 숫자형으로 변형 시켜야 합니다. 그렇다면 어떤 특성의 데이터가 숫자형인지, 그리고 어떤 특성에 데이터가 빈값이 있는지 등을 EDA를 통해 확인해야 합니다.

**우선 타겟값의 분포를 보겠습니다.**

0에 분포가 1보다 훨씬 많은 편향성을 띄고 있습니다.

Python Code - 타겟값 분포확인하기

```py
plt.hist(train["Inspection Fail"])
plt.title("Target Values")
plt.ylabel("Frequency")
plt.show()
```


**이번엔 각 특성들의 Data 종류와 null값, 그리고 unique값을 확인하겠습니다.**

Python Code - EDA Table 만들기

```py
def make_EDA_table(df):
columns = df.columns
d_types = df.dtypes
d_nulls = df.isnull().sum()
d_uniques = pd.Series([len(df[column].unique()) for column in columns], index=df.columns)
d_frame = pd.DataFrame([d_types, d_nulls, d_uniques]).T
d_frame = d_frame.rename(columns={0:"data type", 1:"#_of_null", 2:"#_of_unique"})

return d_frame

make_EDA_table(train)
```


**Feature Engineering**

앞서 말했듯이 **Feature Engineering**은 머신러닝에서 굉장히 중요한 부분입니다. 데이터를 모두 숫자형으로 변환 시켜야 할 뿐 아니라, 각 특성들의 성질을 제대로 이해하고 알맞게 engineering 해야합니다. 그러다면, 특성들을 engineering하기전에 어떤 특성에 어떤 encoder를 사용할지 먼저 생각해야합니다.

올디널인코딩 참고자료: https://pycaret.org/ordinal-encoding/
원핫인코딩 참고자료: https://hackernoon.com/what-is-one-hot-encoding-why-and-when-do-you-have-to-use-it-e3c6186d008f
타겟인코딩 참고자료: https://medium.com/analytics-vidhya/target-encoding-vs-one-hot-encoding-with-simple-examples-276a7e7b3e64

**DBA Name | AKA Name | License #**

위의 EDA Table을 보면, 이 세가지 특성들은 unique한 값들이 너무 많습니다. 따라서 OnehotEncoding을 사용하면 High Cardinality 문제가 생깁니다. 따라서 전 해당 특성들에 Target Encoding을 사용할 예정입니다.

**Facility Type**

Facility Type이라는 특성도 unique한 값이 많은 편입니다. 하지만 해당 특성을 좀더 세부적으로 보겠습니다. 아래 테이블을 보면 Resturant의 분포가 제일 큼니다. 사실상 School 밑으로는 상대적으로 분포가 너무 적습니다. 따라서 feature egineering시 School보다 분포가 적은 데이터는 모두 묶어서 etc로 분류할것입니다. 이렇게 engineering을 하면 unique한 값이 총 4개가 됩니다. 4개의 unique한 값을 대상으로 OneHotEncoding을 진행할것입니다.

```py
train["Facility Type"].value_counts(normalize=True).head(20)
```


**Inspection Type**

Inspection Type도 위와 같은 논리로 상위 6개의 값들만 남겨두고 나머지는 모두 etc로 분류하겠습니다. 그리고 etc를 포함한 총 7개의 unique한 값을 OneHotEncoding할 예정입니다.

```py
train["Inspection Type"].value_counts(normalize=True).head(20)
```


**Inspection Date**
아래 그림을 보면, 연초 연말에 inspection이 가장 많았던 것으로 보입니다. 그리고 년도 그림을 보면 최근에 inspection이 가장 많았던 것으로 보입니다. 하지만 나중에 feature engineering을 할 시에 month 데이터만 쓰고, year 데이터는 쓰지 않겠습니다.

물론 현재 셋업에서는 가지고 있는 데이터에서 어떤 한 부분을 test 셋으로 사용할 예정이지만, 실제로는 새로운, 즉 앞으로 2021년 데이터를 추측할텐데, 우리는 미래에 대한 데이터는 가지고있지 않지 않습니다. 따라서 과거 년도 데이터를 학습 시키더라도 미래 년도 데이터를 예측하기 어렵기때문에 년도 데이터는 쓰지 않겠습니다.

month를 사용하는 이유는 한 해에 특정 달에 inspection pass 와 fail에 관련이 있을수 있다고 판단했기 때문입니다. month는 ordinal data 로 보지않고 label data로 볼것이며, 따라서 target encoding을 할 예정입니다.

Python Code - Inspection Date Histogram

```py
train_inspect_year = pd.to_datetime(train["Inspection Date"]).dt.year
test_inspect_year = pd.to_datetime(test["Inspection Date"]).dt.year

train_inspect_month = pd.to_datetime(train["Inspection Date"]).dt.month
test_inspect_month = pd.to_datetime(test["Inspection Date"]).dt.month

import matplotlib.pyplot as plt

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(15,10))

ax[0,0].hist(train_inspect_year)
ax[0,0].title.set_text("Year - Train Set")

ax[0,1].hist(test_inspect_year)
ax[0,1].title.set_text("Year - Test Set")

ax[1,0].hist(train_inspect_month)
ax[1,0].title.set_text("Month - Train Set")

ax[1,1].hist(test_inspect_month)
ax[1,1].title.set_text("Month - Test Set")

plt.show()
```


**Violations**

```py
train["Violations"][1]
```


Violation 특성은 해당 업체가 위반한 조항들을 나열하고 있습니다. 이 특성은 데이터를 string 형으로 저장하고 있습니다. 이 특성에서 의미있는 데이터를 추출하기 위해선 해당 업체가 몇가지의 조항을 위반하고 있는지를 세어서 학습시킬수 있습니다. 이렇게 하는 이유는 위반한 조항의 갯수가 많을 수록 Inspection Fail에 확률이 높다고 생각하기 때문입니다.

> str()을 사용해서 "|" 으로 쪼갠후 위반한 조항 갯수를 셀 것입니다.

> 각 업체마다 위반 조항 갯수를 세어서 [num of violations] 라는 새로운 컬럼에 저장할것입니다.

> [Violations]라는 특성은 drop 할 것입니다.

> [num of violations] 라는 특성은 순서형 데이터입니다. 특성 그대로를 인코딩하지 않고 학습시킬 예정입니다.

**Risk**

```py
train["Risk"].unique()
```


Risk가 높을수록 inspection을 fail할 확률이 높습니다. 따라서 Risk 특성을 Ordinal Encoding 할 예정입니다.

**아래 code를 실행하여 Feature Engineering을 해보세요!**

Python Code - Feature Engineering Code

```py

# Drop Unnecessary Columns

drop_list = ["Inspection ID", "Address", "City", "State", "Zip", "Latitude", "Longitude", "Location"]

def drop_features(train, test):

train_copy = train.copy()
test_copy = test.copy()

train_copy.drop(columns = drop_list, inplace=True)
test_copy.drop(columns = drop_list, inplace=True)

return train_copy, test_copy

train_dp, test_dp = drop_features(train, test)

# Facility Type | Inspection Type

# 아래 function은 해당 특성에 상위 N개의 category만 남기고 모두 etc로 분류합니다.
# 이 function은 Facility Type 과 Inspection Type 특성에 사용될 것입니다.

def extract_top_N_cat(train, test, feature, N):

train_top_N_cat = train[feature].value_counts(normalize=True)[:N].index.values
test_top_N_cat = test[feature].value_counts(normalize=True)[:N].index.values

def classify_cat_train(element):
if element not in train_top_N_cat:
return "etc"
else:
return element

def classify_cat_test(element):
if element not in test_top_N_cat:
return "etc"
else:
return element

train_copy = train.copy()
test_copy = test.copy()

train_copy[feature] = train_copy[feature].apply(classify_cat_train)
test_copy[feature] = test_copy[feature].apply(classify_cat_test)

return train_copy, test_copy

# Inspection Date

def get_months (train, test):

train_copy = train.copy()
test_copy = test.copy()

train_copy["Inspection Month"] = pd.to_datetime(train_copy["Inspection Date"]).dt.month
test_copy["Inspection Month"] = pd.to_datetime(test_copy["Inspection Date"]).dt.month

train_copy.drop(columns=["Inspection Date"], inplace=True)
test_copy.drop(columns=["Inspection Date"], inplace=True)

return train_copy, test_copy

# Violation

def get_vio_nums(train, test):

def count_num_violations(element):
num_violations = len(str(element).split("|"))

return num_violations

train_copy = train.copy()
test_copy = test.copy()

train_copy["num_violations"] = train_copy["Violations"].apply(count_num_violations)
test_copy["num_violations"] = test_copy["Violations"].apply(count_num_violations)

train_copy.drop(columns=["Violations"], inplace=True)
test_copy.drop(columns=["Violations"], inplace=True)

return train_copy, test_copy

# Aggregating feature engineering functions

def engineering_aggregation(train, test):

train_copy = train.copy()
test_copy = test.copy()

# Facility Type -> extract top 3
train_1, test_1 = extract_top_N_cat(train_copy, test_copy, feature="Facility Type", N=3)

# Inspection Type -> extract top 6
train_2, test_2 = extract_top_N_cat(train_1, test_1, feature="Inspection Type", N=6)

# Inspection Date -> get only months
train_3, test_3 = get_months(train_2, test_2)

# Violations -> get number of violations
train_4, test_4 = get_vio_nums(train_3, test_3)

return train_4, test_4

eng_train, eng_test = engineering_aggregation(train_dp, test_dp)

```

**Holdout Validation Method**

우선 인코딩하기 전에, 나중에 만들 모델을 평가하기위해서 Train Set과 Validation Set을 나누겠습니다. 또한 우리가 학습시킬 특성들과 타켓값도 나누겠습니다. 이번 프로젝트에서 저는 **Holdout Validation Method**를 사용할 예정이기 때문에 Train Set과 Validation Set을 나눕니다.

```py
target = "Inspection Fail"
features = eng_train.columns.drop(target)

train_set = eng_train.sample(frac=0.7, random_state=1)
val_set = eng_train.drop(index=train_set.index)

X_train = train_set[features]
y_train = train_set[target]

X_val = val_set[features]
y_val = val_set[target]

X_test = eng_test[features]
y_test = eng_test[target]
```

**파이프 라인을 활용하여 특성별로 인코딩하기**

```py
from sklearn.pipeline import make_pipeline

from category_encoders import TargetEncoder
from category_encoders import OneHotEncoder
from category_encoders import OrdinalEncoder

target_encoding_cols = ["DBA Name", "AKA Name", "License #", "Inspection Month"]
onehot_encoding_cols = ["Facility Type", "Inspection Type"]

ordinal_encoding_sub_map = {'Risk 1 (High)':1, 'Risk 2 (Medium)':2, 'Risk 3 (Low)':3}
ordinal_encoding_map = [{"col":"Risk", "mapping":ordinal_encoding_sub_map}]

enc_pipeline = make_pipeline(

# 아래 encoder들은 default로 null 값들을 mean 값으로 계산해서 집어넣습니다.
# 만약에 나중에 RandomizedSearchCV를 통해 hyperparameter 튜닝을 원하다면 default setting을 변경시켜야 합니다.
TargetEncoder(cols=target_encoding_cols, smoothing=1000),
OneHotEncoder(cols=onehot_encoding_cols, use_cat_names=True),
OrdinalEncoder(mapping=ordinal_encoding_map)
)

processed_X_train = enc_pipeline.fit_transform(X_train, y_train)
processed_X_val = enc_pipeline.transform(X_val)
```
TargetEncoder API: https://contrib.scikit-learn.org/category_encoders/targetencoder.html
OneHotEncoder API: https://contrib.scikit-learn.org/category_encoders/onehot.html
OrdinalEncoder API: http://contrib.scikit-learn.org/category_encoders/ordinal.html

**인코딩 전후 비교**

인코딩 하기전 X_train 데이터


인코딩 후 X_train 데이터


모든 특성이 숫자형으로 변형되었습니다. 이제 모델을 만들어 학습시킬 차례입니다!

**DecisionTreeClassifier()**

```py
from sklearn.tree import DecisionTreeClassifier

clf_dt = DecisionTreeClassifier(random_state=1)
clf_dt.fit(processed_X_train, y_train);

y_train_pred_dt = clf_dt.predict(processed_X_train)
y_val_pred_dt = clf_dt.predict(processed_X_val)
```
우리는 hyperparameter를 튜닝하지 않은 이 DecisionTreeClassifier()를 baseline 모델로 사용하겠습니다.

# 모델이 분류를 잘 한다는건 무슨 뜻일까요?

아래 그림을 보겠습니다. 아래 그림은 추상의 모델로 데이터를 분류하는 것을 시뮬레이션 해보았습니다. 오렌지색 분포는 타겟값이 positive로 분류되는 데이터이고, 파란색 분포는 타겟값이 negative로 분류되는 데이터 입니다. 하지만 여기서 가상의 모델은 Threshold = 0.5 를 기준으로 0.5보다 확률이 큰 데이터는 positive로 분류하고 0.5보다 확률이 작은 데이터는 negative로 분류하고 있습니다.

이 모델은 분류를 잘 못하는 모델입니다. 실제 Positive가 일정부분 모델에 의해 Negative로 분류되고 있고, 실제 Negative가 일정부분 모델에 의해 Possitive로 분류 되고 있기 때문입니다. 즉 Type I error (False Positive) 과 Type II error (False Negative) 가 그림상에도 많이 보이기 때문입니다.

더 좋은 모델은 Inspection Fail distribution과 Inspection Pass distribution을 더 멀리 떨어뜨립니다. (즉 분류를 더 잘한다는 뜻!) 두 분포가 더 떨어져있으면 Type I error와 Type II error가 더 줄어듭니다.

Python Code

```py
plt.figure(figsize=(8,6))
np.random.seed(1)
neg_data = np.random.normal(0.25,0.19,1000)
pos_data = np.random.normal(0.75,0.19,1000)
plt.hist(neg_data, bins=30, alpha=0.5, range=(0,1), label="negative")
plt.hist(pos_data, bins=30, alpha=0.5, range=(0,1), label="positive")
plt.axvline(x=0.5, color="red", linestyle="--", label="threshold")
plt.xlabel("Probabilities")
plt.ylabel("Frequency")
plt.title("Bad Model")
plt.legend(loc="upper right")
plt.show()
```


다음 그림은 분류를 잘하는 모델입니다. 앞서 말했듯이, 이 모델은 Type I error와 Type II error가 많지 않습니다.

Python Code

```py
plt.figure(figsize=(8,6))
np.random.seed(1)
neg_data = np.random.normal(0.25,0.09,1000)
pos_data = np.random.normal(0.75,0.09,1000)
plt.hist(neg_data, bins=30, alpha=0.5, range=(0,1), label="negative")
plt.hist(pos_data, bins=30, alpha=0.5, range=(0,1), label="positive")
plt.axvline(x=0.5, color="red", linestyle="--", label="threshold")
plt.xlabel("Probabilities")
plt.ylabel("Frequency")
plt.title("Good Model")
plt.legend(loc="upper right")
plt.show()

```


## Biased Target Values (불균등 타겟값)

위의 두 예시 그림은 우리의 target 데이터 (positive 와 negative)가 균등하게 있다고 가정한 것입니다. 하지만 현재 우리가 다루고 있는 데이터처럼 target 데이터가 불균등한 경우가 있습니다. 불균등한 타겟값이 있을때 모델이 성능이 떨어진다면 (즉 분류를 잘 못한다면) 아래와 같은 그림이 나타납니다.

이렇게 모델이 분류를 잘 하지 못한다면 갯수가 적은 타겟값이 갯수가 큰 타겟값 분포에 가려질 수가있습니다.

Python Code

```py
plt.figure(figsize=(8,6))
np.random.seed(1)
neg_data = np.random.normal(0.25,0.19,1700)
pos_data = np.random.normal(0.5,0.19,300)
plt.hist(neg_data, bins=30, alpha=0.5, range=(0,1), label="negative")
plt.hist(pos_data, bins=30, alpha=0.5, range=(0,1), label="positive")
plt.axvline(x=0.58, color="red", linestyle="--", label="threshold")
plt.xlabel("Probabilities")
plt.ylabel("Frequency")
plt.title("Bad Model")
plt.legend(loc="upper right")
plt.show()

```


하지만 불균등한 타겟값을 가지고 있는 데이터이더라도 모델이 분류를 잘 한다면 다음 그림과 같은 좋은 결과를 얻을수 있습니다.

Python Code

```py
plt.figure(figsize=(8,6))
np.random.seed(1)
neg_data = np.random.normal(0.25,0.19,1700)
pos_data = np.random.normal(0.85,0.09,300)
plt.hist(neg_data, bins=30, alpha=0.5, range=(0,1), label="negative")
plt.hist(pos_data, bins=30, alpha=0.5, range=(0,1), label="positive")
plt.axvline(x=0.7, color="red", linestyle="--", label="threshold")
plt.xlabel("Probabilities")
plt.ylabel("Frequency")
plt.legend(loc="upper right")
plt.show()

```


## 우리의 Baseline Model - DecisionTreeClassifier()는 분류를 잘하고 있을까요?

> 실제 데이터가 어떻게 분포되어 있는지 보고

> 현 모델이 분류를 얼마나 잘했는지 봅시다.

Python Code

```py
y_val_pred_proba_dt = clf_dt.predict_proba(processed_X_val)[:,1]
prob_n_actual_comps = pd.DataFrame({"pred_proba":y_val_pred_proba_dt, "target_values":y_val})

actually_fail = prob_n_actual_comps[prob_n_actual_comps["target_values"]==1]
actually_pass = prob_n_actual_comps[prob_n_actual_comps["target_values"]==0]

plt.figure(figsize=(8,6))
plt.title("DecisionTreeClassifier() Model")
plt.yscale('log')
plt.hist(actually_fail["pred_proba"], bins=20, alpha=0.5, label="Inspection Fail")
plt.hist(actually_pass["pred_proba"], bins=20, alpha=0.5, label="Inspection Pass")
plt.axvline(x=0.5, color="red", linestyle="--", label="threshold")
plt.xlabel("Probabilities")
plt.ylabel("Log(Frequency)")
plt.legend()

plt.show()

```


현 모델은 Threshold를 기준으로 확률이 0.5보다 큰것은 Fail로 예측하고 0.5보다 작은것은 Pass로 예측하고 있습니다. 한눈에 볼 수 있듯이 현 모델은 매우 안좋은 모델입니다. Inspection Fail과 Inspection Pass를 제대로 분류하고 있지 못합니다.

실제 Inspection Fail이 많은부분 모델에 의해 pass로 분류되고 있고, 실제 Inspection Pass가 많은부분 모델에 의해 fail로 분류 되고 있습니다. 즉 Type I error (False Positive) 과 Type II error (False Negative) 가 그림상에 많이 보입니다.

지금부터 우리의 목표는 이 baseline model보다 성능이 더 좋은 모델을 만드는 것입니다.

(side note: 사실 우리 데이터의 타겟값은 불균등 하기 때문에 앞서 본것처럼 하나의 분포가 매우 커야되고 다른 하나의 분포는 매우 작아야 합니다. 하지만 현 그림에서는 그렇게 보이지 않고있습니다. 그 이유는 y-axis를 log로 rescaling을 했기 때문입니다. Outlier 값이 존재하여 log로 rescaling 했습니다.)

### Confusion Matrix

위의 그림에 대한 결과를 confusion matrix로 보겠습니다. 위의 그림들에선 시각화로 모델이 분류를 어떻게 하는지 전체적인 청사진의 대한 이해를 돕는 반면 정확한 숫자들을 알기는 어렵습니다. 하지만 Confusion Matrix는 모델이 어떻게 분류했는지 숫자로 나타나기 때문에 몇개의 데이터를 올바르게 예측했는지를 알수있습니다.

Python Code

```py
from sklearn.metrics import plot_confusion_matrix
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
pcm = plot_confusion_matrix(clf_dt, processed_X_val, y_val, cmap=plt.cm.Blues,
ax=ax, labels=[1,0], display_labels = ["Fail(1)","Pass(0)"]);
plt.title(f'Confusion matrix, n = {len(y_val)}', fontsize=15)
plt.show()

```


## 우리가 만든 모델을 평가하는 지표

classification_report를 통해서 우리가 만든 모델로 validation set을 예측했을때 "얼마나 잘 예측했나?"를 평가하는 여러 지표를 보겠습니다.

```py
from sklearn.metrics import classification_report
print(classification_report(y_val, y_val_pred_dt))
```


앞서 말했듯이 시각적으로는 negative의 분포와 positive의 분포가 더 떨어져있을 수록 좋은 모델이라고 말했습니다. 지표로 이야기 하자면 더 좋은 모델은 f1(negative)와 f1(positive)가 크면 클수록 더 좋은 모델이라고 말할수 있습니다. 현재 우리의 Baseline model의 f1(pos)는 0.34이고 f1(neg)는 0.86입니다. 이것들보다 더 높은 f1 score들을 낼수있는 모델을 만드는것이 우리의 목표입니다.

#### ROC AUC Score (ROC Area Under the Curve)

우리는 f1(pos)와 f1(neg)를 함께 보아야 합니다. 물론 더 좋은 모델은 f1(pos)와 f1(neg)가 둘다 baseline model의 f1 score들보다 좋아야합니다. 하지만 어떠한 모델에서는 f1(pos)는 baseline model의 f1(pos)보다 올라갔는데 f1(neg)는 baseline model의 f1(neg)보다 떨어질수 있습니다. 이러한 상황에서는 ROC AUC score라는 지표를 통해 평가할수 있습니다. 더 좋은 모델은 Baseline model의 ROC AUC보다 높은 점수를 냅니다.

우리가 baseline model로 만든 DecisionTreeClassifier()의 ROC AUC를 보겠습니다.

```py
from sklearn.metrics import roc_auc_score

y_val_pred_proba_dt = clf_dt.predict_proba(processed_X_val)[:,1]

roc_auc_score_dt = roc_auc_score(y_val, y_val_pred_proba_dt)
roc_auc_score_dt
```


#### Precision-Recall AUC Score

Target 값의 분포가 균등하다면 ROC AUC만 보고 모델들을 비교하며 평가할수 있습니다. 하지만 불균등하다면 ROC AUC만 보고 판단할수 없습니다. 필요에 따라 Precision-Recall curve AUC도 함께 보아야 합니다.

```py
from sklearn.metrics import average_precision_score as pr_auc_score

pr_auc_score_dt = pr_auc_score(y_val, y_val_pred_proba_dt)
pr_auc_score_dt
```


## RandomForest()

성능이 더 좋다고 판단되는 RandomForest()라는 모델을 만들어보고 학습시키겠습니다. 그리고 Holdout validation으로 예측한 결과를 가지고 Baseline Model보다 더 좋은 모델인지 평가지표를 통해 알아보겠습니다.

```py
from sklearn.ensemble import RandomForestClassifier

clf_rf = RandomForestClassifier(n_jobs=-1, random_state=1)
clf_rf.fit(processed_X_train, y_train);

y_train_pred_rf = clf_rf.predict(processed_X_train)
y_val_pred_rf = clf_rf.predict(processed_X_val)

print(classification_report(y_val, y_val_pred_rf))

```


Randomforest()를 적용해본 결과, f1(pos) & f1(neg) 모두 상승하였습니다. f1(pos) & f1(neg) 둘다 상승했을 경우 ROC AUC도 당연히 상승했을 것입니다. 확인해보겠습니다.

```py
from sklearn.metrics import roc_auc_score

y_val_pred_proba_rf = clf_rf.predict_proba(processed_X_val)[:,1]

roc_auc_score_rf = roc_auc_score(y_val, y_val_pred_proba_rf)
roc_auc_score_rf
```

이번엔 precision-recall AUC score를 보겠습니다.
```py
from sklearn.metrics import average_precision_score as pr_auc_score

pr_auc_score_rf = pr_auc_score(y_val, y_val_pred_proba_rf)
pr_auc_score_rf

```

이것 역시 baseline model의 score보다 상승했습니다. 일단은 RandomForestClassifier()가 DecisionTreeClassifier()보다 성능이 더 좋은 모델로 보입니다.

## RandomForest() with Optimized Hyperparameters (A.K.A Hyperparameter Tuning)

이번에는 RandomForest()를 활용하여 hyperparameter를 튜닝을 통해 튜닝되지 않은 RandomForest()보다 성능이 더 좋은 모델을 만들어 보겠습니다.

Python Code

```py
from sklearn.model_selection import RandomizedSearchCV

rf = RandomForestClassifier()
param_dists = {"max_depth":[5,10,15,30,50,100],
"max_features":['sqrt', 'log2'],
"min_samples_leaf": [1,2,3,4,5,10,20,30,50]}

clf_rf_tuning = RandomizedSearchCV(
estimator = rf,
param_distributions = param_dists,
n_iter = 10,
n_jobs = -1,
cv = 5,
scoring = "roc_auc"
)

clf_rf_tuning.fit(processed_X_train, y_train);

cv_results = pd.DataFrame(clf_rf_tuning.cv_results_)
```

```py
pd.DataFrame(cv_results[cv_results["rank_test_score"]==1]["params"].tolist()[0], index = [0])
```


#### 위에서 RandomizedSearchCV를 통해 얻은 Optimal Hyperparameter들을 기반으로 RandomForest()를 튜닝하겠습니다.

Python Code

```py
from sklearn.ensemble import RandomForestClassifier

clf_rf_tunned = RandomForestClassifier(n_jobs=-1,
random_state=1,
min_samples_leaf=10,
max_features="log2",
max_depth=50
)

clf_rf_tunned.fit(processed_X_train, y_train);

y_train_pred_rf_tunned = clf_rf_tunned.predict(processed_X_train)
y_val_pred_rf_tunned = clf_rf_tunned.predict(processed_X_val)

```

Hyperparameter를 튜닝한 RandomForest 모델을 돌려보니 다음과 같은 결과가 나왔습니다. f1(neg)는 튜닝하지 않은 모델에 비해서 올랐고 f1(pos)는 튜닝하지 않은 모델에 비해서 떨어졌다. 앞서 말했듯이 이러한 경우에서는 필히 ROC AUC score를 보아야합니다.

```py
print(classification_report(y_val, y_val_pred_rf_tunned))
```


앞서 우리가 튜닝하지 않은 RandomForest 모델을 사용했을때 ROC AUC Score는 0.6665가 나왔었습니다.. 아래에서 볼수있듯이 튜닝한 모델의 score가 더 좋게나왔습니다. 튜닝한 RandomForest 모델이 분류를 더 잘하는 모델이라는 뜻입니다.

```py
from sklearn.metrics import roc_auc_score

y_val_pred_proba_rf_tunned = clf_rf_tunned.predict_proba(processed_X_val)[:,1]

roc_auc_score_rf_tunned = roc_auc_score(y_val, y_val_pred_proba_rf_tunned)
roc_auc_score_rf_tunned
```

튜닝한 모델의 Precision-Recall AUC Score도 한번 보겠습니다. 튜닝하지 않았던 모델의 score는 0.355였습니다. 이역시 튜닝한 모델의 score가 더 높게나왔습니다.

```py
from sklearn.metrics import average_precision_score as pr_auc_score

pr_auc_score_rf_tunned = pr_auc_score(y_val, y_val_pred_proba_rf_tunned)
pr_auc_score_rf_tunned
```

### Optimal Threshold

마지막으로 Threshold를 정할차례입니다. 우리가 보아야 할 지표는 ROC Curve와 Precision-Recall Curve입니다. Target의 분포가 균등하다면 ROC Curve만 보면 되지만, 불균등할 경우 상황에 따라 Precision-Recall Curve도 볼수있습니다.

ROC Curve에서의 최적의 Threshold와 Precision-Recall Curve에서 최적의 Threshold는 다릅니다. 따라서 우리는 주어진 문제에 따라서 두 Curve중 하나를 선택하여 최적의 Threshold를 선택해야합니다.

우리의 문제에서는 ROC Curve에서의 최적의 Threshold를 선택하겠습니다.

참고자료: https://towardsdatascience.com/on-roc-and-precision-recall-curves-c23e9b63820c

```py
from sklearn.metrics import roc_curve

fpr, tpr, thresholds = roc_curve(y_val, y_val_pred_proba_rf_tunned)

roc = pd.DataFrame({
'FPR(Fall-out)': fpr,
'TPRate(Recall)': tpr,
'Threshold': thresholds
})

roc
```


#### ROC Curve를 그려보겠습니다.

ROC Curve의 y축은 TPR 그리고 x축은 (1-TNR) 또는 FPR입니다. ROC Curve에서 TPR과 TNR의 합이 가장 큰 부분의 우리의 optimal한 threshold가 될것입니다.

```py
plt.scatter(fpr, tpr)
plt.title('ROC curve')
plt.xlabel('FPR(Fall-out)')
plt.ylabel('TPR(Recall)');
```


#### 최적의 Threshold값 찾기

최적의 threshold는 0.2687로 보입니다.

```py
# threshold 최대값의 인덱스, np.argmax()
optimal_idx = np.argmax(tpr - fpr)
optimal_threshold = thresholds[optimal_idx]

print('idx:', optimal_idx, ', threshold:', optimal_threshold)
```

아래 그림을 보면 Threshold가 0.2687일때 TPR과 TNR의 합이 가장 큽니다.

```py
plt.plot(thresholds, tpr-fpr)
plt.axvline(x=0.26873370568955457, color="red")
plt.xlabel("Thresholds")
plt.ylabel("TPR + TNR")
plt.show()
```

```py
# 우리가 찾아낸 Optimal Threshold를 기반으로 예측값들을 변경시켜줍니다.

y_val_pred_proba_rf_tunned_optimized = y_val_pred_proba_rf_tunned >= optimal_threshold

y_val_pred_rf_tunned_optimized = np.multiply(y_val_pred_proba_rf_tunned_optimized, 1)
```

#### Hyperparameter를 튜닝한 RandomForest() 모델의 Threshold를 optimize 하였습니다.

> 이제 다시 평가지표들을 살펴보겠습니다.

> Threshold를 0.5에서 0.2687로 이동시켰습니다.

```py
print(classification_report(y_val, y_val_pred_rf_tunned_optimized))

```


> f1(neg)는 조금 떨어졌지만 f1(pos)는 많이 올라갔습니다.

> 이 현상을 풀이해 보자면, Threshold를 이동하기 전에는 negative에 대한 값들은 모델이 잘 예측하였습니다. 하지만 positive에 대한 값들은 상대적으로 잘 예측하지 못하였습니다. 하지만 Threshold를 이동시킴으로써 negative를 예측하는 성능을 조금 떨어뜨리면서 positive를 예측하는 성능을 많이 올린것이라고 볼수있습니다.

### 분포그림을 통해 최종모델을 그려보겠습니다.

--> **Threshold:** 빨간 점선

--> **파란색 분포:** Inspection Fail (positive) 또는 (1)

--> **오렌지색 분포:** Inspection Pass (negative) 또는 (0)

> Threshold를 기준으로 오른쪽은 Inspection Fail했다고 예측하는 것이고

> Threshold를 기준으로 왼쪽은 Inspection Pass했다고 예측하는 것입니다.

Python Code

```py
prob_n_actual_comps_rf = pd.DataFrame({"pred_proba":y_val_pred_proba_rf_tunned, "target_values":y_val})

actually_fail_rf = prob_n_actual_comps_rf[prob_n_actual_comps_rf["target_values"]==1]
actually_pass_rf = prob_n_actual_comps_rf[prob_n_actual_comps_rf["target_values"]==0]

plt.figure(figsize=(8,6))
plt.title("RandomForestClassifier() Model")
plt.yscale('log')
plt.hist(actually_fail_rf["pred_proba"], bins=20, alpha=0.5, label="Inspection Fail")
plt.hist(actually_pass_rf["pred_proba"], bins=20, alpha=0.5, label="Inspection Pass")
plt.axvline(x=0.2687, color="red", linestyle="--", label="threshold")
plt.xlabel("Probabilities")
plt.ylabel("Log(Frequency)")
plt.legend()

plt.show()

```


# 결론

조금더 사실적인 그림은 아래 그림입니다. 아래 그림에서 볼 수 있듯이, 현 모델이 Baseline 모델보다는 좋은 모델이지만, 이 모델을 실상에서는 바로 쓸수 없을것 같습니다. 왜냐하면, 잘못 예측하고 있는 부분이 한눈에 봐도 많이 보이기 때문입니다.

Python Code

```py
prob_n_actual_comps_rf = pd.DataFrame({"pred_proba":y_val_pred_proba_rf_tunned, "target_values":y_val})

actually_fail_rf = prob_n_actual_comps_rf[prob_n_actual_comps_rf["target_values"]==1]
actually_pass_rf = prob_n_actual_comps_rf[prob_n_actual_comps_rf["target_values"]==0]

plt.figure(figsize=(8,6))
plt.title("RandomForestClassifier() Model")
#plt.yscale('log')
plt.hist(actually_fail_rf["pred_proba"], bins=20, alpha=0.5, label="Inspection Fail")
plt.hist(actually_pass_rf["pred_proba"], bins=20, alpha=0.5, label="Inspection Pass")
plt.axvline(x=0.2687, color="red", linestyle="--", label="threshold")
plt.xlabel("Probabilities")
plt.ylabel("Frequency")
plt.legend()

plt.show()

```