[프로젝트] 이탈할 고객을 예측하여 수익 감소 방지하기

최근 EDA 역량의 부족함을 느끼고 Kaggle EDA 필사를 하고 있습니다. 타이타닉 EDA 필사를 하며 공부한 내용과 핸즈온 Chapter2를 공부한 후 참고하여 고객 이탈 예측 프로젝트를 진행했습니다. (Kaggle의 데이터 이용)

 

 

고객 이탈 예측

 

순서

 

1. 문제를 정의하고 큰 그림 보기

2. 데이터 가져오기

3. 데이터로부터 통찰을 얻기 위해 탐색하고 시각화

4. 머신러닝 알고리즘을 위한 데이터 준비

5. 모델 선택과 훈련

6. 모델 튜닝

7. 솔루션 제시

8. 론칭, 모니터링, 시스템 유지 보수 

 

💻 사용한 언어: Python

1. 문제를 정의하고 큰 그림 보기

모델을 학습시켜 새로운 데이터가 주어졌을 때 이탈할 고객을 예측하는 것이 목표입니다. 이탈할 것이라고 분류된 고객들의 주요 특성을 파악하고 미래 이탈 방지를 통해 고객 관계 유지 강화 및 수익 감소를 방지하기 위한 프로젝트입니다. 이 문제는 레이블이 있는 지도 학습이고 이탈할 것이지 아닌지를 분류하는 문제입니다. 데이터가 7043개 행, 21개의 변수로 작기 때문에 일반적인 배치 학습을 이용합니다. 실제 이탈하지 않을 고객을 이탈할 것이라고 예측하면 실제 이탈하지 않을 고객에게 투자함으로써 발생하는 손해액이 생기기 때문에 재현율에 초점을 두고 성능을 측정합니다. 

 

2. 데이터 가져오기 

데이터는 어느 한 통신 회사의 데이터로 각 고객이 가입한 서비스(전화, 다중회선, 인터넷, 온라인 보안, 온라인 보안, 장치 보호, 기술 지원, 스트리밍 TV 및 영화), 고객 계정 정보(고객ID, 계약, 결제 방법, 디지털 청구서, 월별 요금 및 총 요금), 고객에 대한 인구 통계 정보(성별, 연령 범위, 파트너 및 부양 가족이 있는지 여부) 변수들이 있습니다. 먼저, 데이터의 구조를 훑어봅니다. head 메서드를 이용하면 21개의 변수, 5개의 행이 나타납니다. 

import pandas as pd

pd.set_option('display.max_columns', 500)  # 모든 컬럼을 보기 위해 설정

df = pd.read_csv("../../data/WA_Fn-UseC_-Telco-Customer-Churn.csv")
df.head()

info 메서드를 이용하면 각 변수별로 결측치 수, 타입을 확인할 수 있습니다. 이 결과를 통해 TotalCharges는 수치형 변수인데 명목형 변수로 인식하였고 SeniorCitizen은 명목형 변수인데 수치형 변수로 되어있습니다. 이 변수들의 타입을 먼저 수정하겠습니다.

df.info()

TotalCharges에 빈 문자열이 있기 때문에 수치형으로 변환하면 에러가 나타납니다. errors = "coerce"를 추가하여 결측치로 바뀌도록 합니다. 

df["SeniorCitizen"] = df["SeniorCitizen"].astype("object")
df["TotalCharges"] = pd.to_numeric(df["TotalCharges"],errors = "coerce")

히스토그램은 각 숫자형 특성 별로 범위(수평축)에 속한 샘플 수(수직축)를 나타내어 데이터를 빠르게 검토할 수 있습니다. hist 메서드를 이용하여 수치형 변수들의 분포를 확인합니다. 이때, hist 메서드는 matplotlib을 사용하고 결국 화면에 그래프를 그리기 위해 사용자 컴퓨터의 그래픽 백엔드를 필요로 합니다. 그래서 사용할 백엔드를 지정해줘야 하는데 %matplotlib inline을 사용하면 주피터 자체의 백엔드를 사용하도록 설정합니다. 

%matplotlib inline
import matplotlib.pyplot as plt

df.hist(figsize = (10,8),bins = 10)

수치형 변수들의 히스토그램

 

각 수치형 변수들의 X축을 보면 스케일이 다른 것을 알 수 있습니다. TotalCharges 변수는 꼬리가 두껍습니다. 머신러닝 알고리즘은 종 모양의 분포가 되면 패턴을 더욱 쉽게 찾을 수 있습니다. 데이터 분포를 정규분포나 가우시안 분포를 따르도록 변형하는 것이 좋습니다.

 

구체적으로 살펴보기 전에 먼저 테스트 데이터를 따로 생성합니다. 이때 테스트 데이터는 절대 들여다보면 안됩니다. 테스트 데이터로 일반화 오차를 추정하며 매우 낙관적인 추정이 되어 시스템을 실제로 론칭했을 때 기대한 성능이 나오지 않을 수 있습니다. 이를 "데이터 스누핑 편향"이라고 합니다. 테스트 데이터를 전체 데이터의 20% 정도로 사용하겠습니다.

 

따로 데이터 분석을 통해 알아낸 사실로 회사와 계약해 온 기간이 짧은 고객이 많이 이탈한다는 것을 참고하여 테스트 데이터를 만들겠습니다. (데이터 분석은 통계 검정을 적용하여 진행했습니다.) 즉, 중요 변수인 고객이 회사와 계약해 온 기간을 계층으로 계층적 샘플링을 하여 테스트 데이터를 생성할 것입니다. 순수한 무작위 샘플링 방식은 큰 데이터셋에서는 일반적으로 괜찮지만, 그렇지 않다면 샘플링 편향이 생길 가능성이 있습니다. 테스트 세트가 전체를 대표하도록 각 계층에서 올바른 수의 샘플을 추출하는 것이 좋습니다. (이때, 전문가를 통해 알게된 중요 변수가 있다면 그 변수를 사용하여 테스트 데이터를 생성하는 것도 좋습니다.) 계층의 중요도를 추정하는데 편향이 발생하지 않기 위해 계층별로 데이터셋에 충분한 샘플 수가 있어야 합니다. 

 

tenure 변수에 0의 값을 가지는 경우가 11건 존재합니다. 이는 계약한 지 1달도 안된 신규 고객이므로 제외하고 모델링하겠습니다. 

 

df = df[df["tenure"]!=0].reset_index(drop=True)
df.shape  # 7032, 22 (11건 삭제됨.)

df["tenure_cat"] = pd.cut(df["tenure"], bins = [0,24,48,60,72], labels = [1,2,3,4])
df["tenure_cat"].hist()

 

고객과 계약해 온 기간을 기반으로 StratifiedShuffleSplit 을 이용하여 계층적 샘플링을 진행합니다. 

from sklearn.model_selection import StratifiedShuffleSplit

split = StratifiedShuffleSplit(n_splits = 1, test_size = 0.2, random_state = 20171490)

for train_index, test_index in split.split(df, df["tenure_cat"]):
    X_train = df.loc[train_index]
    X_test = df.loc[test_index]

샘플링한 후 각 계층별로 샘플 수가 비슷하게 추출되었는지 확인하겠습니다.

X_train["tenure_cat"].value_counts()/len(X_train)
X_test["tenure_cat"].value_counts()/len(X_test)

 

비슷한 비율로 샘플링 된 것을 확인할 수 있습니다. 샘플링한 후 이 카테고리 화한 변수는 제거합니다.

for set_ in (X_train, X_test):
    set_.drop("tenure_cat",axis = 1, inplace = True)

 

3. 데이터로부터 통찰을 얻기 위해 탐색하고 시각화

데이터를 탐색할 때는 훈련 데이터만 이용합니다. 훈련 데이터를 손상시키지 않기 위해 복사본을 만들고 데이터 클래스 분포를 확인해보겠습니다.

from matplotlib import font_manager, rc
font_name = font_manager.FontProperties(fname="c:/Windows/Fonts/malgun.ttf").get_name()
rc('font', family=font_name)

data = X_train.copy()

plt.figure(figsize = (8,6))
splot = sns.countplot(data["Churn"], edgecolor="black",alpha=0.7)  # alpha를 이용하여 색 조정 
plt.title("지난 달에 고객 이탈 유무")


for p in splot.patches:
    splot.annotate(format(p.get_height(), '.0f'), 
                   (p.get_x() + p.get_width() / 2., p.get_height()), 
                   ha = 'center', va = 'center', 
                   size=15,
                   xytext = (0, -12), 
                   textcoords = 'offset points')

class 분포

데이터 불균형이 존재합니다. 데이터 불균형을 해소하지 않으면 모델의 성능이 좋지 않을 것입니다.

 

변수들의 상관관계를 확인합니다. 명목형 변수와 수치형 변수들 사이의 상관관계는 pointbiserialrPoint Biserial Correlation를 이용하여 검사합니다.

from scipy import stats

num_features = list(data.select_dtypes(include=[np.number]))

# Point Biserial Correlation table
pb_list = []
print("Significant Point Biserial Correlations")
print("{0:23} |   {1:4}|  {2:}".format("Factor", "r", "p"))
for i in num_features:
    r, p = stats.pointbiserialr(data['Churn_b'], data[i])
    if p < .01:
        print("{0:23} | {1:5.2f} | {2:.2f}".format(i, r, p))
        pb_list.append(i)

 

수치형 독립변수와 종속변수 간 상관관계

sns.heatmap(data.corr(), annot = True)

Churn과 상관계수값은 tenure이 가장 높고 그다음은 TotalCharges, MonthlyCharges입니다. 수치형 변수들 사이의 상관관계도 heatmap을 이용하여 확인해보겠습니다. 

 

수치형 독립변수들 사이의 상관관계

수치형 독립변수들 간의 상관관계를 확인해보면 TotalCharges와 Tenure가 0.83, MonthlyCharges와 0.65로 꽤 높은 양의 상관관계를 가지고 있습니다. 이러한 경우는 다중공선성이 의심되므로 다중공선성이 발생하는지 확인해봐야 합니다. 

 

from statsmodels.formula.api import ols

model = ols('y ~ tenure + MonthlyCharges+TotalCharges', data)
res = model.fit()
res.summary()

from statsmodels.stats.outliers_influence import variance_inflation_factor

pd.DataFrame({'컬럼': column, 'VIF': variance_inflation_factor(model.exog, i)} 
             for i, column in enumerate(model.exog_names)
             if column != 'Intercept')  # 절편의 VIF는 구하지 않음.

다중공선성 확인 결과

tenure, MonthlyCharges, TotalCharges 변수 모두 통계적으로 유의한 변수이지만 TotalChareges의 경우 다중공선성이 높습니다. 다중공선성은 독립변수들 사이의 상관 정도가 높아 데이터 분석 시 부정적인 영향을 미치는 것을 말합니다. 다중공선성은 변수들 각각의 설명력을 약하게 만듭니다. 다중공선성은 계수가 통계적으로 유의미하다면 다중공선성이 크더라도 특별히 대처하지 않아도 되고 유의미하지 않다면 변수들을 더하거나 빼서 새로운 변수를 만듭니다. 또는 모형에서 해당 변수를 제거하는데 이는 자료의 다양성을 해치고 분석하려던 가설이나 이론에 영향을 미칠 수 있기 때문에 자제합니다. 하지만 (회귀모델처럼) 설명력있는 모델에서는 다중공선성이 중요하지만 머신러닝 모델에서는 예측력이 좋으면 상관 쓰지 않아도 된다는 말도 있더라구요.

 

TotalCharges는 총 요금으로 MonthlyCharges와 비슷한 속성을 갖고 있는 변수라고 생각합니다. 이 변수는 제외하고 사용하겠습니다. (사실 TotalCharges를 변수를 제외해서 사용하는 것이 저의 판단에서 좋을 것 같다고 생각하지만 모델에 적용했을 때 더 좋은 결과가 나올까봐 TotalCharges를 넣은 모델도 돌려봤답니닷...경험하기 전에 더 잘 맞을 것이라고 보장할 수 있는 건 없으니까요!)

 

명목형 변수와 명목형 변수 간의 상관관계는 카이제곱을 이용하여 독립성 검정을 실시합니다. 독립성 검정의 귀무가설은 "두 변수가 서로 독립이다."입니다.

from scipy.stats import chi2_contingency

p_feature = []

for col_ in category_features:
    chi, p, dof, expected = chi2_contingency(pd.crosstab(df[col_], df.Churn))
    
    if p<=0.05:
        p_feature.append(col_)

        print(col_,
              "Chi2 Statistic: {}".format(round(chi,2)),
              "p-value (0.05): {}".format(round(p,2)),
              "degree of freedom: {}".format(dof),
              "",
              "expected value: \n{}".format(pd.DataFrame(np.round_(expected,2), columns=['No','Yes'])),
              "",
              "original value: \n{}".format(pd.DataFrame(pd.crosstab(df[col_], df.Churn))),

              "-----------------------------------", sep = "\n" )
    

만약 유의 확률이 유의수준 0.05보다 높으면 귀무가설을 채택하고 낮으면 기각하고 상관있는 변수라고 판단합니다. 유의확률이 유의 수준보다 높았던 변수는 gender과 PhoneService였습니다. 

 

이를 기반으로 패턴을 잘 찾을 수 있는 파생변수를 만들어보겠습니다. 파생변수는 때때로 특징 공학에서 중요하게 작동합니다. 온라인 보안, 온라인 백업, 장치 보호, 기술지원 서비스 가입 변수들은 모두 인터넷 서비스 사용자들의 보조 서비스입니다. 가입한 보조 서비스 개수 변수를 생성하고 스트리밍 영화와 TV 변수를 합쳐서 스트리밍 서비스 가입 개수 변수를 생성하겠습니다. 

 

4. 머신러닝 알고리즘을 위한 데이터 준비

머신러닝 알고리즘을 위한 데이터를 준비하는 작업은 자동화를 하는 것이 좋습니다. 자동화는 어떤 데이터셋에서도 데이터 변환을 손쉽게 할 수 있고 여러 가지 데이터 변환을 쉽게 시도해볼 수 있고 어떤 조합이 가장 좋은지 편리하게 확인하는데 도움을 줍니다. 

 

먼저 X 데이터와 y 데이터를 분리합니다. 

X = X_train.drop(["Churn"],axis = 1)
y = X_train["Churn"]

y.replace("No",0,inplace = True)
y.replace("Yes",1,inplace = True)

대부분의 머신러닝 알고리즘은 결측치를 다루지 못하므로 이를 처리할 수 있는 함수를 생성합니다. 결측치는 행 또는 열을 제거하거나 다른 값으로 대체할 수 있습니다. 본 데이터는 각 사람들의 데이터이고 TotalCharges에 결측치가 있었는데 이 결측치는 계약 기간이 (1개월도 채 안된) 0개월인 고객들이기 때문에 제거하도록 하겠습니다. 만약 다른 변수에 결측치가 존재하는 새로운 데이터가 주입된다면 알람이 울리도록 만드는 것도 좋을 것 같습니다. 

 

훈련 데이터와 검증 데이터를 나누어 훈련 데이터를 통해 학습하고 검증 데이터를 통해 파라미터를 튜닝한 후 테스트 데이터에 가장 좋은 모델을 적용해보도록 하겠습니다. 

X_train, X_valid, y_train, y_valid = train_test_split(X,y,stratify = y,test_size = 0.2, shuffle = True, random_state = 20171490)

 

앞서 데이터를 탐색하면서 수치형 변수들의 스케일이 달랐습니다. 이상치에 영향을 받지 않는 Robust Scaler를 이용하여 스케일을 맞춰주겠습니다. sklearn의 pipeline 메서드를 사용하면 간편하게 데이터를 전 처리할 수 있습니다. 필요한 변수만 추출하는 class를 만들고 스케일을 맞춰주기 위해 RobustScaler를 적용하는 num_pipeline 을 만듭니다. 

 

class DataFrameSelector(BaseEstimator, TransformerMixin):
    def __init__(self, attribute_names):
        self.attribute_names = attribute_names
    def fit(self, X, y=None):
        return self
    def transform(self, X):
        return X[self.attribute_names]

num_pipeline = Pipeline([
        ("select_numeric", DataFrameSelector(["tenure", "MonthlyCharges"])),
        ('scaler', RobustScaler()),
    ])

다음은 파생변수를 생성하기 위한 class를 정의합니다. 보조 서비스 가입 개수와 스트리밍 서비스 가입 개수 변수를 생성하고 기존의 변수들은 사용하지 않겠습니다. 인터넷 서비스 사용자 중 해당 서비스에 가입했는지 유무보다 보조서비스 가입 개수 변수가 더 의미있을 것이라고 생각했기 때문입니다. 이를 확인하기 위해 class를 정의하기 전에 시각화하여 유의미한 변수일 것인지 추측해보겠습니다. 아래 시각화를 보면 보조서비스 가입 개수가 없는 경우, 스트리밍 가입을 안 한 경우 이탈한 고객들이 많습니다. 스트리밍 서비스는 TV와 영화 모두 가입한 고객들이 둘 중 하나의 서비스만 이용하는 고객보다 더 많이 이탈합니다. 

 

이탈하는 고객들과 그렇지 않은 고객들의 차이가 확연히 보입니다. 이 변수를 사용하겠습니다. 

class CombinedAttributesAdder(BaseEstimator, TransformerMixin):

    def __init__(self, attribute_names): # *args 또는 **kargs 없음
        self.attribute_names = attribute_names
    
    def fit(self, X, y=None):
        return self  
        
    def transform(self, X):
        
        df = X[self.attribute_names]
        df.replace("Yes",1,inplace=True)
        df.replace("No",0,inplace=True)
        df.replace("No internet service",0,inplace=True)
    
        df["SupplementaryService"] = df["OnlineSecurity"]+df["OnlineBackup"]+df["DeviceProtection"]+df["TechSupport"]
        df["Streaming"] = df["StreamingTV"]+df["StreamingMovies"]

        return df.drop(self.attribute_names, axis = 1)

범주형 데이터 중 값이 2개씩 있는 경우는 OneHotEncoding을 적용하면 변수 1개로 나타낼 수 있음에도 불구하고 2개로 나타내게 됩니다. 메모리 절약을 위해 범주형 데이터 중 값이 2개씩 있는 경우는 OrdinalEncoder 클래스를 값이 3개씩 있는 경우는 OneHotEncoder 클래스를 사용하여 인코딩을 적용하겠습니다. 이를 위해 따로 변수를 추출하는 코드를 작성합니다. 이때, "customerID" 변수는 독립변수가 아니기 때문에 제외합니다. 

cate_features2 = []
cate_features_not2 = []

for column in X_train.columns:
    if X_train[column].dtypes==object:
        if len(X_train[column].value_counts())==2:
            cate_features2.append(column)
        else:
            cate_features_not2.append(column)
        
print(cate_features2)
print(cate_features_not2)

cate_features_not2.remove("customerID")

각 목적에 맞는 Pipeline을 만듭니다. 

cate_1hot_pipeline = Pipeline([
        ("select_cate_1hot", DataFrameSelector(['MultipleLines', 'InternetService', 'Contract', 'PaymentMethod'])),
        ('1hot', OneHotEncoder(sparse=False)),  # sparse=False 밀집행렬로 만들어야함. 
    ])

cate_ord_pipeline = Pipeline([
        ("select_cate_ord", DataFrameSelector(cate_features2)),
        ('ordinary', OrdinalEncoder()),
    ])
    
preprocess_pipeline = FeatureUnion(transformer_list=[
        ('attribs_adder', CombinedAttributesAdder(col_names)),
        ("num_pipeline", num_pipeline),
        ("cat_1hot_pipeline", cate_1hot_pipeline),
        ("cat_ord_pipeline", cate_ord_pipeline)])


X_prepared = preprocess_pipeline.fit_transform(X_train)
X_valid_prepared = preprocess_pipeline.fit_transform(X_valid)    

 

5. 모델 선택 및 훈련

훈련 데이터, 검증 데이터에 LogisticRegression, RandomForest, AdaBoost, XGB 분류 모델을 적용하고 성능을 확인합니다. 모델을 적용하고 성능을 확인하는 코드는 이후에도 반복하여 사용할 코드이기 때문에 함수로 정의합니다.

from sklearn.metrics import accuracy_score, precision_score, recall_score, roc_auc_score

def Modeling_result(train_x, train_y, valid_x, valid_y):
    
    models = [LogisticRegression(), RandomForestClassifier(), AdaBoostClassifier(), XGBClassifier()]
    names = ["Logistic","RF","Ada","XGB"]

    train_acc = []
    train_precision = []
    train_recall = []
    train_roc = []

    valid_acc = []
    valid_precision = []
    valid_recall = []
    valid_roc = []


    for model, name in zip(models, names):

        model.fit(train_x,train_y)
        model.predict(valid_x)

        train_acc.append(accuracy_score(train_y,model.predict(train_x)))
        train_precision.append(precision_score(train_y,model.predict(train_x)))
        train_recall.append(recall_score(train_y,model.predict(train_x)))
        train_roc.append(roc_auc_score(train_y, model.predict(train_x)))

        valid_acc.append(accuracy_score(valid_y, model.predict(valid_x)))
        valid_precision.append(precision_score(valid_y,model.predict(valid_x)))
        valid_recall.append(recall_score(valid_y,model.predict(valid_x)))
        valid_roc.append(roc_auc_score(valid_y, model.predict(valid_x)))
        
    return pd.DataFrame({"Model":names, "Train_Accuracy":train_acc, "Train_Precision":train_precision, "Train_Recall":train_recall, "Train_AUC":train_roc,
             "Valid_Accuracy":valid_acc, "Valid_Precision":valid_precision, "Valid_Recall":valid_recall, "Valid_AUC":valid_roc})

gender_right = Modeling_result(X_prepared, y_train, X_valid_prepared, y_valid)
gender_right

RF와 XGB는 과대 적합이 발생했습니다. 과대적합이 발생하면 feature 수를 줄이거나 규제를 줄여서 해결할 수 있습니다. EDA를 통해 gender, PhoneServie는 종속변수와 독립인 변수임을 확인했습니다. 이 변수들을 제외하여 features 수를 줄이는 방안을 택하겠습니다.

cate_ord_no_cor_pipeline = Pipeline([
("select_cate_ord", DataFrameSelector(['SeniorCitizen', 'Partner', 'Dependents','PaperlessBilling'])),
('ordinary', OrdinalEncoder()),
])

preprocess_no_pipeline = FeatureUnion(transformer_list=[
        ('attribs_adder', CombinedAttributesAdder(col_names)),
        ("num_pipeline", num_pipeline),
        ("cat_1hot_pipeline", cate_1hot_pipeline),
        ("cat_ord_pipeline", cate_ord_no_cor_pipeline)])

X_valid_prepared_no = preprocess_no_pipeline.fit_transform(X_valid)
X_prepared_no = preprocess_no_pipeline.fit_transform(X_train)

no_cor_result = Modeling_result(X_prepared_no, y_train, X_valid_prepared_no, y_valid)
no_cor_result

train data에 대해서 성능이 조금 낮아진 것처럼 보이지만 크게 차이가 없습니다. 하지만 성능이 좋지 않습니다. 앞서 말씀드렸다시피 데이터 불균형을 해소하지 않았기 때문으로 보입니다. 데이터 불균형을 해소하기 위해 SMOTE를 이용하여 업샘플링을 진행하겠습니다. 업샘플링에 사용할 데이터는 성능에 크게 차이가 없다면 메모리 절약과 복잡도 감소를 위해 feature수가 더 적은 데이터를 이용할 것입니다. 

 

from imblearn.over_sampling import SMOTE

sm = SMOTE(random_state = 20171490)

X_res_train, y_res_train = sm.fit_sample(X_prepared_no, y_train)
print(y_train.value_counts(), np.bincount(y_res_train))

upsampling_result = Modeling_result(X_res_train, y_res_train, X_valid_prepared_no, y_valid)
upsampling_result

데이터 불균형을 해소했더니 정밀도(Precision)는 낮아졌지만 재현율(Recall)이 올라갔습니다. 다음은 따로 만들어둔 검증 데이터셋이 아닌 훈련 데이터와 검증 데이터셋 전체를 이용하여 교차검증을 적용하겠습니다. 따로 검증 데이터셋을 만들어서 검증을 확인하는 것도 잘 작동하지만 만들어둔 검증 데이터셋에서만 잘 동작하는 모델을 만들 수도 있기 때문에 K-fold cross validation(K 겹 교차검증)을 이용합니다. 이렇듯 평가에 사용되는 데이터 편중을 막기 위해 매번 다른 폴드를 선택해 평가에 사용하고 나머지 폴드들의 데이터를 훈련하는 방법을 사용하여 좀 더 일반화된 모델을 만듭니다. 하지만 모델을 훈련하고 평가하는데 시간이 오래 걸린다는 단점이 있습니다. 

X_res, y_res = sm.fit_sample(preprocess_no_pipeline.fit_transform(X), y)
print(y_train.value_counts(), np.bincount(y_res))

models = [LogisticRegression(), RandomForestClassifier(), AdaBoostClassifier(), XGBClassifier()]
names = ["Logistic","RF","Ada","XGB"]
cv_scores = []

acc = []
precision = []
recall = []
roc = []

for model, name in zip(models, names):
    print(name)
    start = time.time()
    
    # method{'predict’,'predict_proba’, 'predict_log_proba’, 'decision_function'}, default=’predict’
    cv_predict = cross_val_predict(model, X_res, y_res, cv=10)
    cv_scores.append(cross_val_score(model, X_res, y_res, cv=10,scoring='roc_auc'))

    acc.append(accuracy_score(y_res, cv_predict))
    precision.append(precision_score(y_res,cv_predict))
    recall.append(recall_score(y_res,cv_predict))
    roc.append(roc_auc_score(y_res, cv_predict))
    
pd.DataFrame({"Model":names, "Accuracy":acc, "Precision":precision, "Recall":recall, "AUC":roc})

plt.figure(figsize=(8, 4))

i = 1

for cv_score in cv_scores:
    
    plt.plot([i]*10, cv_score, ".")
    
    i+=1
    
plt.boxplot(cv_scores, labels=("Logistic","RF","Ada","XGB"))
plt.ylabel("roc_auc", fontsize=14)
plt.show()

 

6. 모델 튜닝

교차검증을 확인하면 RF가 가장 높고 비슷한 점수를 가지는 경우가 XGB입니다. 박스 플롯으로 그려보면 XGB의 표준편차가 매우 큰 것으로 보입니다. RF 모델로 모델 튜닝을 진행하겠습니다.

param_grid_RF = {
    "n_estimators" :[10,15,20,25], 
    "max_depth" : [5,10,15],
    "min_samples_split":[2,4,8,16],
    "max_features":["sqrt","auto","log2"],
    "class_weight" : ["balanced_subsample","balanced"]}

cv_RF_roc = GridSearchCV(RandomForestClassifier(),param_grid_RF, cv= cv_mo, scoring = "roc_auc",return_train_score=True)
cv_RF_roc.fit(X_res_train, y_res_train)
cv_RF_best = cv_RF_roc.best_estimator_

GridSearch를 통해 가장 좋은 성능을 가지는 파라미터를 적용하여 교차검증을 적용합니다. 

param = cv_RF_roc.best_params_
model_RF = RandomForestClassifier(class_weight = "balanced",max_depth = 15, max_features = 'sqrt', min_samples_split = 2, n_estimators = 25)

# method{'predict’,'predict_proba’, 'predict_log_proba’, 'decision_function'}, default=’predict’
cv_predict = cross_val_predict(model_RF, X_res, y_res, cv=10)

acc = accuracy_score(y_res, cv_predict)
precision = precision_score(y_res,cv_predict)
recall = recall_score(y_res,cv_predict)
roc =roc_auc_score(y_res, cv_predict)

pd.DataFrame([{"Model":"RF_gradient","Accuracy":acc, "Precision":precision, "Recall":recall, "AUC":roc}])

RF 모델을 튜닝한 후 결과를 확인해보면 GridSearch 하기 전보다 Recall이 0.01올랐고 Precision은 0.01 낮아졌습니다. 크게 성능이 달라지지 않았습니다. 

 

훈련 데이터셋과 검증 데이터셋 모두에서 검정 결과가 어떤 지 궁금하기 때문에 이것도 확인해보겠습니다.

param = cv_RF_roc.best_params_
model_RF = RandomForestClassifier(class_weight = "balanced", max_depth = 15, max_features = 'auto', min_samples_split = 2, n_estimators = 25)

model_RF.fit(X_res_train, y_res_train)
print(confusion_matrix(y_res_train,model_RF.predict(X_res_train)))
print(classification_report(y_res_train,model_RF.predict(X_res_train)))

print(confusion_matrix(y_valid,model_RF.predict(X_valid_prepared_no)))
print(classification_report(y_valid,model_RF.predict(X_valid_prepared_no)))

GridSearch 후 모델 결과

모델이 과대적합되었습니다. max_depth(트리의 최대 깊이)을 이용하여 모델을 규제하고 과대적합의 위험을 감소시키겠습니다.

param = cv_RF_roc.best_params_
model_RF = RandomForestClassifier(class_weight = "balanced", max_depth = 8, max_features = 'auto', min_samples_split = 2, n_estimators = 25)

model_RF.fit(X_res_train, y_res_train)
print(confusion_matrix(y_res_train,model_RF.predict(X_res_train)))
print(classification_report(y_res_train,model_RF.predict(X_res_train)))

print(confusion_matrix(y_valid,model_RF.predict(X_valid_prepared_no)))
print(classification_report(y_valid,model_RF.predict(X_valid_prepared_no)))

규제 적용 후 모델 결과

사실 저는 좀 더 좋은 결과가 나왔으면 좋겠습니다. 론칭하기에는 성능이 낮아 보입니다. AutoML을 이용하여 더 좋은 모델이 있는지 확인해보고 싶습니다. 

 

훈련 데이터 X, y를 합친 후 모델에 적용합니다.

from pycaret.classification import *

ml_data = pd.concat([pd.DataFrame(X_res),y_res],axis=1)

setup 메서드를 실행시키면 아래와 같이 자동으로 인식한 데이터셋의 설명서를 제공합니다. 이상이 없다면 엔터를 치고 바로 넘어갑니다. 관련 문서사실 이때 파생 변수를 생성하는 것이 아니라면 범주형 데이터를 수치형 데이터로 바꾸는 것, 결측치를 대체하는 것과 같은 전처리 과정은 하지 않아도 자동으로 해줍니다. 저는 이미 전처리한 데이터가 있고 파생변수를 만들기 때문에 파이프라인을 적용한 데이터를 사용하겠습니다.

clf = setup(data = ml_data, target = "Churn")

AUC를 기준으로 3개의 모델을 저장하겠습니다. compare_models 메서드를 적용하면 아래 결과와 같이 적용 모델들을 비교하여 출력해줍니다. 순위 정렬 기준을 AUC, 최종 선택될 모델 수를 3개로 선택하여 모델을 확인해보겠습니다. setup 설정에 맞춰 각 모델의 평균점수를 비교하여 가장 우수한 모델부터 순차적으로 보여줍니다. best로 선택된 모델은 catboost, xgboost, lightgbm 모델입니다. SVM 모델과 Ridge 모델의 AUC값은 0이라고 나오는데 왜 인지 이해가 되지 않네요ㅜㅜ...

 

직접 만든 RF 모형과 비교했을 때보다 더 안 좋은 성능을 가지네요. 

best_3 = compare_models(sort = 'AUC', n_select = 3)  # 모델 성능 비교 

AUC 값이 가장 좋은 모델 3개를 대상으로 tune_model을 이용하여 모델 튜닝을 통해 성능이 더 좋아지도록 만들겠습니다. optimize 파라미터를 이용하여 최적화하고자 하는 대상의 성능 지표를 선정할 수 있습니다. 

tuned_top3 = [tune_model(i) for i in best_3]
tuned_top3

모델 튜닝 결과 직접 만들었던 랜덤포레스트 모형보다 더 좋은 성능을 보이지 않습니다. 앙상블 모델과 스태킹 모델을 통해 보다 좋은 성능을 나타내도록 만들 수 있는지 확인해보겠습니다.  

blender = blend_models(estimator_list = tuned_top3,fold = 10, method = 'soft') 
blender

튜닝한 모델들로 앙상블 soft 방식을 적용했더니 Recall은 비슷하지만 Precision이 0.80으로 0.01올랐네요! 아래 시각화를 보면 훈련 데이터 점수가 점점 낮아지는 시각화를 보입니다. 

 

evaluate_model 메서드를 사용하면 Plot type을 클릭하여 원하는 플롯을 볼 수 있습니다. 정말 간편하게 잘 만든 라이브러리입니다. 이때, Feature importance, Decision Tree plot은 모델이 트리 기반이어야 사용 가능합니다. 

evaluate_model(estimator=blender)

stacker = stack_models(estimator_list=tuned_top3, fold = 10)
plot_model(estimator = stacker, plot = 'learning')

투표 기반 분류기와 스태킹 모델은 직접 만들었던 랜덤 포레스트 모델보다 더 좋은 성능을 보이지 않기 때문에 최종 모델로 랜덤포레스트를 사용하겠습니다.  

 

가장 처음 분리했던 text 데이터에 전처리 파이프라인을 적용합니다.

y_test = X_test["Churn"]
X_test = X_test.drop(["Churn"],axis = 1)

X_prepared_test = preprocess_no_pipeline.fit_transform(X_test)

y_test.replace("Yes",1,inplace = True)
y_test.replace("No",0,inplace = True)

model_RF = RandomForestClassifier(class_weight = "balanced",max_depth = 8, max_features = 'auto', min_samples_split = 2, n_estimators = 25, random_state = 42)
model_RF.fit(X_res, y_res)
y_pred = model_RF.predict(X_prepared_test)

acc = accuracy_score(y_test, y_pred)
precision = precision_score(y_test,y_pred)
recall = recall_score(y_test,y_pred)
roc =roc_auc_score(y_test, y_pred)

pd.DataFrame([{"Model":"RF_test","Accuracy":acc, "Precision":precision, "Recall":recall, "AUC":roc}])

테스트 셋 모델 적용 결과

 테스트 셋 혼동행렬

테스트 셋에서 성능 수치를 좋게 하려고 하이퍼파라미터를 튜닝하면 새로운 데이터에 일반화되기 어렵습니다. 교차 검증을 사용하여 측정한 것보다 조금 성능이 낮은 것이 보통이기는 하지만 성능이 꽤 차이가 나기 때문에 아쉬운 결과입니다. Precision이 낮은게 정말 아쉽네요,,

 

다른 파생변수를 사용하거나 더 적절한 모델을 선택해야 할 것 같습니다. 만약 데이터와 관련해서 예측하는데 중요한 사항들을 인지하고 있다면 보다 좋은 성능을 낼 수 있도록 전처리할 수 있지 않을까 하는 아쉬움이 있습니다.

 

(만약 이탈 가능성이 높은 고객을 추정하고 싶다면 predict_proba 메서드를 이용해볼 수도 있습니다.)

 

feature_names = cate_1hot_pipeline.named_steps["1hot"]
feature_names.categories_

cat_one_hot_features = []

i = 0

for feature in ['MultipleLines', 'InternetService', 'Contract', 'PaymentMethod']:
        
    cat_one_hot_features.append([feature+"_"+string for string in list(feature_names.categories_[i])])    

    i += 1
    
cat_one_hot_features = sum(cat_one_hot_features,[])
print(cat_one_hot_features)

features =  ["SupplementaryService", "Streaming"] + ["tenure","MonthlyCharges"] + cat_one_hot_features + ['SeniorCitizen', 'Partner', 'Dependents','PaperlessBilling']
print()
print(len(features))

feature_importances = pd.DataFrame(model_RF.fit(X_res,y_res).feature_importances_,index = features,columns=['importance']).sort_values('importance',ascending=False)
feature_importances = feature_importances.head(10)
feature_importances

sns.barplot(feature_importances["importance"],feature_importances.index)

중요한 변수

변수의 중요도를 확인한 결과 중요한 변수들은 Contract_Month-to-Month, tenure, InternetService_Fiber_optic 변수입니다. 

7. 솔루션 제시

론칭하기 전에 학습한 것, 한 일과하지 않은 일, 수립한 가정, 시스템 제한 사항 등을 강조하면서 설루션과 문서를 출시하고 깔끔한 도표와 기억하기 쉬운 제목으로 자료를 만들어놓습니다. 

 

8. 론칭, 모니터링, 시스템 보수

현재 모델은 론칭하기에 적절한 성능을 보인다고 하기 어렵습니다. 만약 훈련된 모델이 론칭 허가를 받았다면 시스템에 적용하기 위한 준비를 합니다. 코드를 정리하고 문서와 테스트 케이스를 작성합니다. 상용 환경에 배포하기 위해 전체 전처리 파이프라인과 예측 파이프라인이 포함된 모델을 joblib을 이용하여 저장하는 것이 좋습니다.

 

import joblib

joblib.dump(model_RF,"./model_RF.pkl")

 

배포한 후 일정 간격으로 시스템의 성능을 실시간으로 체크하고 성능이 떨어졌을 때 알람을 받을 수 있는 모니터링 코드를 작성해야 합니다. 데이터가 계속 변화하면 성능이 바뀔 수 있기 때문에 데이터셋을 업데이트하고 모델을 정기적으로 다시 훈련할 수 있도록 하이퍼파라미터를 자동으로 세부 튜닝하는 스크립트를 작성하는 것이 좋습니다. 또한 모델이 실패했을 때 무엇을 할지 정의하고 어떻게 대비할지 관련 프로세스를 모두 준비합니다. 

 

 

느낀점

이번 프로젝트는 평소에 하던 방법과 비슷하지만 좀 더 깊게 생각했던 프로젝트였습니다. 핸즈온 머신러닝을 참고했다는 부분이 컸던 것 같습니다. 데이터 분석 방법을 보다 체계적으로 다룬 느낌이랄까요? 이 포스팅에서는 임계값을 조절하여 정밀도와 재현율을 훈련시키는 부분은 작성되지 않았지만 핸즈온 Chapter2를 공부하며 알게되었습니다. 또한, 자동화의 중요성을 느끼고 파이프라인에 흥미를 가져 재밌는 시간이었습니다. 이제 데이터 분석을 할 때 파이프라인을 애용해보려고 합니다! AutoML을 처음 시도해봤었는데 데이터에 먼저 AutoML을 적용해보고 직접 분석하고 모델링하는 과정을 거쳐도 좋을 것 같다는 생각이 들었습니다. 정형데이터라면 AutoML에서 데이터 전처리를 자동으로 해주니까요! 직접 모델링 하기 전에 AutoML을 먼저 적용한 후 파생 변수를 만드는 것과 같은 특징공학, 더 좋은 모델을 선택하는 과정을 거칠 수 있도록 Pycaret 라이브러리에 대해 공부를 해봐야겠습니다.

 

참고

[1] 고객 이탈 예측 그림

[2] 핸즈온 머신러닝 2판 (도서)

[3] EDA 코드

[4] 모델링 코드