로지스틱 회귀분석

4. 로지스틱 회귀분석

로지스틱 회귀분석(logistic regression analysis)은 일반화 선형모형(generalized linear model, GLM)이라 불리는 큰 범주의 통계모형 모델링 방법에 속하는 방법입니다. 우선 GLM의 특징만 간단히 훑어보고 로지스틱 회귀모형에 대해 다루겠습니다.

GLM(Generalized Linear Model)

• GLM은 문자 그대로 선형적이지 않은 대상(비선형)을 선형적으로 '일반화'시킨 모형입니다.

• 선형화 시키는 이유는 여러 가지가 있을 수 있지만, 가장 대표적으로 선형모형에서만 사용할 수 있는 모형의 해석, 확장, 수정 등의 방법을 사용하기 위해서입니다.

• 비선형모형의 경우는 모형을 다루는 방법이 많이 제한될 뿐만 아니라 새로운 데이터에 민감하기 때문에 선형모형에 비해 덜 선호되는 경향이 있습니다.

그렇다면 로지스틱 회귀분석에서 선형화 시키는 대상은 무엇일까요? 바로 관심 범주에 속할 확률입니다. 확률은 일반적으로 예측자들에 따라 비선형하게 분포되어 있고, 형태는 S자 형태가 되는 경우가 많습니다. 물론 처음부터 확률에 선형라인을 적합시켜 선형회귀분석을 진행할 수도 있지만 이는 몇 가지 제약이 따릅니다.

library(ggplot2)

k = runif(n = 500,min = 2 , max = 4)

K = c()
P = c()
P1 = c()

for(i in k){

  p1 = exp(-21.3 + 6.74*i) /(1 + exp(-21.3 + 6.74*i)) 
  # error = runif(n = 1, min = -0.1, max = 0.1)
  error = 0 
  p = p1 + error
  p = ifelse(p < 0 ,0,p)
  p = ifelse(p > 1, 1,p)


  K = c(K,i)
  P = c(P,p)
  P1 = c(P1,p1)

}

ggplot(NULL) +
  geom_point(aes(x = K, y = P),col = 'royalblue', alpha = 0.2) +
  theme_bw() +
  scale_x_continuous(breaks = seq(2,4,by = 0.5)) +
  xlab("Predictor") + ylab("Prob")

위 그래프는 어떤 연속형 예측자에 따른 관심범주에 속할 확률입니다. 예를 들어 알콜 섭취량과 비만일 확률은 이런 식으로 분포할 수 있습니다. 알콜 섭취가 많을 수록 비만일 확률은 높아지지만 완전 선형적이라기 보다는 약간의 커브가 존재하는 것입니다. 우선여기에 일반적인 선형회귀 라인을 적합시켜 보겠습니다.

ggplot(NULL) +
  geom_point(aes(x = K, y = P),col = 'royalblue', alpha = 0.2) +
  geom_smooth(aes(x = K, y = P),method = 'lm', col = 'grey20') +
  geom_hline(yintercept = c(0,1), linetype = 'dashed') +
  scale_x_continuous(breaks = seq(2,4,by = 0.5)) +
  theme_bw() + 
  xlab("Predictor") + ylab("Prob")

선형 회귀선이 아주 비합리적이진 않습니다. 연속형 예측자가 증가함에 따라 확률 역시 증가하기 때문이죠. 또한 이런 선형식의 경우 위에서 배운 분산분석을 통한 모형의 적합도 검정부터 시작해서 기존에 회귀분석에서 사용하던 모든 분석을 진행할 수 있습니다. 그렇지만 이런 경우에는 구조적인 문제가 발생합니다. 바로 확률은 0과1 사이의 값을 갖는 데 반해 회귀선은 경우에 따라 음의 값이나 1을 초과하는 예측값를 제시할 수도 있다는 것입니다. 영역을 표시하자면 아래 동그라미와 같은 구간입니다.

ggplot(NULL) +
  geom_point(aes(x = K, y = P),col = 'royalblue', alpha = 0.2) +
  geom_smooth(aes(x = K, y = P),method = 'lm', col = 'grey20') +
  geom_hline(yintercept = c(0,1), linetype = 'dashed') +
  geom_point(aes(x = c(2.25,3.9), y = c(-0.05,1.05)),size = 10,shape =  1, col = 'red') +
  scale_x_continuous(breaks = seq(2,4,by = 0.5)) +
  theme_bw() + 
  xlab("Predictor") + ylab("Prob")

세로축은 확률 값을 표현한 것인데, 빨간 동그라미는 확률이 가질 수 없는 값을 추정하고 있습니다. 위 데이터에서는 -4 이하의 점에서 -확률 값으로 추정되고 있지만, 이 회귀선을 이용하여 새 데이터를 예측해 보면 이보다 훨씬 많은 포인트가 잘못된 값으로 추정될 가능성이 있습니다. 또한 이러한 구조적인 문제로 확률의 특성을 이용한 기타 추가 분석 역시 불가능하게 됩니다.

그렇다면 이제 비선형적인 모형을 적합시켜보겠습니다. 아래 그래프는 S자 형태를 띠는 비선형적인 회귀곡선입니다.

ggplot(NULL) +
  geom_point(aes(x = K, y = P),col = 'royalblue', alpha = 0.2) +
  geom_line(aes(x = K, y = P1),col = 'red', size = 2.5, alpha  = 0.5) +
  geom_hline(yintercept = c(0,1), linetype = 'dashed') +
  scale_x_continuous(breaks = seq(2,4,by = 0.5)) +
  theme_bw() + 
  xlab("Predictor") + ylab("Prob")

한눈에 보기에도 선형회귀선보다는 훨씬 데이터에 잘 적용됨을 알 수 있고 구조적인 문제도 발생하지 않습니다. 그렇지만 위에서 말씀드린대로 비선형 모델은 여러 가지 추가 분석에 제약이 있습니다. 이러한 제약을 극복하기 위해서 약간의 변환을 통해 이를 선형화하는 것이 로지스틱 회귀분석의 기본 아이디어입니다.

예측자가 주어졌을 때, 관심범주에 속할 확률을 선형으로 예측하기 이해서는 $0 \leq \pi_x \leq1$의 제약을 $-\infty<log(\frac{\pi_x}{1-\pi_x})<\infty$의 형태로 변환을 해주어야 합니다. $\pi_x$를 $\frac{\pi_x}{1-\pi_x}$로 변환, 즉 $odds$로 변환을 해주면, $\frac{\pi_x}{1-\pi_x}$는 $0 < \frac{\pi_x}{1-\pi_x}<\infty$의 범위를 가지게 됩니다. 그 후 $log$변환을 취해주면, 로그 함수의 특성에 따라 $-\infty<log(\frac{\pi_x}{1-\pi_x})<\infty$의 범위를 가지게 되는 것입니다. 우리는 관심범주의 확률을 선형모형으로 예측하기 위해 순수 $\pi_x$를 활용하는 것이 아니라, $log(\frac{\pi_x}{1-\pi_x})$를 활용할 것입니다. 이 변환을 $logit$변환이라 부르며, $logit$변환 값을 반응변수롤 하고 선형모델을 적합시킨 것을 로지스틱회귀분석이라고 합니다.

선형화를 위해 아래의 식을 고려해보겠습니다. 여기서 $\pi_x$는 설명 예측자가 주어졌을 때의 관심 범주에 속할 확률입니다.

$$
\pi_x = \frac {exp(\beta_0 + \beta_1 x )} {1+exp(\beta_0 + \beta_1 x )}
$$

이는 일반 선형회귀에서의 $\pi_x = \beta_0 + \beta_1 x$ 와 같이 직선의 관계를 나타내는 함수에 대응하는 것으로, 선형적인 관계가 아닌 성공확률 $\pi_x$ 와 설명 예측자와의 관계를 S자 모양으로 표현해주는 함수입니다. 또한 선형화를 위해 위의 함수를 이용하여 선형식 모양으로 정리를하면 다음과 같습니다.

$$
logit(\pi_x) = log \left( \frac {\pi_x} {1-\pi_x} \right) = \beta_0 +\beta_1x
$$

$log \left(\frac{\pi_x}{1-\pi_x} \right)$ 와 같은 형태를 $\pi_x$ 에 대한 로짓함수라고 부르며 $logit(\pi_x)$ 와 같이 표현합니다. 이 로짓함수를 이용하면 우변을 선형식으로 표현할 수 있습니다. 이것이 처음에 말씀드렸던 확률과 예측자의 비선형적 관계를 선형식으로 일반화하는 과정입니다. 또한 로지스틱 회귀모형과 같은 GLM에서는 절편과 기울기를 일반적인 회귀분석과 다르게 최소제곱 추정법이 아닌 최대가능도도추정법(method of maximum likelihood estimation)이라 불리는 방법에 의해 추정합니다. 이는 다음 장에서 다루겠습니다

이제 구체적으로 로지스틱 모형을 살펴봅시다.

기본적인 로지스틱 모형은 반응 값이 이항변수임을 가정합니다. 잘 생각해보면, 위에서 계속 다루었던 확률은 주어진 값이 될 수 없다는 것을 떠올리실 수 있을 겁니다. 확률은 측정될 수 있는 것이 아니기 때문이죠. 설명을 위해 보여드린 그래프의 확률 값 역시 마찬가지입니다.

로지스틱 모형에서의 확률은 이항 반응 값으로 측정된 질적 자료를 이용해 추정됩니다. 예를 들어 알콜 섭취량과 비만일 확률의 관계를 알기 위한 조사를 할 때, '비만일 확률' 자체는 측정할 수 없고 '비만 여부' 라는 이항 반응 값을 측정해 그 확률을 추정한다는 것이죠.

이런 이항 반응 값(관심 범주 = 1)와 예측자간의 산점도를 그리면 다음과 같습니다.

PG = ifelse(P1 < 0.5,0, 1)
K1 = K + runif(n = length(K),min = -0.5, max = 0.5)

ggplot(NULL) +
  geom_point(aes(x = K1, y = PG),col = 'royalblue') +
  geom_line(aes(x = K, y = P1),col = 'red', size = 2.5, alpha  = 0.5) +
  theme_bw() +
  xlab("Predictor") + ylab("Prob") +
  scale_x_continuous(breaks = seq(2,4,by = 0.5)) 

관심범주 아니면 비 관심범주로 조사된 세로 축 반응 변수는 0아니면 1의 값만 갖습니다. 관심있는 범주면 1을 부여하고 그렇지 않으면 0을 부여해 구별해 놓은 것이지요. 그에 비해 예측자는 연속적으로 분포함을 알 수 있습니다. 산점도를 보면 직관적으로 예측자 값이 큰 값을 가지면 관심범주에 속할 확률이 크고 작으면 비 관심 범주에 속할 확률이 크다는 것을 확인할 수 있습니다. 문제는 예측자가 겹치는 부분입니다. 하늘색 선을 기준으로 가장 왼쪽, 오른쪽 영역은 겹치는 부분이 없어 어렵지 않게 범주를 예측할 수 있지만 두 범주가 공존하는 가운데 부분은 예측자가 겹치게됩니다. 이 부분이 확률을 추정하여 로지스틱 회귀곡선을 그렸을 경우 변곡점이 생기는 영역이고 추정된 확률이 0.5에 가까운 애매한 값으로 표현되는 곳입니다. 즉, 예측하기 어려운 영역이라고 할 수 있습니다. 로지스틱 회귀곡선을 적합시켜 보겠습니다.

ggplot(NULL) +
  geom_point(aes(x = K1, y = PG),col = 'royalblue') +
  geom_vline(xintercept = c(2.7,3.6),linetype = 'dashed', col = 'red') +
  geom_text(aes(x = c(2,4), y = c(0.5,0.5)),label = c("비관심범주","관심 범주"),col = 'red',
            size = 5) +
  geom_text(aes(x = 3.1, y = 0.5), label = "공존", col = 'royalblue',
            size = 5) +
  theme_bw() +
  xlab("") + ylab("") +
  scale_x_continuous(breaks = seq(1,5,by = 0.5)) 

적합된 회귀곡선은 다음과 같습니다.

$$
\hat \pi_x = \frac {exp(-21.3+6.74x)} {1+exp(-21.3+6.74x)} \ \Leftrightarrow\ log\left( \frac {\hat \pi_x} {1-\hat \pi_x} \right) = -21.3+6.74x
$$

선형적으로 추정된 $logit(\hat \pi)$ 는 꼭 연속형 예측자뿐 아니라 범주형 예측자도 가능하고 혼합된 형태도 가능합니다. 또한 모형에 대한 해석은 선형적으로 바라보아도 상관은 없으나 로짓 값이 직관적이지 않아 일반적으로 오즈비를 통해 하게됩니다.

예를 들어 위 식에서 예측자의 기울기 6.74는 '예측자가 한 단위(1) 증가했을 때의 성공할 오즈 '가 됩니다. 이것은 두 모형의 차를 통해 확인할 수 있습니다.

예측자가 $x$ 인 모형과 한 단위(1) 증가하여 $x+1$ 인 모형의 차를 봅시다. 아래 식에서는 추정된 확률의 구별을 위해 예측자가 $x$ 인 경우에는 $\hat{\pi_{x}}$, $x+1$ 인 경우는 $\hat{\pi_{x+1}}$ 로 표현하겠습니다.

$$
\left[; -21.3+6.74(x+1) ;\right] - \left[; -21.3+6.74(x) ;\right] = 6.74\
$$

$$
\ \quad \
\Leftrightarrow \ \left[log\left( \frac {\hat \pi_{x+1}} {1-\hat \pi_{x+1}} \right) \right]-\left[log\left( \frac {\hat \pi_x} {1-\hat \pi_x} \right) \right] =log \left[ \left(\frac {\hat \pi_{x+1}} {1-\hat \pi_{x+1}} \right)/\left( \frac {\hat \pi_x} {1-\hat \pi_x} \right)\right]= log(\widehat {o.r} )
$$

즉, 로지스틱 회귀모형에서의 기울기는 예측자가 1 증가했을 때의 로그오즈비 추정량이고 $exp[기울기]$ 는 예측자가 한 단위 증가했을 때의 오즈비 추정량입니다. 위의 예에서 적용해보면 예측자가 한 단위 증가할 때 관심범주에 속할 오즈가 $exp(6.74) = 846$ 배가 된다는 것을 알 수 있습니다. 이를 응용하면 단위를 조정해서 확인할 수도있습니다. 예를 들어 위에서 다룬 예제는 예측자의 범위가 매우 작으므로 1이 아닌 0.1이 증가했을 때의 오즈비를 확인하고 싶을 수 있습니다. 로그오즈비는 $\left[\ -21.3+6.74(x+0.1) \right] - \left[\ -21.3+6.74(x) \right] = 0.674$ 이 되겠고 오즈는 $exp(0.674)=1.96$ 이 됩니다. 즉, 예측자가 0.1 증가하면 관심 범주에 속할 오즈가 1.96배 됨을 알 수 있습니다.