로지스틱 회귀분석

• 로지스틱(Logistic) 회귀분석은 회귀분석이라는 명칭과 달리 회귀분석 문제와 분류문제 모두에 사용할 수 있다. 로지스틱 회귀분석 모형에서는 종속변수가 이항분포를 따르고 그 모수 $ \theta $ 가 독립변수 $ x $ 에 의존한다고 가정한다.

$$p(y \mid x) = \text{Bin} (y; \theta (x), N)$$

• 위 식에서 보듯이 로지스틱 함수는 y의 값이 특정한 구간내의 값 (0 ~ N)만 가질 수 있기 때문에 종속변수가 이러한 특성을 가진 경우에 회귀분석 방법으로 사용할 수 있다.

• 또는 이항 분포의 특별한 경우 $ (N=1) $ 로 $ y $ 가 베르누이 확률분포인 경우도 있을 수 있다. 여기에서는 베르누이 확률분포를 따르는 로지스틱 회귀분석만 고려하기로 한다.

$$p(y \mid x) = \text{Bern} (y; \theta (x) )$$

• 종속변수 $ y $ 가 $ 0 $ 또는 $ 1 $ 인 분류 예측 문제를 풀 때는 $ x $ 값을 이용하여 $ \theta ( x ) $ 를 예측한 후 다음 기준에 따라 $ \hat{y} $ 값을 출력한다.

$$\hat{y} = \begin{cases} 1 & \text{ if } \theta (x) \geq 0.5 \\ 0 & \text{ if } \theta (x) < 0.5 \end{cases}$$

• 회귀분석을 할 때는 $ \hat{y} $ 으로 $ y = 1 $ 이 될 확률값 $ \theta ( x ) $ 을 직접 사용한다.

$$\hat{y} = \theta (x)$$

시그모이드 함수

• 로지스틱 회귀모형에서는 베르누이 확률분포의 모수 $ \theta $ 가 $ x $ 의 함수라고 가정한다. $ \theta (x) $ 는 $ x $ 에 대한 함수를 0부터 1 사이의 값만 나올 수 있도록 시그모이드 함수(sigmoid function)라는 함수를 사용하여 변형한 것을 사용한다.

• 시그모이드 함수는 종속변수의 모든 실수 값에 대해 유한한 구간 $ ( a, b ) $ 사이의 한정된 값을 가지고

$$a < f(x) < b$$

• 항상 양의 기울기를 가지는 단조증가하는

$$a > b \; \rightarrow \; f(a) > f(b)$$

• 함수의 집합을 말한다. 실제로는 다음과 같은 함수들이 주로 사용된다.

• 로지스틱(Logistic)함수

$$\text{logitstic}(z) = \sigma(z) = \dfrac{1}{1+\exp{(-z)}}$$

• 하이퍼볼릭 탄젠트(Hyperbolic tangent)함수

$$\tanh(z) = \frac{\sinh z}{\cosh z} = \frac{(e^z - e^{-z})/2}{(e^z + e^{-z})/2} = 2 \sigma(2z) - 1$$

• 오차(Error)함수

$$\text{erf}(z) = \frac{2}{\sqrt\pi}\int_0^z e^{-t^2}\,dt$$

• 하이퍼볼릭탄젠트 함수는 로지스틱함수를 위아래 방향으로 2배 늘리고 좌우 방향으로 1/2로 축소한 것과 같다.

---

xx = np.linspace(-5, 5, 1000)
plt.figure(figsize=(12,8))
plt.plot(xx, 1/(1+np.exp(-xx)), 'r-', label="로지스틱함수")
plt.plot(xx, sp.special.erf(0.5*np.sqrt(np.pi)*xx), 'g:', label="오차함수")
plt.plot(xx, np.tanh(xx), 'b--', label="하이퍼볼릭탄젠트함수")
plt.ylim([-1.1, 1.1])
plt.legend(loc=2)
plt.xlabel("x")
plt.savefig('sigmoid_functions_ploting')
plt.show()

---

로지스틱 함수

• 로지스틱함수는 음의 무한대부터 양의 무한대까지의 실수값을 0부터 1사이의 실수값으로 1대 1 대응시키는 시그모이드함수이다. 보통 시그모이드함수라고 하면 로지스틱함수를 의미한다.

• 위의 추정식에서 가장 오른쪽에 있는 로그 오즈비(log odds ratio)를 먼저 설명하자면, 로그 안에 들어간 오즈비(odds ratio)는 베르누이 시도에서 1이 나올 확률 $ \theta (x) $와 0이 나올 확률 $ 1 - \theta (x) $의 비율(ratio)을 의미한다.

$$odds ratio = \frac{ \theta (x) }{1 - \theta (x)}$$

0부터 1사이의 값만 가지는 확률값인 $ \theta (x) $ 를 오즈비로 변환하면 0부터 양의 무한대까지의 값을 가질 수 있다. 오즈비를 로그변환한 것이 위에서 언급한 Logit function이다. 이로써 로지트 함수의 값은 로그변환에 의해 음의 무한대 $(- \infty)$부터 양의 무한대$(\infty)$까지의 값을 가질 수 있다.

$$y = logit(odds ratio) = log ( \frac{ \theta (x) }{ 1 - \theta (x) } )$$

• 로지스틱함수(Logistic function)은 로지트 함수의 역함수이다. 즉 음의 무한대부터 양의 무한대 까지의 값을 가지는 입력변수를 0부터 1사이의 값을 가지는 출력변수로 변환한 것이다.

$$logistic(z) = \theta (z) = \frac{1}{1 + exp(-z)}$$

선형 판별 함수

• 또한, 로지스틱 회귀분석에서는 판별함수 수식으로 선형함수를 사용하므로 판별 경계면 또한 선형이 됨을 유의해야한다.

• 로지스틱함수 $ \sigma(z) $ 를 사용하는 경우에는 $ z $ 과 $ \theta $ 값은 다음과 같은 관계가 있다.

  • $ z = 0 $ 일 때 $ \theta = 0.5 $
  • $ z > 0 $ 일 때 $ \theta > 0.5; \rightarrow \hat{y} = 1 $
  • $ z < 0 $ 일 때 $ \theta < 0.5; \rightarrow \hat{y} = 0 $

• 즉, $ z $ 가 분류 모형의 판별 함수(decision function)의 역할을 한다. 로지스틱 회귀분석에서는 판별함수 수식으로 선형함수를 사용한다. 따라서 판별 경계면도 선형이 된다.

$$z = w^{T} x$$

• 위의 그림에서 알 수 있듯이, X의 범위가 [$-\infty$, $+\infty$]인 경우에 Y의 범위를 [0,1]로 만들어주는 함수이다.

• 위의 그림에서 좌측그림의 적합시킨 회귀선(파란선)을 보면 예측할 확률값이 500미만이면 음수를 갖을 수 있게 되는 것을 확인할 수 있다. 이는 확률값을 예측하는 모델을 만든다는 가정 자체에 위반하는 것이므로 이전에 학습했었던 회귀모형들과는 다른 방법의 회귀모형을 만들어 적합시켜야 할 것임이 분명해졌다.

• 일반적으로 decision boundary(threshold)를 0.5로 하지만 class가 imbalanced 한 경우는 조절하여보면서 적합시킬 수 있다.

• 기존의 일반 선형 회귀이듯, 다중 선형 회귀는 최소제곱법을 통해 해를 추정해 낼 수 있지만, 로지스틱 회귀는 비선형 방정식이 되어 동일하게 최소제곱법으로 회귀계수를 추정해내기 힘들다. 그러므로, 위의 그림처럼 MLE를 통해 회귀계수들을 추정한다.

그레디언트 벡터가 영벡터가 되는 모수의 값이 로그가능도를 최대화하는 값이다. 하지만 그레디언트 벡터 수식이 $ w $ 에 대한 비선형 함수이므로 선형모형과 같이 간단하게 그레디언트가 0이 되는 모수 $ w $ 값에 대한 수식을 구할 수 없으며 수치적인 최적화 방법(numerical optimization)을 통해 반복적으로 최적 모수 $ w $ 의 값을 구해야 한다.

수치적 최적화

• 위의 MLE 방식을 통해 parameter의 값을 업데이트하는데, 아래와 같은 방식으로 수치적 최적화를 진행한다. 위의 로그 가능도함수 $log(l(\beta_{0}, \beta_{1}))$을 편의상 LL이라하자. 로그가능도함수 LL을 최대화하는 것은 다음 목적함수를 최소화하는 것과 같다.

$$J = - LL$$

• SGD(Steepest Gradient Descent)방식을 사용(Stochastic Gradient Descent아님!!)하여 아래와 같은 그레디언 벡터를 구할 수 있다.

$$g_{k} = \frac{d}{dw} (-LL)$$

• 위에서 구한 그레디언트 방향으로 step size $(n_{k})$만큼 이동한다.

$$\begin{eqnarray} w_{k+1} &=& w_{k} - \eta_k g_k \\ &=& w_{k} + \eta_k \sum_{i=1}^N \big( y_i - \mu(x_i; w_k) \big) x_i\\ \end{eqnarray}$$

• 위의 그림에서 요약된 다중 로지스틱 회귀분석의 결과를 살펴보면, 맨 마지막 설명과 같이 해석할 수 있지만, 좀 더 자연스러운 해석하고 싶을 경우 아래 그림과 같이 Logit 보다 exp를 취해줘 Odds로 변환하여 해석하는 것을 권한다.

StatsModels 패키지의 로지스틱 회귀

• 다음과 같은 1차원 독립변수를 가지는 분류문제를 풀어보자.

---

from sklearn.datasets import make_classification

X0, y = make_classification(n_features=1, n_redundant=0, n_informative=1,
                            n_clusters_per_class=1, random_state=4)

plt.figure(figsize=(12,8))
plt.scatter(X0, y, c=y, s=100, edgecolor="k", linewidth=2)
sns.distplot(X0[y == 0, :], label="y = 0", hist=False)
sns.distplot(X0[y == 1, :], label="y = 1", hist=False)
plt.ylim(-0.2, 1.2)
# plt.savefig('logistic_regression_example_datas')
plt.show()

---

• StatsModels 패키지는 베르누이 분포를 따르는 로지스틱 회귀 모형 Logit를 제공한다. 사용방법은 OLS 클래스 사용법과 동일하다. 종속변수와 독립변수 데이터를 넣어 모형을 만들고 fit 메서드로 학습을 시킨다. fit 메서드의 disp=0 인수는 최적화 과정에서 문자열 메세지를 나타내지 않는 역할을 한다.

---

X = sm.add_constant(X0)
logit_mod = sm.Logit(y, X)
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())

결과

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                          Logit   Df Residuals:                       98
Method:                           MLE   Df Model:                            1
Date:                Mon, 08 Jun 2020   Pseudo R-squ.:                  0.7679
Time:                        00:12:48   Log-Likelihood:                -16.084
converged:                       True   LL-Null:                       -69.295
Covariance Type:            nonrobust   LLR p-value:                 5.963e-25
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.2515      0.477      0.527      0.598      -0.683       1.186
x1             4.2382      0.902      4.699      0.000       2.470       6.006
==============================================================================

---

• 결과 객체에서 summary메서드를 사용하여 리포트를 출력할 수 있다. 결과 리포트에서 판별함수의 수식이 다음과 같다는 것을 알 수 있다.

$$\theta (x) = \sigma(4.2382x + 0.2515)$$

• 따라서 $ z $ 값의 부로를 나누는 기준값은 $ 4.2382x + 0.2515 = 0.5 $ 가 되는 $ x $ 값 즉, $ (0.5-0.2515)/4.2382 $ 이다. predict 메서드를 사용하면 $ \theta (x) $ 값을 출력한다.

• 유의확률을 감안하면 상수항의 값은 0과 마찬가지이므로 $ \theta (x) $ 가 다음과 같다고 볼 수도 있다.

$$\mu(x) = \sigma(4.2382x)$$

• 이렇게 생각하면 $ z $ 값의 부호를 나누는 기준값은 실질적으로는 $ 0.5/4.2382 = 0.118 $ 이다. 실제 데이터가 네모 포인트로 되어있고, 원형의 데이터가 예측 데이터이다.

---

xx = np.linspace(-3, 3, 100)
mu = logit_res.predict(sm.add_constant(xx))

plt.figure(figsize=(12, 8))
plt.plot(xx, mu, lw=3)
plt.scatter(X0, y, c=y, s=100, edgecolor="k", lw=2)
plt.scatter(X0, logit_res.predict(X), label=r"$\hat{y}$", marker='s', c=y,
            s=100, edgecolor="k", lw=1)
plt.vlines(0.118, ymin=0, ymax=1, lw=1, linestyles='--', colors="red")
plt.xlim(-3, 3)
plt.xlabel("x")
plt.ylabel(r"$\mu$")
plt.title(r"$\hat{y} = \mu(x)$")
plt.legend()
# plt.savefig('predicted_value_y_hat_logistic_regression_example')
plt.show()

---

판별함수

• Logit 모형의 결과 객체에는 fittedvalues라는 속성으로 판별함수 $ z = w^{T}x $ 값이 들어가 있다. 이 값을 이용하여 분류문제를 풀 수도 있다.

---

plt.figure(figsize=(12,8))
plt.scatter(X0, y, c=y, s=100, edgecolor="k", lw=2, label="데이터")
plt.plot(X0, logit_res.fittedvalues, label="판별함수값")
plt.legend()
# plt.savefig('decision_function_with_logit_fitted_value')
plt.show()

---

• 위의 $ z $ 값을 통해 직접 시그모이드 함수(로지스틱함수)를 취하여 그려보면 다음과 같이 위에서 로지스틱함수를 그려 예측한 그림과 동일해 지는 것을 확인할 수 있다.

---

plt.figure(figsize=(12,8))
plt.scatter(X0, y, c=y, s=100, edgecolor="k", lw=2, label="데이터")
plt.scatter(X0, 1/(1+np.exp(-logit_res.fittedvalues)), label="로지스틱함수값")
plt.legend()
# plt.savefig('logistic_fuction_using_logit_fitted_value')
plt.show()

---

로지스틱 회귀 성능 측정

• 로지스틱 회귀 성능은 맥파든 의사결정계수(McFadden pseudo R square) 값으로 측정한다.

$$R^2_{\text{pseudo}} = 1 - \dfrac{G^2}{G^2_0}$$

• $ G^{2} $ 는 이탈도(deviance)라고 하는 양으로 다음과 같이 정의된다. 여기에서 $ \hat{y} $ 은 $ y=1 $ 일 확률 $ \theta $ 를 뜻한다.

$$G^2 = 2\sum_{i=1}^N \left( y_i\log\dfrac{y_i}{\hat{y}_i} + (1-y_i)\log\dfrac{1-y_i}{1-\hat{y}_i} \right)$$

• 이탈도는 모형이 100% 정확한 경우에는 0이 되고 모형의 성능이 나빠질수록 값이 커진다. 또한, 이탈도는 로그가능도에 음수를 취한 값과 같다.

$$G^2 = -LL$$

• $ G^{2} $ 는 현재 이탈도이고 $ G_{0}^{2} $ 는 귀무모형(null model)으로 측정한 이탈도이다.

• 귀무모형이란 모든 $ x $ 가 $ y $ 를 예측하는데 전혀 영향을 미치지 않는 모형을 말한다. 즉, 무조건부 확률 $ p(y) $ 에 따라 $ x $ 에 상관없이 동일하게 $ y $ 를 예측하는 모형을 말한다. 결국 우리가 만들 수 있는 가장 성능이 나쁜 모형이 된다.

$$\mu_{\text{null}} = \dfrac{\text{number of $Y=1$ data}}{\text{number of all data}}$$

• 따라서 맥파든 의사결정계수는 가장 성능이 좋을 때는 1이 되고 가장 성능이 나쁠 때는 0이된다.

• Scikit-learn 패키지의 metric 서브패키지에는 로그 손실을 계산하는 log_loss 함수가 있다. 위 예제에서 최적 모형의 로그 손실은 약 16.08로 계산된다.

---

from sklearn.metrics import log_loss

y_hat = logit_res.predict(X)
log_loss(y, y_hat, normalize=False)

결과

16.084355200413036

---

• 귀무 모형의 모수값을 구하면 0.51이고 이 값으로 로그 손실을 계산하면 약 69이다.

---

mu_null = np.sum(y) / len(y)
mu_null

결과

0.51

---

---

y_null = np.ones_like(y) * mu_null
log_loss(y, y_null, normalize=False)

결과

69.29471672244784

---

• 두 값을 이용하여 맥파든 의사 결정계수 값을 계산할 수 있다.

---

1 - (log_loss(y, y_hat) / log_loss(y, y_null))

결과

0.7678848264170398

---

Scikit-Learn 패키지의 로지스틱 회귀

• Scikit-Learn 패키지는 로지스틱 회귀 모형 LogisticRegression 를 제공한다.

---

from sklearn.linear_model import LogisticRegression

model_sk = LogisticRegression().fit(X0, y)

xx = np.linspace(-3, 3, 100)
mu = 1.0/(1 + np.exp(-model_sk.coef_[0][0]*xx - model_sk.intercept_[0]))
plt.plot(xx, mu)
plt.scatter(X0, y, c=y, s=100, edgecolor="k", lw=2)
plt.scatter(X0, model_sk.predict(X0), label=r"$\hat{y}$", marker='s', c=y,
            s=100, edgecolor="k", lw=1, alpha=0.5)
plt.xlim(-3, 3)
plt.xlabel("x")
plt.ylabel(r"$\mu$")
plt.title(r"$\hat{y}$ = sign $\mu(x)$")
plt.legend()
plt.show()

---

연습문제

1. 붓꽃 분류문제에서 클래스가 세토사와 베르시칼라 데이터만 사용하고 (setosa=0, versicolor=1) 독립변수로는 꽃받침 길이(Sepal Length)와 상수항만 사용하여 StatsModels 패키지의 로지스틱 회귀모형으로 결과를 예측하고 보고서를 출력한다. 이 보고서에서 어떤 값이 세토사와 베르시칼라를 구분하는 기준값(threshold)으로 사용되고 있는가?
2. 위 결과를 분류결과표(confusion matrix)와 분류결과보고서(classification report)로 나타내라.
3. 이 모형에 대해 ROC커브를 그리고 AUC를 구한다. 이 때 Scikit-Learn의 LogisticRegression을 사용하지 않고 위에서 StatsModels로 구한 모형을 사용한다.

• 1번에 대한 문제 풀이는 아래와 같다.

---

from sklearn.datasets import load_iris

iris=load_iris()
X=pd.DataFrame(iris.data[iris.target <= 1], columns=iris.feature_names)
X=X.loc[:, 'sepal length (cm)']
y=pd.Series(iris.target[iris.target <= 1], name="target")

import statsmodels.api as sm

X = sm.add_constant(X0)
logit_mod = sm.Logit(y, X)
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())

결과

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                 target   No. Observations:                  100
Model:                          Logit   Df Residuals:                       98
Method:                           MLE   Df Model:                            1
Date:                Mon, 08 Jun 2020   Pseudo R-squ.:                0.003818
Time:                        01:50:59   Log-Likelihood:                -69.050
converged:                       True   LL-Null:                       -69.315
Covariance Type:            nonrobust   LLR p-value:                    0.4669
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0025      0.201     -0.012      0.990      -0.396       0.391
x1             0.1287      0.177      0.726      0.468      -0.219       0.476
==============================================================================

---

• Setosa와 versicolor의 구분하는 기준값은 $ 1/1+exp(z)=0.5 $ 인 지점을 찾아야하므로 $ z= -0.0025+0.1287 x $ 이므로 $ x = 0.019425019425019424 $ 에 대한 수직선이 판별함수 경계선이 될 것이다.

---

from sklearn.metrics import classification_report, confusion_matrix

pred_y=np.array([1 if x >= 0.5 else 0 for x in list(logit_res.predict(X))])
y=np.array(y)

print(confusion_matrix(y, pred_y))
print()
print(classification_report(y, pred_y))

결과

[[27 23]
 [23 27]]

              precision    recall  f1-score   support

           0       0.54      0.54      0.54        50
           1       0.54      0.54      0.54        50

    accuracy                           0.54       100
   macro avg       0.54      0.54      0.54       100
weighted avg       0.54      0.54      0.54       100

---

---

from sklearn.metrics import roc_auc_score
print(roc_auc_score(y,logit_res.predict(X)))

결과

0.5327999999999999

---

---

from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y, logit_res.predict(X))
plt.figure(figsize=(12,8))
plt.plot(fpr, tpr)
plt.show()

---

로지스틱 회귀를 사용한 이진 분류의 예

• 다음 데이터는 미국 의대생의 입학관련 데이터이다. 데이터의 의미는 다음과 같다.

• Acceptance : 0이면 불합격, 1이면 합격

• BCPM : Bio/Chem/Physics/Math 과목의 학점 평균
• GPA : 전체과목 학점 평균
• VR : MCAT Verbal reasoning 과목 점수
• PS : MCAT Physical sciences 과목 점수
• WS : MCAT Writing sample 과목 점수
• BS : MCAT Biological science 과목 점수
• MCAT : MCAT 총점
• Apps : 의대 지원 횟수

---

data_med = sm.datasets.get_rdataset("MedGPA", package="Stat2Data")
df_med = data_med.data
df_med.tail()

---

• 일단 학점(GPA)과 합격여부의 관계를 살펴보자.

---

plt.figure(figsize=(12,8))
sns.stripplot(x="GPA", y="Acceptance", data=df_med,
              jitter=True, orient='h', order=[1, 0])
plt.grid(True)
# plt.savefig('MCAT_data_target')
plt.show()

---

• 로지스틱 회귀분석을 실시한다. MCAT = VR + PS + WS + BS 이므로 MCAT은 독립 변수에서 제외해야 한다. 데이터 행렬이 선형독립이어야 하기 때문이다.

---

model_med = sm.Logit.from_formula("Acceptance ~ Sex + BCPM + GPA + VR + PS + WS + BS + Apps", df_med)
result_med = model_med.fit()
print(result_med.summary())

결과

Optimization terminated successfully.
         Current function value: 0.280736
         Iterations 9
                           Logit Regression Results                           
==============================================================================
Dep. Variable:             Acceptance   No. Observations:                   54
Model:                          Logit   Df Residuals:                       45
Method:                           MLE   Df Model:                            8
Date:                Mon, 08 Jun 2020   Pseudo R-squ.:                  0.5913
Time:                        02:56:41   Log-Likelihood:                -15.160
converged:                       True   LL-Null:                       -37.096
Covariance Type:            nonrobust   LLR p-value:                 6.014e-07
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    -46.6414     15.600     -2.990      0.003     -77.216     -16.067
Sex[T.M]      -2.2835      1.429     -1.597      0.110      -5.085       0.518
BCPM          -6.1633      6.963     -0.885      0.376     -19.811       7.484
GPA           12.3973      8.611      1.440      0.150      -4.479      29.274
VR             0.0790      0.311      0.254      0.799      -0.530       0.688
PS             1.1673      0.539      2.164      0.030       0.110       2.225
WS            -0.7784      0.396     -1.968      0.049      -1.554      -0.003
BS             1.9184      0.682      2.814      0.005       0.582       3.255
Apps           0.0512      0.147      0.348      0.728      -0.237       0.340
==============================================================================

---

• 예측 결과와 실제 결과를 비교하면 다음과 같다.

---

df_med["Prediction"] = result_med.predict(df_med)
plt.figure(figsize=(12, 8))
sns.boxplot(x="Acceptance", y="Prediction", data=df_med)
# plt.savefig('box_plot_predict_value')
plt.show()

---

• 위 분석 결과를 토대로 유의하지 않은 변수들을 제외하고 PS와 BS 점수만을 이용하여 다시 회귀분석하면 다음과 같다.

---

model_med = sm.Logit.from_formula("Acceptance ~  PS + BS", df_med)
result_med = model_med.fit()
print(result_med.summary())

결과

Optimization terminated successfully.
         Current function value: 0.460609
         Iterations 7
                           Logit Regression Results                           
==============================================================================
Dep. Variable:             Acceptance   No. Observations:                   55
Model:                          Logit   Df Residuals:                       52
Method:                           MLE   Df Model:                            2
Date:                Mon, 08 Jun 2020   Pseudo R-squ.:                  0.3315
Time:                        03:04:00   Log-Likelihood:                -25.333
converged:                       True   LL-Null:                       -37.896
Covariance Type:            nonrobust   LLR p-value:                 3.503e-06
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    -15.5427      4.684     -3.318      0.001     -24.723      -6.362
PS             0.4798      0.316      1.518      0.129      -0.140       1.099
BS             1.1464      0.387      2.959      0.003       0.387       1.906
==============================================================================

---

• 위 결과를 바탕으로 다음 점수가 $ 15.5427 + 0.5 $ 보다 크면 합격이라고 예측할 수 있다.

$$0.4798 \text{PS} + 1.1464 \text{BS}$$

로지스틱 회귀를 사용한 회귀분석

• 로지스틱 회귀는 분류문제뿐만 아니라 종속변수 $ y $ 가 $ 0 $ 부터 $ 1 $ 까지 막혀있는 회귀분석 문제에도 사용할 수 있다. 이 때는 다음처럼 $ \theta $ 값을 종속변수 $ y $ 의 예측값으로 사용한다.

$$\hat{y} = \mu(x)$$

• 만약 실제 $ y $ 의 범위가 0 부터 1이 아니면 스케일링을 통해 바꿔야 한다.

예제
다음 데이터는 1974년도에 “여성은 가정을 보살피고 국가를 운영하는 일은 남자에게 맡겨두어야 한다”라는 주장에 대한 찬성, 반대 입장을 조사한 결과이다. 각 열은 다음을 뜻한다.

• education : 교육 기간
• sex : 성별
• agree : 찬성 인원
• disagree : 반대 인원
• ratio : 찬성 비율

---

data_wrole = sm.datasets.get_rdataset("womensrole", package="HSAUR")
df_wrole = data_wrole.data
df_wrole["ratio"] = df_wrole.agree / (df_wrole.agree + df_wrole.disagree)
df_wrole.tail()

---

• 교육을 많이 받는 사람일수록 찬성 비율이 감소하는 것을 볼 수 있다.

---

plt.figure(figsize=(12, 8))
sns.scatterplot(x="education", y="ratio", style="sex", data=df_wrole)
plt.grid(True)
plt.savefig('education_ratio_up')
plt.show()

---

• 분석 결과는 다음과 같다.

---

model_wrole = sm.Logit.from_formula("ratio ~ education + sex", df_wrole)
result_wrole = model_wrole.fit()
print(result_wrole.summary())

결과

Optimization terminated successfully.
         Current function value: 0.448292
         Iterations 6
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                  ratio   No. Observations:                   41
Model:                          Logit   Df Residuals:                       38
Method:                           MLE   Df Model:                            2
Date:                Mon, 08 Jun 2020   Pseudo R-squ.:                  0.3435
Time:                        03:13:11   Log-Likelihood:                -18.380
converged:                       True   LL-Null:                       -27.997
Covariance Type:            nonrobust   LLR p-value:                 6.660e-05
===============================================================================
                  coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------
Intercept       2.0442      0.889      2.299      0.022       0.302       3.787
sex[T.Male]    -0.1968      0.736     -0.267      0.789      -1.640       1.247
education      -0.2127      0.071     -2.987      0.003      -0.352      -0.073
===============================================================================

---

• 성별은 유의하지 않다는 것을 알게되었으므로 성별을 제외하고 다시 모형을 구한다.

---

model_wrole2 = sm.Logit.from_formula("ratio ~ education", df_wrole)
result_wrole2 = model_wrole2.fit()
print(result_wrole2.summary())

결과

Optimization terminated successfully.
         Current function value: 0.449186
         Iterations 6
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                  ratio   No. Observations:                   41
Model:                          Logit   Df Residuals:                       39
Method:                           MLE   Df Model:                            1
Date:                Mon, 08 Jun 2020   Pseudo R-squ.:                  0.3422
Time:                        03:14:51   Log-Likelihood:                -18.417
converged:                       True   LL-Null:                       -27.997
Covariance Type:            nonrobust   LLR p-value:                 1.202e-05
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.9345      0.781      2.478      0.013       0.405       3.464
education     -0.2117      0.071     -2.983      0.003      -0.351      -0.073
==============================================================================

---

---

plt.figure(figsize=(12, 8))
sns.scatterplot(x="education", y="ratio", data=df_wrole)
xx = np.linspace(0, 20, 100)
df_wrole_p = pd.DataFrame({"education": xx})
plt.plot(xx, result_wrole2.predict(df_wrole_p), "r-", lw=4, label="예측")
plt.legend()
plt.show()

---

실습

---

# 분석에 필요한 패키지 불러오기
import os
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn import metrics
from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score, roc_auc_score, roc_curve
import statsmodels.api as sm
import matplotlib.pyplot as plt
import itertools
import time

---

data info

• Experience : 경력
• Income : 수입
• Famliy : 가족단위
• CCAvg : 월 카드사용량
• Education : 교육수준 (1: undergrad; 2, Graduate; 3; Advance )
• Mortgage : 가계대출
• Securities : account 유가증권계좌유무
• CD account : 양도예금증서 계좌 유무
• Online : 온라인계좌유무
• CreidtCard : 신용카드유무

---

# Personal Loan 데이터 불러오기
ploan = pd.read_csv('../data/Personal Loan.csv')
ploan.head()

결과

---

# 의미없는 변수 제거 ID, zip code 제외
ploan_processed = ploan.dropna(axis=0).drop(['ID', 'ZIP Code'], axis=1, inplace=False)

# 상수항 추가
ploan_processed = sm.add_constant(ploan_processed, has_constant="add")
ploan_processed.head()

결과

설명변수(X), 타켓변수(Y) 분리 및 학습데이터와 평가데이터

---

# 대출여부: 1 or 0
feature_columns = ploan_processed.columns.difference(['Personal Loan'])
X = ploan_processed[feature_columns]
y = ploan_processed['Personal Loan']

train_x, test_x, train_y, test_y = train_test_split(X, y, stratify=y,train_size=0.7,test_size=0.3,random_state=42)
print(train_x.shape, test_x.shape, train_y.shape, test_y.shape)

결과

(1750, 12) (750, 12) (1750,) (750,)

---

로지스틱회귀모형 모델링 y = f(x)

---

## 로지스틱 모형 적합
model = sm.Logit(train_y, train_x)
results = model.fit(method='newton')

결과

Optimization terminated successfully.
         Current function value: 0.131055
         Iterations 9

---

---

print(results.summary())

결과

Logit Regression Results                           
==============================================================================
Dep. Variable:          Personal Loan   No. Observations:                 1750
Model:                          Logit   Df Residuals:                     1738
Method:                           MLE   Df Model:                           11
Date:                Sat, 13 Jun 2020   Pseudo R-squ.:                  0.6030
Time:                        15:20:24   Log-Likelihood:                -229.35
converged:                       True   LL-Null:                       -577.63
Covariance Type:            nonrobust   LLR p-value:                2.927e-142
======================================================================================
coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------
Age                    0.0245      0.102      0.240      0.810      -0.175       0.224
CCAvg                  0.0985      0.063      1.562      0.118      -0.025       0.222
CD Account             4.3726      0.568      7.703      0.000       3.260       5.485
CreditCard            -1.2374      0.337     -3.667      0.000      -1.899      -0.576
Education              1.5203      0.190      7.999      0.000       1.148       1.893
Experience            -0.0070      0.102     -0.069      0.945      -0.206       0.192
Family                 0.7579      0.128      5.914      0.000       0.507       1.009
Income                 0.0547      0.004     12.659      0.000       0.046       0.063
Mortgage              -0.0001      0.001     -0.144      0.885      -0.002       0.002
Online                -0.4407      0.263     -1.674      0.094      -0.957       0.075
Securities Account    -1.8520      0.561     -3.299      0.001      -2.952      -0.752
const                -13.9203      2.773     -5.021      0.000     -19.354      -8.486
======================================================================================

---

---

#회귀계수 출력
results.params

결과

Age                    0.024471
CCAvg                  0.098468
CD Account             4.372577
CreditCard            -1.237447
Education              1.520329
Experience            -0.007032
Family                 0.757911
Income                 0.054695
Mortgage              -0.000133
Online                -0.440746
Securities Account    -1.852006
const                -13.920298
dtype: float64

---

회귀 계수 해석

• 나이가 한살 많을수록록 대출할 확률이 1.024 높다.
• 수입이 1단위 높을소룩 대출할 확률이 1.05배 높다
• 가족 구성원수가 1많을수록 대출할 확률이 2.13배 높다
• 경력이 1단위 높을수록 대출할 확률이 0.99배 높다(귀무가설 채택)
• Experience, Mortgage는 제외할 필요성이 있어보임

---

np.exp(results.params)

결과

Age                   1.024773e+00
CCAvg                 1.103479e+00
CD Account            7.924761e+01
CreditCard            2.901239e-01
Education             4.573729e+00
Experience            9.929928e-01
Family                2.133814e+00
Income                1.056218e+00
Mortgage              9.998665e-01
Online                6.435563e-01
Securities Account    1.569221e-01
const                 9.005163e-07
dtype: float64

---

---

## y_hat 예측
pred_y = results.predict(test_x)

def cut_off(y,threshold):
    Y = y.copy() # copy함수를 사용하여 이전의 y값이 변화지 않게 함
    Y[Y>threshold]=1
    Y[Y<=threshold]=0
    return(Y.astype(int))

pred_Y = cut_off(pred_y,0.5)

# confusion matrix
cfmat = confusion_matrix(test_y, pred_Y)
print(cfmat)

결과

[[661  12]
 [ 28  49]]

---

---

## confusion matrix accuracy계산하기
accuracy = (cfmat[0,0]+cfmat[1,1])/cfmat.sum()
accuracy

def acc(cfmat):
    acc = (cfmat[0,0]+cfmat[1,1])/cfmat.sum()
    return acc

결과

0.9466666666666667

---

임계값(cut-off)에 따른 성능지표 비교

• threshold에 따라 정확도가 달라지므로 가장 좋은 성능을 보여주는 0.5~0.6사이의 threshold를 사용하는 것이 좋을 것이다.

---

threshold = np.arange(0,1,0.1)
table = pd.DataFrame(columns=['ACC'])
for i in threshold:
    pred_Y = cut_off(pred_y,i)
    cfmat = confusion_matrix(test_y, pred_Y)
    table.loc[i] = acc(cfmat)
table.index.name='threshold'
table.columns.name='performance'
table

결과

---

# sklearn ROC 패키지 제공
fpr, tpr, thresholds = metrics.roc_curve(test_y, pred_y, pos_label=1)

# Print ROC curve
plt.plot(fpr,tpr)

# Print AUC
auc = np.trapz(tpr,fpr)
print('AUC:', auc)

plt.show()

결과

변수선택법

---

feature_columns = list(ploan_processed.columns.difference(["Personal Loan"]))
X = ploan_processed[feature_columns]
y = ploan_processed['Personal Loan'] # 대출여부: 1 or 0

train_x, test_x, train_y, test_y = train_test_split(X, y, stratify=y,train_size=0.7,test_size=0.3,random_state=42)
print(train_x.shape, test_x.shape, train_y.shape, test_y.shape)

결과

(1750, 12) (750, 12) (1750,) (750,)

---

---

def processSubset(X,y, feature_set):
            model = sm.Logit(y,X[list(feature_set)])
            regr = model.fit()
            AIC = regr.aic
            return {"model":regr, "AIC":AIC}

'''
전진선택법
'''
def forward(X, y, predictors):
    # 데이터 변수들이 미리정의된 predictors에 있는지 없는지 확인 및 분류
    remaining_predictors = [p for p in X.columns.difference(['const']) if p not in predictors]
    tic = time.time()
    results = []
    for p in remaining_predictors:
        results.append(processSubset(X=X, y= y, feature_set=predictors+[p]+['const']))
    # 데이터프레임으로 변환
    models = pd.DataFrame(results)

    # AIC가 가장 낮은 것을 선택
    best_model = models.loc[models['AIC'].argmin()] # index
    toc = time.time()
    print("Processed ", models.shape[0], "models on", len(predictors)+1, "predictors in", (toc-tic))
    print('Selected predictors:',best_model['model'].model.exog_names,' AIC:',best_model[0] )
    return best_model

def forward_model(X,y):
    Fmodels = pd.DataFrame(columns=["AIC", "model"])
    tic = time.time()
    # 미리 정의된 데이터 변수
    predictors = []
    # 변수 1~10개 : 0~9 -> 1~10
    for i in range(1, len(X.columns.difference(['const'])) + 1):
        Forward_result = forward(X=X,y=y,predictors=predictors)
        if i > 1:
            if Forward_result['AIC'] > Fmodel_before:
                break
        Fmodels.loc[i] = Forward_result
        predictors = Fmodels.loc[i]["model"].model.exog_names
        Fmodel_before = Fmodels.loc[i]["AIC"]
        predictors = [ k for k in predictors if k != 'const']
    toc = time.time()
    print("Total elapsed time:", (toc - tic), "seconds.")

    return(Fmodels['model'][len(Fmodels['model'])])


'''
후진소거법
'''
def backward(X,y,predictors):
    tic = time.time()
    results = []

    # 데이터 변수들이 미리정의된 predictors 조합 확인
    for combo in itertools.combinations(predictors, len(predictors) - 1):
        results.append(processSubset(X=X, y= y,feature_set=list(combo)+['const']))
    models = pd.DataFrame(results)

    # 가장 낮은 AIC를 가진 모델을 선택
    best_model = models.loc[models['AIC'].argmin()]
    toc = time.time()
    print("Processed ", models.shape[0], "models on", len(predictors) - 1, "predictors in",
          (toc - tic))
    print('Selected predictors:',best_model['model'].model.exog_names,' AIC:',best_model[0] )
    return best_model


def backward_model(X, y):
    Bmodels = pd.DataFrame(columns=["AIC", "model"], index = range(1,len(X.columns)))
    tic = time.time()
    predictors = X.columns.difference(['const'])
    Bmodel_before = processSubset(X,y,predictors)['AIC']
    while (len(predictors) > 1):
        Backward_result = backward(X=train_x, y= train_y, predictors = predictors)
        if Backward_result['AIC'] > Bmodel_before:
            break
        Bmodels.loc[len(predictors) - 1] = Backward_result
        predictors = Bmodels.loc[len(predictors) - 1]["model"].model.exog_names
        Bmodel_before = Backward_result['AIC']
        predictors = [ k for k in predictors if k != 'const']

    toc = time.time()
    print("Total elapsed time:", (toc - tic), "seconds.")
    return (Bmodels['model'].dropna().iloc[0])


'''
단계적 선택법
'''
def Stepwise_model(X,y):
    Stepmodels = pd.DataFrame(columns=["AIC", "model"])
    tic = time.time()
    predictors = []
    Smodel_before = processSubset(X,y,predictors+['const'])['AIC']
    # 변수 1~10개 : 0~9 -> 1~10
    for i in range(1, len(X.columns.difference(['const'])) + 1):
        Forward_result = forward(X=X, y=y, predictors=predictors) # constant added
        print('forward')
        Stepmodels.loc[i] = Forward_result
        predictors = Stepmodels.loc[i]["model"].model.exog_names
        predictors = [ k for k in predictors if k != 'const']
        Backward_result = backward(X=X, y=y, predictors=predictors)
        if Backward_result['AIC']< Forward_result['AIC']:
            Stepmodels.loc[i] = Backward_result
            predictors = Stepmodels.loc[i]["model"].model.exog_names
            Smodel_before = Stepmodels.loc[i]["AIC"]
            predictors = [ k for k in predictors if k != 'const']
            print('backward')
        if Stepmodels.loc[i]['AIC']> Smodel_before:
            break
        else:
            Smodel_before = Stepmodels.loc[i]["AIC"]
    toc = time.time()
    print("Total elapsed time:", (toc - tic), "seconds.")
    return (Stepmodels['model'][len(Stepmodels['model'])])

---

---

Forward_best_model = forward_model(X=train_x, y= train_y)
Stepwise_best_model = Stepwise_model(X=train_x,y=train_y)
Backward_best_model = backward_model(X=train_x,y=train_y)

pred_y_full = results.predict(test_x) # full model
pred_y_forward = Forward_best_model.predict(test_x[Forward_best_model.model.exog_names])
pred_y_backward = Backward_best_model.predict(test_x[Backward_best_model.model.exog_names])
pred_y_stepwise = Stepwise_best_model.predict(test_x[Stepwise_best_model.model.exog_names])

pred_Y_full= cut_off(pred_y_full,0.5)
pred_Y_forward = cut_off(pred_y_forward,0.5)
pred_Y_backward = cut_off(pred_y_backward,0.5)
pred_Y_stepwise = cut_off(pred_y_stepwise,0.5)

cfmat_full = confusion_matrix(test_y, pred_Y_full)
cfmat_forward = confusion_matrix(test_y, pred_Y_forward)
cfmat_backward = confusion_matrix(test_y, pred_Y_backward)
cfmat_stepwise = confusion_matrix(test_y, pred_Y_stepwise)

print(acc(cfmat_full))
print(acc(cfmat_forward))
print(acc(cfmat_backward))
print(acc(cfmat_stepwise))

결과

0.9466666666666667
0.944
0.944
0.944

---

• 변수선택법을 통한 성능이나 full model이나 성능차이가 거의 없으므로 full model을 사용할 것이다.