9  로지스틱 회귀분석 (logistic regression)

Acknowledgement: 본 절의 구성에는 다음 교재를 발췌하였다.

주 출처

보조 출처

9.1 반응변수가 0/1인 경우

  • 질적 변수
  • 입시, 입사 등 선발과정: ‘잘함(good)’ vs ‘못함(poor)’
  • 건강 연구: 반응변수는 암 발병 여부(yes vs. no), 설명변수는 나이, 성별, 흡연, 다이어트, 가족의 의료기록 등
  • 금융: 반응변수는 회사의 지불능력(파산 vs. 지불능력 있음), 예측변수들은 회사와 관련된 다양한 금융 정보들
  • \(P(Y=1 \mid {\bf X}={\bf x}) =: p({\bf x})\)에 대한 확률 모형이 필요함.

9.2 로지스틱 회귀분석의 확률모형

9.2.1 로지스틱 선형모형

  • \(Y \in\{0,1\}\) 및 \({\bf X} \in \mathbb{R}^{p}\) 가정하고, \(p({\bf x})= P(Y=1 \mid {\bf X}={\bf x})\) 라 쓰면,

\[ \log \left(\frac{p({\bf x})}{1-p({\bf x})}\right)=\beta_{0}+{\bf x}^{T} \boldsymbol{\beta} \tag{1} \] \[ \text { 또는 } p({\bf x})=\frac{\exp \left(\beta_{0}+ {\bf x}^{T} \boldsymbol{\beta}\right)}{1+\exp \left(\beta_{0}+{\bf x}^{T} \boldsymbol{\beta}\right)} \]

  • 로짓함수 (logit function): \(\operatorname{logit}(t):=\log \{t /(1-t)\}\)
  • \(0<t<1,-\infty<\operatorname{logit}(t)<\infty, t\) 의 단조증가함수
  • \(\operatorname{logit}^{-1}(t)=1 /\{1+\exp (-t)\}=\exp (t) /\{1+\exp (t)\}\)

9.2.2 웬 비선형 변환?

  • 먼저 \(Y=0,1\)인 경우 \(E(Y) = P(Y=1)\)임을 기억하자. 마찬가지로 \(E(Y|{\bf X}={\bf x}) = P(Y=1|{\bf X}={\bf x}) = p({\bf x})\)이다.
  • 따라서, 지난 7,8장에서 다룬 선형회귀분석 모형을 적용하자면 \(p({\bf x}) = \beta_0 + {\bf x}^{T} \boldsymbol{\beta}\) 식의 모형을 고려할 수 있고, 이 모형에서의 계수추정은 이미 7,8장에서 다룬 바 있다. 이 모형을 이진반응변수에 적합하면 안 되는 걸까?
  • 결론적으로 “완벽히 안 된다”라고 말할 수는 없지만, 여러 가지 문제가 발생한다. 가장 큰 문제는 조건부확률의 적합값에서 발생한다. \(P(Y=1|{\bf X}={\bf x}) = \beta_0 + {\bf x}^{T} \boldsymbol{\beta}\)의 좌변은 확률이니 \({\bf x}\)의 값에 상관없이 늘 0과 1사이에 있어야 한다. 그런데 \(\beta_0 + {\bf x}^{T} \boldsymbol{\beta}\) 항은 \({\bf x}\)의 값에 따라 0보다 작을 수도 있고 1보다 클 수도 있다. 즉, 선형회귀분석 모형에서는 \(p({\bf x})\)의 적합값이 0보다 작거나 1보다 클 수도 있다는 문제점이 있다. (아래 그림의 첫번째)
  • \(p({\bf x})\)의 적합값이 반드시 0과 1사이에 오게 하는 방법이 있다. 어떤 알려진 변환 \(g\)에 대하여 \(g( p({\bf x}) )=\beta_0 + {\bf x}^{T} \boldsymbol{\beta}\) 라고 가정하는 것이다. 이때 \(g\) 함수는 \(p_{\mathbf{x}}\) 를 선형함수인 \(\beta_0 + \mathbf{x}^T \boldsymbol{\beta}\) 와 연결시키는 함수로서 ‘연결함수’(link function)라 한다.
  • 로지스틱 회귀분석모형의 link function은 logit function \(\log \{t /(1-t)\}\) 이다.
  • 다른 대표적 link function은 probit function 으로, \(t \mapsto \int_{- \infty}^t \phi (u) du\)의 역함수로 정의된다 (단, \(\phi(u)\)는 표준정규분포의 확률밀도함수). 계량경제학에서 가격에 따른 선택을 연구할 때 많이 사용한다.
  • 로짓 함수와 프로빗 함수 모두 정의역은 \((0,1)\)이고 치역은 모든 실수이다.
  • 본 장에서는 logit link function을 중점적으로 다룬다.

FIGURE 4.2 in [ISLR]. Left: 파란색 선은 선형회귀분석(즉 최소제곱법)으로 추정한 확률값이다. 추정된 확률값이 음수도 된다! 오렌지색 점들은 자료에서의 \(y\)값 (\(0 / 1\)) 들이다. Right: 로지스틱 회귀분석으로 추정한 확률값이다. 모든 확률값들이 0과 1사이에 있다.

9.2.3 로지스틱 선형모형에서 계수의 의미

  • 오즈(odds) (또는 \(y=1\) 의 오즈): 클래스 0 에 속할 확률에 대한 클래스 1 에 속할 확률의 비율로 정의된다. 어떤 질병이 걸릴 위험을 묘사할 때 절대위험차이(absolute risk difference), 상대위험비(relative risk ratio)와 함께 대표적으로 사용되는 지표이다.

\[ \begin{aligned} \text { odds }: & =\frac{p}{1-p} \\ \operatorname{odds}(Y=1 \mid X & = {\bf x}):=\frac{p({\bf x})}{1-p({\bf x})} \end{aligned} \]

  • 로지스틱 다중회귀모형에서 \(\beta_j\)의 해석: (1)에서 \(X_{j}\) 가 한 단위 증가하면 \(y=1\) 의 오즈는 \(\exp \left(\beta_{j}\right)\) 배만큼 증가한다 (다른 독립변수들이 모두 일정하다는 가정 하에). 즉 \(\beta\)값을 \(Y\)의 평균적 증분으로 해석할 수 있는 선형회귀분석(7,8장)보다는 조금 우회적으로 해석이 된다.

  • 요컨대,

    • 로지스틱 단순선형회귀모형에서 \(\exp(\beta_1)\)은 unadjusted odds ratio
    • 로지스틱 다중선형회귀모형에서 \(\exp(\beta_1)\)은 adjusted odds ratio (adjusted for \(X_2, \ldots, X_p\))가 된다.

\[ \begin{aligned} & \frac{\operatorname{odds}\left(Y=1 \mid X=\left(x_{1}, \ldots x_{j-1}, x_{j}+1, x_{j+1} \ldots, x_{p}\right)^{T}\right)}{\operatorname{odds}(Y=1 \mid X=x)} \\ = & \frac{\exp \left\{\beta_{0}+\left(x_{1}, \ldots x_{j-1}, x_{j}+1, x_{j+1} \ldots, x_{p}\right) \beta\right\}}{\exp \left\{\beta_{0}+\left(x_{1}, \ldots x_{j-1}, x_{j}, x_{j+1} \ldots, x_{p}\right) \beta\right\}}=\exp \left(\beta_{j}\right) \end{aligned} \]

9.2.4 예제

Coefficient Std. error Z-statistic P-value
Intercept -10.8690 0.4923 -22.08 \(<0.0001\)
balance 0.0057 0.0002 24.74 \(<0.0001\)
income 0.0030 0.0082 0.37 0.7115
student[Yes] -0.6468 0.2362 -2.74 0.0062

TABLE 4.3 in [ISLR]. 로지스틱 회귀 모델의 추정 계수들의 예시. 이 모델의 설명변수는 계좌 잔액(balance), 소득(income), 그리고 학생여부(student yes)이다. 반응변수는 채무 불이행(default)이다. 소득은 천 달러 단위로 측정되었다.

  • 카드 청구예정금액(balance)이 \(\$ 1,500\) 이고 연간수입(income) 이 \(\$ 40,000\) 인 학생 (student[Yes])에 대한 채무불이행(Default) 확률 추정값은?

\[ \widehat{p}=\frac{e^{-10.869+0.00574 \times 1500+0.003 \times 40-0.6468 \times 1}}{1+e^{-10.869+0.000574 \times 1500+0.003 \times 40-0.6468 \times 1}}=0.058 \]

  • Exercises
  1. 카드 청구예정액 및 연간수입이 위와 동일한 비-학생의 채무불이행(디폴트) 확률 추정값은?
  2. 비-학생에 대한 학생의 디폴트 확률 오즈비의 추정값은?

9.3 모형 적합과 예측

  • 모형 계수 \((\beta_{0}, \boldsymbol{\beta})\) 의 추정에는 최대가능도법 (maximum likelihood)이 사용된다.
    • 최대가능도법은 확률모형에서 가장 보편적으로 사용되는 모수 추정법이다. 간단히 말해서, 최대가능도법으로 얻어진 모수의 추정량은 주어진 자료를 관찰할 확률을 가장 크게 한다. 7,8장의 최소제곱추정법은 (\(Y\)가 정규분포임을 가정할 때) 최대가능도법과 동치이다.
  • R에서의 문법
glm(y ~ x, family = binomial(link="logit")) 
glm(y ~ x, family = binomial(link="probit"))
# glm은 일반화선형모형(generalized linear model)의 줄임말이다.
  • 최대가능도법으로 얻어진 계수의 추정량을 \((\widehat{\beta}_{0}, \widehat{\boldsymbol{\beta}})\)라 쓰자.
  • 새로운 피처 \({\bf x}\) 에 대하여, 로짓모형은 먼저 \(Y=1\) 일 조건부확률을 아래와 같이 추정할 수 있다.

\[ \widehat{p}({\bf x}):=\widehat{P}(Y=1 \mid {\bf X}={\bf x}):=\frac{\exp \left(\widehat{\beta}_{0}+ {\bf x}^{T} \widehat{\boldsymbol{\beta}}\right)}{1+\exp \left(\widehat{\beta}_{0}+{\bf x}^{T} \widehat{\boldsymbol{\beta}}\right)} \]

  • 사전 지정된 cutoff value \(c\)에 따라 \(Y\)의 예측값을 1/0으로 가를 수 있다.

\[ \widehat{Y}:= \begin{cases}1 & \text { if } \widehat{p}({\bf x}) \geq c \\ 0 & \text { if } \widehat{p}({\bf x})<c\end{cases} \]

  • e.g. \(c=0.5\), \[ \widehat{Y}:= \begin{cases}1 & \text { if } \widehat{p}({\bf x}) \geq 0.5 \\ 0 & \text { if } \widehat{p}({\bf x})<0.5\end{cases} \]

9.4 예제: 한치록 예제 17.2

  • Mroz가 1987년 Econometrica 논문에서 사용한 자료(미국의 결혼한 여성 753 명의 경제활동에 관한 자료)
#install.packages("sampleSelection")
data(Mroz87, package="sampleSelection")
names(Mroz87)
 [1] "lfp"      "hours"    "kids5"    "kids618"  "age"      "educ"    
 [7] "wage"     "repwage"  "hushrs"   "husage"   "huseduc"  "huswage" 
[13] "faminc"   "mtr"      "motheduc" "fatheduc" "unem"     "city"    
[19] "exper"    "nwifeinc" "wifecoll" "huscoll"
  • 변수 lfp는 1975년에 경제활동에 참가한지 여부를 나타내는 이진변수이다. hours는 wife 가 일한 시간, kids5는 가구 내 6세 미만 자녀의 수, kids618은 가구 내 6-18세 자녀의 수, age는 wife의 나이, educ은 wife의 교육수준(연수), wage는 wife의 추정임금(1975 년 달러), faminc은 가구소득(1975년 달러)이다. exper는 wife의 근로경력(연수)이다. nwifeinc는 wife의 소득을 제외한 가구소득(1천 달러)이다.
str(Mroz87)
'data.frame':   753 obs. of  22 variables:
 $ lfp     : int  1 1 1 1 1 1 1 1 1 1 ...
 $ hours   : int  1610 1656 1980 456 1568 2032 1440 1020 1458 1600 ...
 $ kids5   : int  1 0 1 0 1 0 0 0 0 0 ...
 $ kids618 : int  0 2 3 3 2 0 2 0 2 2 ...
 $ age     : int  32 30 35 34 31 54 37 54 48 39 ...
 $ educ    : int  12 12 12 12 14 12 16 12 12 12 ...
 $ wage    : num  3.35 1.39 4.55 1.1 4.59 ...
 $ repwage : num  2.65 2.65 4.04 3.25 3.6 4.7 5.95 9.98 0 4.15 ...
 $ hushrs  : int  2708 2310 3072 1920 2000 1040 2670 4120 1995 2100 ...
 $ husage  : int  34 30 40 53 32 57 37 53 52 43 ...
 $ huseduc : int  12 9 12 10 12 11 12 8 4 12 ...
 $ huswage : num  4.03 8.44 3.58 3.54 10 ...
 $ faminc  : int  16310 21800 21040 7300 27300 19495 21152 18900 20405 20425 ...
 $ mtr     : num  0.722 0.661 0.692 0.781 0.622 ...
 $ motheduc: int  12 7 12 7 12 14 14 3 7 7 ...
 $ fatheduc: int  7 7 7 7 14 7 7 3 7 7 ...
 $ unem    : num  5 11 5 5 9.5 7.5 5 5 3 5 ...
 $ city    : int  0 1 0 0 1 1 0 0 0 0 ...
 $ exper   : int  14 5 15 6 7 33 11 35 24 21 ...
 $ nwifeinc: num  10.9 19.5 12 6.8 20.1 ...
 $ wifecoll: Factor w/ 2 levels " TRUE","FALSE": 2 2 2 2 1 2 1 2 2 2 ...
 $ huscoll : Factor w/ 2 levels " TRUE","FALSE": 2 2 2 2 2 2 2 2 2 2 ...

일하는 여성 \((\mathrm{lfp}=1)\) 이 428 명이다( 1 아니면 0 인 값들의 합은 1 의 개수와 동일하다).

sum(Mroz87$lfp)
[1] 428
  • 다음은 Mroz87 자료를 사용한, 여성의 경제활동 참여에 관한 로지스틱 회귀 추정 결과이다.
formula <- lfp ~ nwifeinc + educ + exper + I(exper^2) + age + kids5 + kids618
model <- glm(formula, data=Mroz87, family=binomial(link="logit"))
summary(model)

Call:
glm(formula = formula, family = binomial(link = "logit"), data = Mroz87)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.425452   0.860365   0.495  0.62095    
nwifeinc    -0.021345   0.008421  -2.535  0.01126 *  
educ         0.221170   0.043439   5.091 3.55e-07 ***
exper        0.205870   0.032057   6.422 1.34e-10 ***
I(exper^2)  -0.003154   0.001016  -3.104  0.00191 ** 
age         -0.088024   0.014573  -6.040 1.54e-09 ***
kids5       -1.443354   0.203583  -7.090 1.34e-12 ***
kids618      0.060112   0.074789   0.804  0.42154    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1029.75  on 752  degrees of freedom
Residual deviance:  803.53  on 745  degrees of freedom
AIC: 819.53

Number of Fisher Scoring iterations: 4
  • 여기서 model은 glm객체이다. (class() 함수로 확인 가능). 적합된 확률값은 다음 명령어들로 얻을 수 있다.
predict(model, type="response")
# 또는
fitted(model)
# 또는
model$fitted
  • 각 변수에 대한 adjusted OR의 추정값은, 적합계수와 신뢰구간에 \(\exp()\) 함수를 씌워주면 된다.
# adjusted OR
exp(coef(model)) 
(Intercept)    nwifeinc        educ       exper  I(exper^2)         age 
  1.5302825   0.9788810   1.2475360   1.2285929   0.9968509   0.9157386 
      kids5     kids618 
  0.2361344   1.0619557 
# adjusted OR의 신뢰구간
exp(confint(model))
Waiting for profiling to be done...
                2.5 %    97.5 %
(Intercept) 0.2829827 8.2907776
nwifeinc    0.9625447 0.9949193
educ        1.1473738 1.3607465
exper       1.1540225 1.3092803
I(exper^2)  0.9948707 0.9988906
age         0.8894633 0.9418229
kids5       0.1567171 0.3485513
kids618     0.9175642 1.2306237