다운로드
작성자: admin 작성일시: 2016-06-10 16:25:07 조회수: 15621 다운로드: 670
카테고리: 머신 러닝 태그목록:

6.1 로지스틱 회귀분석

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

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

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

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

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

종속변수 $y$가 0또는 1인 분류 예측 문제를 풀 때는 $x$ 값에 따라 $\mu(x)$를 예측한 후 다음 기준에 따라 y를 예측한다.

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

또는 $\hat{y}$로 $y=1$이 될 확률값 $\mu(x)$를 직접 출력할 수도 있다.

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

시그모이드함수

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

$$ \mu = f(w^Tx) $$

모수 $\mu$는 일반적인 회귀분석의 종속변수와 달리 0 부터 1까지의 실수값만 가질 수 있기 때문에 시그모이드함수이라 불리는 특별한 형태의 함수 $f$를 사용해야 한다. 시그모이드함수는 종속변수의 모든 실수 값에 대해

  • 유한한 구간 $(a,b)$ 사이의 한정된(bounded) 값과
  • 항상 양의 기울기를 가지는 (단조증가)

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

  • 로지스틱(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 $$
In [1]:
xx = np.linspace(-5, 5, 1000)
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.show()

로지스틱함수

로지스틱함수는 무한대의 실수값을 0부터 1사이의 실수값을 1대1 대응시키는 시그모이드함수다. 다음 과정을 통해 정의되었다.

베르누이 시도에서 1이 나올 확률 $\mu$ 과 0이 나올 확률 $1-\mu$ 의 비(ratio)는 다음과 같은 수식이 되며 이를 승산비(odds ratio)라고 한다.

$$ \text{odds ratio} = \dfrac{\mu}{1-\mu} $$

0부터 1사이의 값만 가지는 $\mu$를 승산비로 변환하면 0부터 $\infty$의 값을 가질 수 있다.

승산비를 로그 변환한 것이 로지트함수(Logit function)다.

$$ z = \text{logit}(\text{odds ratio}) = \log \left(\dfrac{\mu}{1-\mu}\right) $$

로지트함수의 값은 로그 변환에 의해 $-\infty$부터 $\infty$까지의 값을 가질 수 있다.

로지스틱함수(Logistic function)는 로지트함수의 역함수이다. 즉 $-\infty$부터 $\infty$까지의 값을 가지는 변수를 0부터 1사의 값으로 변환한 결과이다.

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

선형 판별함수

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

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

즉 $z$가 분류 모형의 판별함수(decision function)의 역할을 한다.

$$ z = w^Tx $$

에서 로지스틱 회귀모형의 영역 경계면은 선형이라는 것을 알 수 있다.

로지스틱 회귀분석 모형의 모수 추정

로지스틱 회귀분석 모형은 일종의 비선형 회귀모형이지만 다음과 같이 최대가능도(Maximum Likelihood Estimation, MLE)방법으로 모수 $w$를 추정할 수 있다.

$$ p(y \mid x) = \text{Bern} (y;\mu(x;w) ) = \mu(x;w)^y ( 1 - \mu(x;w) )^{1-y} $$

$\mu$를 로지스틱함수 형태로 표현하면 다음과 같다.

$$ \mu(x;w) = \dfrac{1}{1 + \exp{(-w^Tx)}} $$

즉,

$$ \begin{eqnarray} p(y \mid x) &=& \left( \dfrac{1}{1 + \exp{(-w^Tx)}} \right) ^y \left( 1 - \dfrac{1}{1 + \exp{(-w^Tx)}} \right) ^{1-y} \\ &=& \left( \dfrac{1}{1 + \exp{(-w^Tx)}} \right) ^y \left( \dfrac{\exp{(-w^Tx)}}{1 + \exp{(-w^Tx)}} \right) ^{1-y} \\ \end{eqnarray} $$

데이터 표본이 $\{ x_i, y_i \}_{1:N}$일 경우 로그가능도 $\text{LL}$ 를 구하면 다음과 같다.

베르누이 확률분포의 정의에서

$$ \begin{eqnarray} \text{LL} &=& \log \prod_{i=1}^N \mu(x_i;w)^{y_i} (1-\mu(x_i;w))^{1-y_i} \\ &=& \sum_{i=1}^N \left( y_i \log\mu(x_i;w) + (1-y_i)\log(1-\mu(x_i;w)) \right) \\ &=& \sum_{i=1}^N \left( y_i \log\left(\dfrac{1}{1 + \exp{(-w^Tx_i)}}\right) + (1-y_i)\log\left(\dfrac{\exp{(-w^Tx_i)}}{1 + \exp{(-w^Tx_i)}}\right) \right) \\ \end{eqnarray} $$

가 된다.

로그가능도를 최대화하는 $w$ 값을 구하기 위해 다음처럼 미분을 한다.

$$ \dfrac{\partial\text{LL}}{\partial w} = \sum_{i=1}^N \dfrac{\partial\text{LL}}{\partial \mu(x_i;w)} \dfrac{\partial\mu(x_i;w)}{\partial w} $$

우선 $\mu$를 $w$로 미분하면

$$ \dfrac{\partial \mu(x_i;w)}{\partial w} = \dfrac{\partial}{\partial w} \dfrac{1}{1 + \exp{(-w^Tx_i)}} \ = \dfrac{\exp{(-w^Tx_i)}}{(1 + \exp{(-w^Tx_i)})^2} x_i \ = \mu(x_i;w)(1-\mu(x_i;w)) x_i $$

LL을 $\mu$로 미분하면

$$ \dfrac{\partial\text{LL}}{\partial \mu(x_i;w)} = \left( y_i \dfrac{1}{\mu(x_i;w)} - (1-y_i)\dfrac{1}{1-\mu(x_i;w)} \right) $$

두 식을 곱하면

$$ \begin{eqnarray} \dfrac{\partial \text{LL}}{\partial w} &=& \sum_{i=1}^N \left( y_i \dfrac{1}{\mu(x_i;w)} - (1-y_i)\dfrac{1}{1-\mu(x_i;w)} \right) \mu(x_i;w)(1-\mu(x_i;w)) x_i \\ &=& \sum_{i=1}^N \big( y_i (1-\mu(x_i;w)) - (1-y_i)\mu(x_i;w) \big) x_i \\ &=& \sum_{i=1}^N \big( y_i - \mu(x_i;w) \big) x_i \\ \end{eqnarray} $$

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

수치적 최적화

로그 가능도함수 $LL$을 최대화하는 것은 다음 목적함수를 최소화하는 것과 같다.

$$ J = -LL $$

최대경사도(Steepest Gradient Descent)방법을 사용하자.

그레디언트 벡터는 $$ g_k = \dfrac{d}{dw}(-LL) $$

이고, 이 방향으로 스텝사이즈 $\eta_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} $$

정규화

로지스틱 회귀에서도 과최적화를 방지하기 위해 릿지, 라소, 일레스틱넷 방식의 정규화 페널티를 목적 함수인 로그가능도함수에 추가할 수 있다. 예를 들어 릿지 정규화를 하면 다음과 같은 목적 함수를 최소화하는 것과 같다.

$$ J = -\text{LL} + \lambda \| w \|^2 $$

StatsModels 패키지의 로지스틱 회귀

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

In [2]:
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)
X = sm.add_constant(X0)

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.show()

StatsModels 패키지는 베르누이 분포를 따르는 로지스틱 회귀 모형 Logit 를 제공한다. 사용방법은 OLS 와 동일하다.

In [3]:
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:                Tue, 23 Jul 2019   Pseudo R-squ.:                  0.7679
Time:                        20:28:11   Log-Likelihood:                -16.084
converged:                       True   LL-Null:                       -69.295
                                        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
==============================================================================
In [4]:
xx = np.linspace(-3, 3, 100)
mu = logit_res.predict(sm.add_constant(xx))
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.xlim(-3, 3)
plt.xlabel("x")
plt.ylabel(r"$\mu$")
plt.title(r"$\hat{y} = \mu(x)$")
plt.legend()
plt.show()

판별함수

Logit 모형의 결과 객체에는 fittedvalues라는 속성으로 판별함수 $z=w^Tx$ 값이 들어가 있다.

In [5]:
plt.scatter(X0, y, c=y, s=100, edgecolor="k", lw=2, label="데이터")
plt.plot(X0, logit_res.fittedvalues * 0.1, label="판별함수값")
plt.legend()
plt.show()

성능 측정

로지스틱 회귀 성능은 McFadden pseudo R square 값으로 측정한다.

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

$G^2$는 deviance 라고 하는 양으로 다음과 같이 정의된다. 로그 손실(log loss)이라고도 한다.

$$ 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) $$

여기에서 $\hat{y}$는 다음처럼 $y=1$일 확률을 뜻한다.

$$ \hat{y}_i = \mu(x_i) $$

deviance는 모형이 100% 정확한 경우에는 0이 되고 모형의 성능이 나빠질 수록 값이 커진다.

이 값은 로그 가능도의 음수값과 같다.

$$ G^2 = - LL $$

$G^2$는 현재 deviance이고 $G^2_0$는 귀무모형(null model)으로 측정한 deviance이다.

귀무무형이란 모든 $x$가 $y$를 예측하는데 전혀 영향을 미치지 않는 모형을 말한다. 즉, 무조건부 확률 $p(y)$에 따라 $x$에 상관없이 동일하게 $y$를 예측하는 모형을 말한다.

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

scikit-learn 패키지의 metric 서브패키지에는 로그 손실을 계산하는 log_loss 함수가 있다. normalize=False로 놓으면 위와 같은 값을 구한다. normalize 인수의 디폴트 값은 True이고 이 때는 로그 손실의 평균값을 출력한다.

위 예제에서 최적 모형의 로그 손실은 약 16이다.

In [6]:
from sklearn.metrics import log_loss

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

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

In [7]:
mu_null = np.sum(y) / len(y)
mu_null
Out:
0.51
In [8]:
y_null = np.ones_like(y) * mu_null
log_loss(y, y_null, normalize=False)
Out:
69.29471672244784

두 값을 이용하여 McFadden pseudo R square 값을 계산할 수 있다.

In [9]:
1 - (log_loss(y, y_hat) / log_loss(y, y_null))
Out:
0.7678848264170398

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

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

In [10]:
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()

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

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

  • Acceptance: 0이면 불합격, 1이면 합격
  • BCPM: Bio/Chem/Physics/Math 과목의 학점 평균
  • GPA: 전체과목 학점 평균
  • VR: MCAT Verbal reasoning 과목 점수
  • PS: MCAT Physical sciences 과목 점수
  • WS: MCAT Writing sample 과목 점수
  • BS: MCAT Biological sciences 과목 점수
  • MCAT: MCAT 촘점
  • Apps: 의대 지원 횟수
In [11]:
data_med = sm.datasets.get_rdataset("MedGPA", package="Stat2Data")
df_med = data_med.data
df_med.tail()
Out:
Accept Acceptance Sex BCPM GPA VR PS WS BS MCAT Apps
50 D 0 M 2.41 2.72 8 8 8.0 8 32 7
51 D 0 M 3.51 3.56 11 8 6.0 9 34 6
52 A 1 F 3.43 3.48 7 10 7.0 10 34 14
53 D 0 M 2.61 2.80 7 5 NaN 6 18 6
54 D 0 M 3.36 3.44 11 11 8.0 9 39 1

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

In [12]:
sns.stripplot(x="GPA", y="Acceptance", data=df_med,
              jitter=True, orient='h', order=[1, 0])
plt.grid(True)
plt.show()

로지스틱 회귀분석을 실시한다. MCAT = VR + PS + WS + BS이므로 이 MCAT은 독립 변수에서 제외해야 한다.

In [13]:
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:                Tue, 23 Jul 2019   Pseudo R-squ.:                  0.5913
Time:                        20:28:16   Log-Likelihood:                -15.160
converged:                       True   LL-Null:                       -37.096
                                        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
==============================================================================
In [14]:
df_med["Prediction"] = result_med.predict(df_med)
sns.boxplot(x="Acceptance", y="Prediction", data=df_med)
plt.show()

연습 문제 1

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

연습 문제 2

  1. 붓꽃 분류문제에서 클래스가 베르시칼라(versicolor)와 버지니카(virginica) 데이터만 사용하여(versicolor=0, virginica=1) 로지스틱 회귀모형으로 결과를 예측하고 보고서를 출력한다. 독립변수는 모두 사용한다. 이 보고서에서 세토사와 베르시칼라를 구분하는 경계면의 방정식을 찾아라.
  2. 위 결과를 분류결과표와 분류결과보고서로 나타내라.
  3. 이 모형에 대해 ROC커브를 그리고 AUC를 구하라. 이 때 Scikit-Learn의 LogisticRegression을 사용하지 않고 위에서 StatsModels로 구한 모형을 사용한다.

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

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

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

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

예제

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

  • education: 교육 기간
  • sex: 성별
  • agree: 찬성 인원
  • disagree: 반대 인원
  • ratio: 찬성 비율
In [15]:
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()
Out:
education sex agree disagree ratio
37 16 Female 13 115 0.101562
38 17 Female 3 28 0.096774
39 18 Female 0 21 0.000000
40 19 Female 1 2 0.333333
41 20 Female 2 4 0.333333

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

In [16]:
sns.scatterplot(x="education", y="ratio", style="sex", data=df_wrole)
plt.grid(True)
plt.show()