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

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

주 출처

• 한치록 (2024). 계량경제학 강의 (5판). 박영사. – 17장 “이항반응모형”
• [ISLR] James, Witten, Hastie, and Tibshirani (2021). An Introduction to Statistical Learning: with Applications in R (2nd Edition). Springer. – 4장 “Classification”
• 김기영, 전명식, 강현철, 이성건 번역 (2015). Chatterjee and Hadi 지음. 예제를 통한 회귀분석 (5판). 자유아카데미.

보조 출처

• Wooldridge (2016). Introductory Econometrics: A Modern Approach (6th Edition). Cengage Learning.

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

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

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

9.2.1 로지스틱 선형모형

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

$$\begin{array}{}\text{(1)}& \log \left(\frac{p(\mathbf{x})}{1-p(\mathbf{x})}\right)={\beta }_0+{\mathbf{x}}^T\mathit{\beta }\end{array}$$ $$\text{ 또는 }p(\mathbf{x})=\frac{\exp ({\beta }_0+{\mathbf{x}}^T\mathit{\beta })}{1+\exp ({\beta }_0+{\mathbf{x}}^T\mathit{\beta })}$$

• 로짓함수 (logit function): $\mathrm{logit}(t):=\log \{t/(1-t)\}$
• $0<t<1,-\mathrm{\infty }<\mathrm{logit}(t)<\mathrm{\infty },t$ 의 단조증가함수
• ${\mathrm{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|\mathbf{X}=\mathbf{x})=P(Y=1|\mathbf{X}=\mathbf{x})=p(\mathbf{x})$이다.
• 따라서, 지난 7,8장에서 다룬 선형회귀분석 모형을 적용하자면 $p(\mathbf{x})={\beta }_0+{\mathbf{x}}^T\mathit{\beta }$ 식의 모형을 고려할 수 있고, 이 모형에서의 계수추정은 이미 7,8장에서 다룬 바 있다. 이 모형을 이진반응변수에 적합하면 안 되는 걸까?
• 결론적으로 “완벽히 안 된다”라고 말할 수는 없지만, 여러 가지 문제가 발생한다. 가장 큰 문제는 조건부확률의 적합값에서 발생한다. $P(Y=1|\mathbf{X}=\mathbf{x})={\beta }_0+{\mathbf{x}}^T\mathit{\beta }$의 좌변은 확률이니 $\mathbf{x}$의 값에 상관없이 늘 0과 1사이에 있어야 한다. 그런데 ${\beta }_0+{\mathbf{x}}^T\mathit{\beta }$ 항은 $\mathbf{x}$의 값에 따라 0보다 작을 수도 있고 1보다 클 수도 있다. 즉, 선형회귀분석 모형에서는 $p(\mathbf{x})$의 적합값이 0보다 작거나 1보다 클 수도 있다는 문제점이 있다. (아래 그림의 첫번째)
• $p(\mathbf{x})$의 적합값이 반드시 0과 1사이에 오게 하는 방법이 있다. 어떤 알려진 변환 $g$에 대하여 $g(p(\mathbf{x}))={\beta }_0+{\mathbf{x}}^T\mathit{\beta }$ 라고 가정하는 것이다. 이때 $g$ 함수는 $p_{\mathbf{x}}$ 를 선형함수인 ${\beta }_0+{\mathbf{x}}^T\mathit{\beta }$ 와 연결시키는 함수로서 ‘연결함수’(link function)라 한다.
• 로지스틱 회귀분석모형의 link function은 logit function $\log \{t/(1-t)\}$ 이다.
• 다른 대표적 link function은 probit function 으로, $t\mapsto {\int }_{-\mathrm{\infty }}^t\varphi (u)du$의 역함수로 정의된다 (단, $\varphi (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{array}{rl}\text{ odds }:& =\frac{p}{1-p}\\ \mathrm{odds}(Y=1\mid X& =\mathbf{x}):=\frac{p(\mathbf{x})}{1-p(\mathbf{x})}\end{array}$$

• 로지스틱 다중회귀모형에서 ${\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,\dots ,X_p$)가 된다.

$$\begin{array}{rl}& \frac{\mathrm{odds}(Y=1\mid X={(x_1,\dots x_{j-1},x_j+1,x_{j+1}\dots ,x_p)}^T)}{\mathrm{odds}(Y=1\mid X=x)}\\ =& \frac{\exp \{{\beta }_0+(x_1,\dots x_{j-1},x_j+1,x_{j+1}\dots ,x_p)\beta \}}{\exp \{{\beta }_0+(x_1,\dots x_{j-1},x_j,x_{j+1}\dots ,x_p)\beta \}}=\exp \left({\beta }_j\right)\end{array}$$

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)이 $\mathrm{\$}1,500$ 이고 연간수입(income) 이 $\mathrm{\$}40,000$ 인 학생 (student[Yes])에 대한 채무불이행(Default) 확률 추정값은?

$$\hat{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,\mathit{\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)의 줄임말이다.

• 최대가능도법으로 얻어진 계수의 추정량을 $({\hat{\beta }}_0,\hat{\mathit{\beta }})$라 쓰자.
• 새로운 피처 $\mathbf{x}$ 에 대하여, 로짓모형은 먼저 $Y=1$ 일 조건부확률을 아래와 같이 추정할 수 있다.

$$\hat{p}(\mathbf{x}):=\hat{P}(Y=1\mid \mathbf{X}=\mathbf{x}):=\frac{\exp ({\hat{\beta }}_0+{\mathbf{x}}^T\hat{\mathit{\beta }})}{1+\exp ({\hat{\beta }}_0+{\mathbf{x}}^T\hat{\mathit{\beta }})}$$

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

$$\hat{Y}:=\{\begin{array}{ll}1& \text{ if }\hat{p}(\mathbf{x})\ge c\\ 0& \text{ if }\hat{p}(\mathbf{x})<c\end{array}$$

• e.g. $c=0.5$, $$\hat{Y}:=\{\begin{array}{ll}1& \text{ if }\hat{p}(\mathbf{x})\ge 0.5\\ 0& \text{ if }\hat{p}(\mathbf{x})<0.5\end{array}$$

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