의학 통계 코드 저장소

[이론] 가우스-마르코프 정리 (Gauss-Markov Theorem)

이전 최소제곱법 글(2023.06.15 - [통계 이론] - [이론] 보통최소제곱법 (Ordinary Least Squares))에서 회귀 계수를 구하는 방법에 대해 알아보았다. 가우스-마르코프 정리는, 최소제곱법에서 추정한 회귀 계수 값들이 특정한 조건 하에서 가장 좋은 추정량이 된다는 정리다. 그 조건들은 다음과 같다.

(1) 선형성 (Linearity) : Linearity in parameters

회귀계수가 선형이어야 한다. 즉, $Y=aX+b$와 같은 관계여야 한다. 여기에서 중요한 것은 a지 X가 아니다. a가 선형이어야 하지, X는 선형이 아니어도 된다. 즉, $Y=aX^2+b$와 같은 관계도 선형이다. 심지어 $Y=af(X)+b$도 회귀계수에 대해서는 선형이다. 

(2) 독립변수의 비확률성(Nonstochasticity) : Nonstochasticity of regressors (independent variables)

 독립변수는 미리 정해져 있어야 한다. 즉, 독립변수는 확률변수면 안 된다. 풀어 말하면, 독립변수의 값이 구해져 있는 상태에서 종속변수의 값을 측정해야 한다. 하지만, 일반적인 연구를 할 때에는 독립변수와 종속변수를 동시에 측정한다. 즉, 독립변수가 미리 정해져 있지 않고, 달리 말하면 확률변수라고 말할 수 있다. 이런 경우에는 다음과 같은 조건으로 조금 완화하기도 한다. "독립변수의 값과 잔차의 값은 관련이 없어야 한다." 

위 그림에서 $x_i$값과 오렌지색 선분의 길이는 상관이 없어야 한다. 만약, $x_i$값이 커짐에 따라 오렌지색 선분의 길이가 길어진다면 그것은 관련이 있다는 뜻이다. 예를 들어, 나이와 지식수준의 관련성을 확인하는 상황을 생각해 보면, 3살짜리 아이들의 지식수준은 비슷비슷한데(오렌지색 선분이 짧음), 성인의 차이가 많이 나므로 비상관성이 성립하지 않을 수 있다. (이는 표준화 등의 방법으로 해결할 수 있긴 하다.) 

(3) 오차의 평균=0

  관측치들은 회귀선 위에도, 아래에도 고르게 분포하여 오차의 평균은 0이 되어야 한다.

$$ \forall_i, E \left( \epsilon_j | x_i \right)=0$$

즉, 독립변수가 특정 값인 관측치들만 따로 뽑아 종속변수의 평균을 구하면 회귀선 위에 위치한다는 뜻이다.

$$ \forall_i, E\left(y_j | x_i \right) =  \alpha + \beta x_i $$

(4) 등분산성 (Homoscedasticity): 오차의 등분산성

 대표적인 대푯값으로는 평균과 표준편차가 있다. 오차의 평균은 위에서 다루었으니, 표준편차도 다루어보도록 한다. 표준편차의 제곱이 분산이므로 표준편차를 다루는 것은 분산을 다루는 것과 같다.

 독립변수의 값에 관계없이, 종속변수가 퍼져있는 정도는 비슷해야 한다. 즉, 오차가 퍼져있는 정도는 독립변수의 값에 관계없이 일정해야 한다. 따라서 적절한 상수 $\sigma$에 대해 다음 수식이 성립해야 한다. 그리고, $\sigma$를 표준편차라 한다.

$$\forall_i, var \left( \epsilon_j | x_i \right) = \sigma^2 $$

 만약 등분산성이 성립하지 않는다면, 즉 $var \left( \epsilon_j | x_i \right)=\sigma_i ^2$이라면, 이분산성 (heteroscedasticity)을 보인다고 한다.

(5) 비자기상관성 (Nonautocorrelation)

 한 오차로 다른 오차값을 예상할 수 있으면 안 된다. 각 개인이 10살일 때, 20살일 때의 키에 대한 데이터가 있다고 하자. 10살일 때 키가 평균보다 컸던 사람(오차가 큰 사람)은 20살일 때에도 클(오차가 클) 확률이 높다. 즉, 한 오차값을 알면 다른 오차값을 예측할 수 있게 된다. 이런 경우 자기상관(autocorrelation)이 있다고 하며, 이런 경우는 적절하지 않다. 이를 수식으로 표현하면 다음과 같다.

 $$\forall_{i,k}, cov\left(\epsilon_j, \epsilon_l | x_i , x_k \right)=0$$

(6) 다중공선성 없음 (No multicollinearity)

 마지막으로, 다중공선성이 없어야 한다. 보통 다중공선성은 공선성의 개념을 포함하는 것처럼 통용되므로, 본 글에서 다중공선성이 없다는 것은, 공선성도 없다는 것을 내포하기로 한다. 먼저, 공선성과 다중공선성에 대해서 알아보기로 한다.

 공선성이란 "모델에 있는 두 독립변수 간에 긴밀한 상관관계가 있는 것"으로 정의할 수 있다. 예를 들어, 종속변수가 "삶의 질"일 때, 독립변수로 "월 수입"과 "연 수입" 두 개를 넣는 상황을 상정해 보자. 이때, "연 수입 = 월 수입 * 12"라는 정확한 상관관계가 존재하므로, 공선성이 존재한다고 할 수 있고, 두 개의 변수를 모델에 모두 넣으면 안 된다.

 다중공선성이란 "모델에 있는 세 개 이상의 독립변수 간에 일종의 선형 관계가 있는 것"으로 정의할 수 있다. 예를 들어, 독립변수로 "출생 연도", "나이", "현재 년도" 세 가지를 쓰는 상황을 생각해 보자. 이때 "나이=현재 년도 - 출생 연도"라는 정확한 선형 관계가 존재하므로, 다중공선성이 존재한다고 할 수 있다. 따라서 이 세 변수를 모두 독립변수로 사용하려 하면 안 된다. (이런 경우에 세 변수를 모두 독립변수로 사용하기 위한 방법으로 Age-Period-Cohort Analysis가 존재하고, 이 다중공선성을 해결하기 위해 여러 방법론들이 제시되었다. 다음 링크에서 내용을 확인할 수 있다. 링크 추가 예정)

여기까지가 가우스 마르코프 정리의 내용이다. 통계 공부를 좀 해보신 분들이라면 알 수 있겠지만, 이는 선형 회귀 분석의 가정으로 이어진다. 그리고 이런 모델을 Classical Linear Regression Model이라고 한다.

SPSS : 링크 추가 예정

SAS : 링크 추가 예정

R: 2022.12.22 - [선형 회귀 분석/R] - [R] 선형 회귀 분석의 전제 조건 - lm(), plot(), gvlma()

 여기서 의문이 든다. 그렇다면 왜 선형 회귀 분석의 전제조건으로 "잔차의 정규성"을 요구할까?

이는 가설 검정을 통해 유의성, 즉 p-value를 산출하기 위함이다. 정확히 말하면, 잔차의 정규성을 만족하는 순간, 보통 최소 제곱법은 최대가능도와 같은 것이 되므로, 최소제곱법 대신 최대가능도 방법을 사용할 수 있다. 각각에 대한 내용은 다음 링크를 확인하길 바란다.

최소제곱법: 2023.06.15 - [통계 이론] - [이론] 보통최소제곱법 (Ordinary Least Squares)

최대가능도: 2022.11.25 - [통계 이론] - [이론] 로지스틱 회귀분석에서 회귀 계수를 구하는 방법 - Maximum likelihood estimation

  그리고, 잔차의 정규성까지 만족하는 모델을 Classical Normal Linear Regression Model이라고 한다.

[이론] 가우스-마르코프 정리 (Gauss-Markov Theorem) 정복 완료!

작성일: 2023.06.21.

최종 수정일: 2023.06.21.

[이론] 보통최소제곱법 (Ordinary Least Squares)

 의학통계를 비롯한 많은 통계에서 모수(parameter)를 추정할 때 가장 많이 쓰이는 방법 중 하나가 보통최소제곱법(Ordinary Least Squares)이다. 이에 대해 알아보고자 한다. 보통최소제곱법 (Ordinary Least Squares)은 최소제곱법 (Least Squares)와 다른 것이지만 많이들 구분하지 않고 사용하고 있다. 본 글에서는 철저히 구분하여 사용할 것이다.

최소제곱법 (Least Squares)

  최소제곱법의 분류를 이해하기 위해서는 최소제곱법이 무엇인지 알아야 한다. 

 네 개의 표본으로부터 얻은 데이터를 다음과 같이 좌표평면 위에 표시했다고 하자.

누군가가 "이 데이터는 적절한 선형 관계로 표시할 수 있을 것 같아."라고 생각하였고, 적절한 선을 다음과 같이 그렸다고 하자.

 다른 누군가는 이 직선이 아니라, 기울기와 $y$절편이 조금은 다른 직선이 더 적절하다고 주장할 수도 있다. 그렇다면 어떤 기준으로 직선을 정해야 다수가 받아들일 수 있을까? 그에 대한 답 중 하나가 최소제곱법(Least Squares)이다.

최소 제곱법은 직선 관계에서만 답을 내어줄 수 있는 것은 아니다. 또 다른 누군가는 직선이 아니라 아래와 같은 이차함수가 더 적절하다고 주장할 수 있고, 다른 누군가는 오른/왼쪽으로 조금 이동한 이차함수가 더 적절하다고 주장할 수 있다.

 최소제곱법 (Least Squares)은 이런 문제들에 대한 답을 내어준다. 직선 관계로 표현하기로 마음먹었다면, 어떤 직선이 데이터를 가장 잘 설명하는지 답을 내어줄 수 있다. 이차 함수 관계로 표현하기로 했다면, 어떤 이차 함수가 가장 잘 설명하는지 알려줄 수 있다. 

 요약하면, 최소제곱법은 잔차를 제곱한 것들의 합이 최소화되는 순간을 찾는 테크닉이다. 방법은 다음과 같다.

1) 데이터를 잘 설명할 것으로 예상되는 임의의 직선을 그린다. 

2) 그 직선 위 $x$좌표가 데이터의 $X_1 (x_1,y_1)$의 x좌표$(x_1)$와 같은 점$(\hat{X_1})$의 $y$값$(\hat{y_1})$은 직선 관계를 사용하여 구한 예측치다.

3) 실제 $y$값$(x_1)$과 예측치$(\hat{y_1})$의 차이를 잔차라고 한다. (오렌지색 선분)

4) 오렌지색 선분들이 짧으면 짧을 수록, 직선이 데이터를 잘 설명한다고 할 수 있다. 그런 목적으로 잔차의 합을 구하면 문제가 발생한다. 데이터들은 직선 위에도, 아래에도 있으므로 잔차는 +, - 두 가지가 모두 존재한다. 따라서 단순 합산은, 오렌지색 선분이 아무리 길어도 잘 상쇄되기만 하면 0이 될 수도 있다. 그러므로 그냥 더하는 것이 아니라 제곱을 하여 더하게 된다. 이는 통계에서 분산이 편차의 합이 아닌, 편차의 제곱의 합인 이유와 같다.

5) 편차의 제곱의 합을 최소로 만드는 기울기와 $y$절편을 구하면 된다.

최소제곱법 (Least Squares)의 분류

사실 위 설명은 최소제곱법 중 하나인 보통최소제곱법 (Ordinary Least Squares, OLS)에 대한 설명이다. 최소제곱법의 분류는 다음과 같다.

1. Linear Least Squares : 선형

 - Ordinary Least Squares

 - Weighted Least Squares : 등분산성이 보장되지 않을 때

 - Generalized Lieast Squares : 등분산성뿐만 아니라 자기상관 등 문제가 있을 

2. Non-Linear Least Squares

보통최소제곱법 (Ordinary Least Squares, OLS)의 계산법

 지금부터는 실제로 잔차 제곱의 합을 최소화하는 기울기와 $y$절편을 구해보고자 한다.

직선의 방정식은 다음과 같다고 하자

$$y=ax+b$$

데이터는 $n$개의 표본으로 이루어져 있으며, 각 데이터는 다음과 같이 표현된다 

$$X_{i} (x_i, y_i), where \; i\in\mathbb{N}\; AND\; 1\leq i\leq n $$

그럼 $X_{i} (x_i, y_i)$의 잔차 $d_i$는 다음과 같이 계산된다.

$$d_i=y_i - (ax_i +b) $$

잔차의 제곱의 합은 다음과 같다.

$$\sum_{i=1}^{n} {\left(y_i -ax_i-b \right)^2}$$

이 값이 최소화되는 지점에서는 $a$에 대해서든, $b$에 대해서든 편미분을 하면 0이 된다.

$$ \frac {\partial} {\partial a} \sum_{i=1}^{n} {\left(y_i -ax_i-b \right)^2} =-2\sum_{i=1}^{n} {x_i \left(y_i - ax_i -b \right) } = 0 \cdots (1) $$

$$ \frac {\partial} {\partial b} \sum_{i=1}^{n} {\left(y_i -ax_i-b \right)^2} =-2\sum_{i=1}^{n} {\left(y_i - ax_i -b \right)}=0 \cdots (2) $$

(2)식에서 다음을 얻는다.

$$ \begin{align} b&= \frac {\sum_{i=1}^{n} y_i} {n} - a \frac {\sum_{i=1}^{n} x_i}{n} \\ &=\bar{Y} - a\bar{X} \\&&\end{align} $$

이를 (1) 식에 대입하면 다음을 얻는다.

$$ \begin{align} a &= \frac {\sum_{i=1}^{n}x_i \left( y_i - \frac {\sum_{i=1}^{n} y_i } {n} \right)} {\sum_{i=1}^{n} x_i \left(x_i - \frac {\sum_{i=1}^{n} x_i} {n} \right)} \\ &= \frac {\sum_{i=1}^{n} x_i \left( y_i - \bar{Y} \right)} {\sum_{i=1}^{n} x_i \left( x_i - \bar{X} \right)} \cdots(3) \\&&\end{align}$$

한편, 

$$ \sum_{i=1}^{n} \left( \bar{X}^2 - x_i \bar{X} \right)= \bar{X} \sum_{i=1}^{n} \left( \bar{X} - x_i \right)=0$$

$$ \sum_{i=1}^{n} \left( \bar{X}y_i - \bar{X}\bar{Y} \right)=\bar{X} \sum_{i=1}^{n} \left( y_i - \bar{Y} \right)=0$$

이므로 (3)식의 분자, 분모는 다음과 같이 쓸 수 있다.

$$\sum_{i=1}^{n} x_i \left(y_i - \bar{Y} \right) = \sum_{i=1}^{n} \left[ x_i \left(y_i - \bar{Y}\right)- \left(\bar{X} y_i - \bar{X}\bar{Y} \right)\right] = \sum_{i=1}^{n} \left(x_i - \bar{X} \right) \left( y_i - \bar{Y} \right) = (n-1) cov \left(X,Y \right) $$

$$\sum_{i=1}^{n} x_i \left( x_i - \bar{X} \right)=\sum_{i=1}^{n} \left[ x_i  \left( x_i - \bar{X} \right) + \left( -x_i\bar{X} + \bar{X}^2 \right) \right]=\sum_{i=1}^{n} \left( x_i - \bar{X} \right)^2 =(n-1) var(X)$$

따라서 다음을 얻는다.

$$a= \frac {cov \left(X,Y \right)} {var(X)}$$

$$b=\bar{Y} - \bar{X} \frac {cov \left(X,Y \right)} {var(X)}$$

값이 하나라는 점에서도 알 수 있지만, Local minima가 존재하지 않음을 알 수 있고, 이는 바로 global minima임을 알 수 있다.

이렇게, 선형 회귀분석에서 적절한 회귀 계수들을 찾는 방법인 최소제곱법에 대해 알아보았다. '최소제곱법을 통해 구한 값이 정말 최적일까?'라는 질문에 대한 답변은 '가우스-마르코프 정리 Gauss-Markov Theorem'을 통해 얻을 수 있다. 해당 내용은 다음 링크에서 알아볼 수 있다. 링크 추가 예정 

[이론] 보통최소제곱법 (Ordinary Least Squares) 정복 완료!

작성일: 2023.06.15.

최종 수정일: 2023.06.15.

[이론] 로지스틱 회귀분석에서 회귀 계수를 구하는 방법 - Maximum likelihood estimation

 선형 회귀 분석에서 회귀 계수를 구하는 방법은 최소 제곱법으로 비교적 직관적이다. 데이터를 가장 잘 설명할 수 있는 직선을 잘 그리면 되는 것이다. 즉 아래 그림에서 오렌지색 선분 길이 제곱의 합이 최소화되는 직선을 찾는 것이다.

 하지만 로지스틱 회귀분석 (logistic regression)에서 회귀 계수를 같은 방법으로 구하는 것은 적절하지 않다. 예를 들어 혈중 ALT에 따른 고혈압 여부에 대한 로지스틱 회귀분석을 실시한다고 하자. 이때 로지스틱 회귀분석에서 결과변수는 이분형 (예시: 고혈압 or 정상)이므로 그래프를 그리면 두 개의 선으로 나타나게 된다. 이런 데이터를 가장 잘 설명하는 하나의 선을 구하는 것은 적절하지 않아 보인다. 따라서 로지스틱 회귀분석에는 maximum likelihood estimation (MLE)라는 방법으로 회귀 계수를 구한다.

 1. Maximum likelihood estimation이란 무엇인가?

상황 A: 앞면이 나올 확률이 $p$인 동전을 10회 던졌는데 6회 앞면이 나왔다고 하자.

이런 상황이 발생했을 확률은 다음과 같이 표현할 수 있다.

$$ L=p^6 (1-p)^4$$

  이 식은 상황 A가 발생할 확률이자 가능도 (likelihood)를 의미한다. 우리는 상황 A를 표현한 위 식에서 $p$값을 유추하고자 한다. 어떻게 해야 할까? 만약 특정 $p$값일 때 10번 중 6번 앞면이 나올 확률(가능도)이 최대가 된다면 그 $p$값으로 유추하는 것이 적절할 것이다. 따라서 가능도를 최대화시키는 과정이 필요하다. 따라서 $ p^6 (1-p)^4$값이 최대가 되어야 하며, 이는 $\frac {\partial L} {\partial p}=0$이 되는 $p$를 찾는 것과 동일하다.  따라서 $p= \frac {6} {10}$이 도출된다. 이것이 Maximum likelihood estimation이다. 

2. 로지스틱 회귀분석의 회귀식

 위 내용을 로지스틱 회귀분석에 적용해볼 것이다. 이전에 로지스틱 회귀분석의 회귀식을 먼저 알아야 한다. 회귀식은 다음과 같이 나타난다. 

$$log(odds) = \alpha + \sum_i \beta_i x_i$$

그런데, $odds= \frac {p} {1-p}$이므로 위 식은 다음과 같이 변형될 수 있다.

$$p= \frac {1} {1+ e^{-\alpha-\sum_i \beta_i x_i}}$$

즉 위 식은 어떤 개인에게 이벤트(고혈압 여부)가 존재할 확률을 $\alpha$, $\beta$, $x$로 표현할 수 있음을 의미한다.

3. 단순한 예시

계산을 단순화하기 위해서 위 $p= \frac {1} {1+ e^{-\alpha-\sum_i \beta_i x_i}}$식에서 $n=1$인 상황을 생각해보자. 즉 독립변수는 한 개만 존재한다고 생각한다. 그리고 그 독립변수는 이분형 변수라고 생각하자. 예를 들면 음주 여부 (음주자 vs 비음주자)가 있을 수 있다. $x_1$는 음주 여부를 의미하고, 0은 비음주자, 1은 음주자를 의미한다고 하자.

비음주자: $x_1 = 0$

음주자: $x_1 = 1$

그리고 $p$는 고혈압일 확률을 의미한다고 하자. 그러면 다음과 같은 분할표를 작성할 수 있고, 각각 해당하는 사람의 수는 $a$, $b$, $c$, $d$라고 하자.

4. 음주 여부에 따른 고혈압일 확률

어떤 누군가가 고혈압일 확률은  $p= \frac {1} {1+ e^{-\alpha- \beta_1 x_1}}$이다. 비음주자는 $x_1 = 0$, 음주자는 $x_1 = 1$이므로 다음을 구할 수 있다.

비음주자가 고혈압에 걸릴 확률: $ \frac {1} {1+ e^{-\alpha}}$

음주자가 고혈압에 걸릴 확률: $ \frac {1} {1+ e^{-\alpha-\beta_1}}$

비음주자 (혹은 음주자)가 고혈압에 안 걸릴 확률은 1에서 위 확률을 빼주면 되므로 다음을 구할 수 있다.

비음주자가 정상인일 확률: $1- \frac {1} {1+ e^{-\alpha}}$

음주자가 정상인일 확률: $1- \frac {1} {1+ e^{-\alpha-\beta_1}}$

5. 가능도 (likelihood) 계산하기

 위 분할표와 같이 딱 $a$, $b$, $c$, $d$이 위 표에 존재할 확률(가능도)은 얼마일까?

1) 먼저  음주자가 $a$명이 고혈압에 걸릴 확률을 계산해 보겠다.

음주자 1명이 고혈압에 걸릴 확률은 $ \frac {1} {1+ e^{-\alpha-\beta_1}}$이다.

음주자 $a$명이 고혈압에 걸릴 확률은 $ \frac {1} {1+ e^{-\alpha-\beta_1}}$을 $a$번 곱한 $\left( \frac {1} {1+ e^{-\alpha-\beta_1}}\right)^a$이다.

 (연구대상들이 고혈압에 걸리는 사건들은 모두 독립이라는 가정이 필요하다. 동전이 앞면이 나오면서 주사위에 1이 나올 확률을 구할 때 그저 $\frac{1}{2}$에 $\frac{1}{6}$을 곱할 수 있는 이유는 동전을 던지는 것과 주사위를 던지는 것이 서로 아무런 영향을 주지 않는 '독립'이기 때문이다. 음주자 $a$명을 다룰 때에도 $a$번 곱하여 확률을 구할 수 있다는 데에는 고혈압에 걸리는 사건들이 전제조건이 필요하다. 따라서 유전 정보를 공유하는 가족이나 쌍둥이 등이 연구대상에 있다면 가정을 만족하지 않을 수 있다.)

2) 음주자 $b$명이 고혈압에 걸리지 않을 확률은 $\left( 1- \frac {1} {1+ e^{-\alpha-\beta_1}} \right)$을 $b$번 곱한 $\left( 1- \frac {1} {1+ e^{-\alpha-\beta_1}}\right)^b$이다.

3) 비음주자 $c$명이 고혈압에 걸릴 확률은 $ \frac {1} {1+ e^{-\alpha}}$을 $c$번 곱한 $\left( \frac {1} {1+ e^{-\alpha}}\right)^c$이다. 

4) 비음주자 $d$명이 고혈압에 걸리지 않을 확률은 $\left( 1- \frac {1} {1+ e^{-\alpha}}\right) $을 $d$번 곱한 $\left( 1- \frac {1} {1+ e^{-\alpha}}\right)^d$이다. 

가능도

 그렇다면, 가능도($L$)는 위에서 계산한 값을 모두 곱하여 구할 수 있다.

$$L=\left( \frac {1} {1+ e^{-\alpha-\beta_1}}\right)^a  \left( 1- \frac {1} {1+ e^{-\alpha-\beta_1}}\right)^b \left( \frac {1} {1+ e^{-\alpha}}\right)^c \left( 1- \frac {1} {1+ e^{-\alpha}}\right)^d $$

6. 가능도 최대화 하기 : Maximum likelihood estimation

 가능도가 최대화되어 있다면, 가능도는 극대점에 있을 것이므로, 미분하였을 때 0이 될 것이다.

따라서 다음 식이 성립한다.

$$\frac {\partial L} {\partial \alpha}=0$$

$$\frac {\partial L} {\partial \beta_1}=0$$

위 두 식이 성립하는 $\alpha$와 $\beta_1$값을 구하는 것이 목적이다.

그런데 이 식은 계산을 하기가 성가시다. 대신에, $L$에 로그를 씌운 $LL=log(L)$에 대해 다음을 만족하는 $\alpha$와 $\beta_1$값을 구한다고 하자.

$$\frac {\partial LL} {\partial \alpha}=0$$

$$\frac {\partial LL} {\partial \beta_1}=0$$

아래 두 식으로 구한 $\alpha$와 $\beta_1$값과 위 두 식으로 구한 $\alpha$와 $\beta_1$값은 정확히 일치하는데, 아래 두 식은 계산하기가 훨씬 수월하다. 

연립 방정식을 통하면 다음을 알 수 있다.

$$\begin{align} \alpha&=\log\frac{c}{d} \\ \beta_1&=\log\frac{ad}{bc}=\log\left(OR\right)\\&&\end{align}$$

복잡한 식도 비슷하게 구하면 회귀계수들을 계산해낼 수 있다.

 선형 회귀분석은 손으로 직접 회귀계수를 구하는 방법들이 많이 언급되어 있던데, 로지스틱 회귀분석의 경우 언급된 경우가 많지 않아 본 포스팅을 작성해 보았다. 본 글에서 확인할 수 있듯이, 여기에서는 변수의 정규성 등을 가정하지 않았다. 따라서 선형 회귀분석과는 달리 로지스틱 회귀분석에서는 잔차의 정규성을 전제하지 않는다.

[이론] 로지스틱 회귀분석에서 회귀 계수를 구하는 방법 - Maximum likelihood estimation 정복 완료!

작성일: 2022.11.25.

최종 수정일: 2023.05.15.

[이론] 콕스 회귀 모델에서 생존 함수의 표현

콕스 회귀 모델에서 생존 함수는 다음과 같이 표현하곤 한다. 

$$S(t,\mathbf{X}) = \left[ S_0 (t) \right] ^{e^{\sum_{i=1}^{p} \beta_i X_i}} $$

왜 이렇게 표현되는지를 이해해야 추후에 나오는 log-log plot 등을 이해할 수 있기에 여기에서 소개한다.

1. Notation

먼저 notation에 대해 소개하고자 한다.

1)$\mathbf{X}$

$\mathbf{X}$는 $X$들의 집합체이며, $X$는 모델을 설명하는 변수를 의미한다. 예를 들면, 특정 위험 요소에 노출 여부, 성별, 나이, 인종, 음주 여부, 흡연 여부 등이 있을 수 있다. 이 모든 것들을 적기에는 귀찮고 공간 낭비이므로 $\mathbf{X}$로 표현한다.

2)$S(t,\mathbf{X})$

$S(t,\mathbf{X})$는 특정 $\mathbf{X}$을 가진 사람의 시간 $t$에서의 생존 확률이다. 예를 들면, 남자, 53세, 백인, 음주자, 비흡연자의 12개월에서의 생존 확률인 것이다.

3) $S_0(t)$

$S_0 (t)$는 모델을 설명하는 변수들의 값이 기본값(reference)인 사람의 시간 $t$에서의 생존 확률을 의미한다. 예를 들어 모델의 기본값을 남성, 0세, 흑인, 비음주자, 비흡연자로 잡았다면, 이런 사람의 12개월에서의 생존 확률을 의미한다.

4) $e^{\sum_{i=1}^{p} \beta_i X_i}$

이 식에서 모델을 설명하는 변수는 $p$개 였음을 알 수 있고, $\beta_i$는 각 변수에 대한 회귀계수를 의미한다.

2. 생존 함수와 위험 함수의 관계

이전 포스팅에서 다음과 같은 관계를 얻어냈다. 2022.11.10 - [통계 이론] - [이론] 생존 함수와 위험 함수의 관계

$$S(t) = \exp \left[ - \int_{0}^{t} h(u) du  \right] = e^{- \int_{0}^{t} h(u) du}$$

3. 콕스 회귀 모델에서 생존 함수의 표현

추가적인 notation을 두 개만 다루고 넘아가겠다.

1)$h_{\mathbf{X}}(t)$

특정 $\mathbf{X}$을 가진 사람의 시간 $t$에서의 위험 함수를 나타낸다.

2)1)$h_0(t)$

모델을 설명하는 변수들의 값이 기본값(reference)인 사람의 시간 $t$에서의 위험 함수를 나타낸다.

그렇다면 다음을 알 수 있다.

$$\begin{align} S(t,\mathbf{X})&=\exp\left[-\int_{0}^{t}h_{\mathbf{X}}(u)du\right]\\&=\exp\left[-\int_{0}^{t}h_{0}(u)e^{\sum_{i=1}^{p}\beta_{i}X_i}du\right]\\&=\exp\left[-\int_{0}^{t}h_{0}(u)du\times e^{\sum_{i=1}^{p}\beta_{i}X_i}\right] \\&=\left[S_{0}(t)\right]^{\exp\left(\sum_{i=1}^{p}\beta_{i}X_i \right)} \\&& \end{align}$$

[이론] 콕스 회귀 모델에서 생존 함수의 표현 정복 완료!

작성일: 2022.11.10.

최종 수정일: 2022.11.10.

생존 함수와 위험 함수의 관계

 생존 분석에 대해 공부를 하다 보면 생존 함수, 위험 함수가 나오게 되고, 그들의 관계를 나타내는 식이 등장한다.

1. 생존 함수와 위험 함수의 정의

먼저 생존 함수($S(t)$), 위험 함수($h(t)$)에 대해 설명하고 넘어가고자 한다.

1) 생존 함수

생존 함수 $S(t)$는 시간 $t$일 때, 살아있는 사람의 분율 (proportion)이다. 즉, 시간 $t$까지 살아남을 확률을 의미한다.

- 처음 ($t=0$)에는 모든 사람이 살아있으므로 $S(0)=1$이다.

- 시간이 무한히 흐르고 나면 모든 사람이 죽으므로 $S(\infty)=0$이다.

2) 위험 함수

위험 함수 $h(t)$는 그 개념이 조금 더 복잡하다. 거칠게 설명하면 다음과 같다.

"시간 $t$까지 살아남았을 때, 그 순간에 단위시간당 사망할 확률"

이 설명은 이해하기 어려우므로 예시를 들어 이해해보고자 한다.

(0) 상황

아침 9시에 200명으로 연구를 시작했다.

하지만 12시가 되었을 때 100명이 살아남았다.

(1) 1분간 관찰해보기

12시부터 12시 1분까지 10명이 죽었다고 하자. 즉, 사망할 확률은 10%(소수점으로 나타내면 0.1)이다

그렇다면, 60초동안 사망할 확률은 0.1이므로 위험 함수의 값은 다음과 같이 정해진다.

$$h(12시)= \frac {0.1} {60} = 0.00166667$$

그런데 위에서 위험 함수는 "그 순간에 단위시간당 사망할 확률"이라고 하였으므로 60초라는 간격은 너무 길다. 60초가 아니라 10초 동안 관찰해보자

(2) 10초간 관찰해보기

12시부터 12시 0분 10초까지 5명이 죽었다고 하자. 즉, 사망할 확률은 5%(소수점으로 나타내면 0.05)이다

그렇다면, 10초동안 사망할 확률은 0.05이므로 위험 함수의 값은 다음과 같이 정해진다.

$$h(12시)= \frac {0.05} {10} = 0.005$$

이렇게 시간 간격을 무한히 줄이기 위해 0에 가까워지는 극한을 적용하면 위험 함수의 값을 구할 수 있다.

$t$이상 살 확률을 시각 변수 $T$에 대해 $P(T\geq t)$라고 표현한다면, 위험 함수 $h(t)$는 다음과 같이 나타낼 수 있다.

$$h(t) = \lim_{\Delta t \rightarrow 0} \frac {P (t \leq T \leq t+ \Delta t \vert T \geq t)} {\Delta t} $$

2. 생존 함수와 위험 함수의 관계

보통 다음과 같은 관계로 나타난다.

$$h(t) = -\frac{S'(t)} {S(t)} = - \frac {dS(t)/dt} {S(t)}$$

또한 이는 일종의 미분방정식이므로 양변을 적분하고 정리하면 다음을 얻을 수 있다.

$$S(t) = \exp \left[ - \int_{0}^{t} h(u) du  \right] = e^{- \int_{0}^{t} h(u) du}$$

어떻게 이런 관계가 도출되는지 설명하고자 한다.

시간 $[0,1]$ 사이를 굉장히 큰 수 $n$에 대해 $n$개의 구간으로 나눈다.

처음 ($t=0$)에는 모든 사람이 살아있으므로 $S(0)=1$이다. 이번에는 $S \left( \frac {1} {n} \right)$을 구해볼 것이다.

그렇다면 우리가 관심을 갖는 시간의 구간은 $\left[ 0, \frac{1}{n} \right]$이다.

 위험 함수는 "그 순간에 단위시간당 사망할 확률"이므로 매 순간마다 다른 값을 가지는 것이 당연하지만, $n$이 매우 큰 수이므로 구간 $\left[ 0, \frac{1}{n} \right]$은 매우 짧은 찰나의 순간일 것이고, 그 사이에는 값이 거의 변하지 않는다고 봐도 무방하다. 이 구간의 $h(t)$는 $h \left( \frac {1} {n} \right)$으로 고정되어 있다고 가정하자.

 그런데 위험함수는 단위 시간당 사망할 확률이므로 기준이 구간 $[0,1]$이다. 그런데 우리가 관심을 갖는 구간은 $\left[ 0, \frac{1}{n} \right]$이므로 이 구간에서 사망할 확률은 $h \left( \frac {1} {n} \right)$에 $\frac{1}{n}$을 곱해주어야 한다.

$$\left[ 0, \frac{1}{n} \right]에~사망할~확률=h \left( \frac {1} {n} \right)\frac{1}{n}$$

이 구간에서 생존할 확률은 이 값을 1에서 빼주면 된다.

$$\left[ 0, \frac{1}{n} \right]에~생존할~확률=1-h \left( \frac {1} {n} \right)\frac{1}{n}$$

그렇다면 $S \left( \frac {1} {n} \right)$은 처음 ($t=0$)에 생존한 사람의 분율인 $S(0)$에 $\left( 1-h \left( \frac {1} {n} \right)\frac{1}{n} \right)$을 곱하여 구할 수 있다.

$$S \left( \frac {1} {n} \right) =S(0) \times \left( 1-h \left( \frac {1} {n} \right)\frac{1}{n} \right)=1-h \left( \frac {1} {n} \right)\frac{1}{n} $$

$\left[ 0, \frac{1}{n} \right]$에서의 생존할 확률을 구한 것과 같은 논리로 이를 일반화하면 $\left[ \frac{j-1}{n}, \frac{j}{n} \right]$에 생존할 확률은 다음과 같다.

$$ \left[ \frac{j-1}{n}, \frac{j}{n} \right]에~생존할~확률=1-h \left( \frac {j} {n} \right)\frac{1}{n} $$

한편 어떤 수 $k$에 대해  $\frac {k} {n}$ 까지 살아남을 확률을 구하는 것은

$\left[0,\frac{1}{n}\right]$에 생존한 사람이

$\left[\frac{1}{n},\frac{2}{n}\right]$에 생존한 사람이

$\left[\frac{2}{n},\frac{3}{n}\right]$에 생존한 사람이

$...$

$\left[\frac{k-1}{n},\frac{k}{n}\right]$에 생존할 확률을 구하는 것과 같다.

따라서 $\frac {k} {n}$ 까지 살아남을 확률인 $S \left( \frac{k}{n}\right)$은 다음과 같이 구할 수 있다.

$$\begin{align}S \left( \frac{k}{n}\right) &= S(0) \times\left(1-h \left( \frac {1} {n} \right)\frac{1}{n} \right)\times\left(1-h \left( \frac {2} {n} \right)\frac{1}{n} \right)\cdots\times\left(1-h \left( \frac {k} {n} \right)\frac{1}{n} \right)\\&=\prod_{j=1}^{k} \left( 1-h \left( \frac{j}{n} \right) \frac{1}{n} \right)\\&& \end{align}$$

그러므로 다음 등식이 성립한다.

$$S \left( \frac{k}{n} \right) \div S \left( \frac{k-1}{n} \right) = 1-h \left( \frac{k}{n} \right) \frac{1}{n}$$

이 식을 정리하면 다음을 얻게 되고

$$h \left( \frac{k}{n} \right) = - \frac{S \left( \frac{k}{n} \right) - S \left( \frac{k-1}{n} \right)}{ \frac{1}{n} S \left( \frac{k}{n} \right)}$$

$\lim_{ n \rightarrow \infty, \frac{k}{n} \rightarrow t}$의 극한을 씌우게 되면 다음 등식이 성립한다.

$$h(t) = -\frac{S'(t)} {S(t)} = - \frac {dS(t)/dt} {S(t)}$$

따라서 위 수식을 얻게 된다.

[이론] 생존 함수와 위험 함수의 관계 정복 완료!

작성일: 2022.11.10.

최종 수정일: 2022.11.10.

[이론] p-value에 관한 고찰

 P-value란 무엇인가?

 의학 및 보건학 논문을 읽다 보면 빠지지 않고 나오는 숫자가 p-value다. 연구를 하는 사람들도 그저 p-value는 0.05보다 작기만을 바라는 경향이 있다. 하지만, p-value에 관한 의미를 정확히 이해하지 못한다면 엉뚱한 결론을 짓는 잘못을 저지를 수도 있다. 따라서 본 포스팅에서는 p-value가 무엇인지 알아보고자 한다.

왜 다들 p-value를 이해하려 하지 않을까?

 논문을 쓰는 저자들도, 읽는 독자들도 p-value의 의미를 이해하지 않고 사용하는 경향이 있다. 이는 아마 p-value의 개념 자체가 꽤 복잡하기 때문일 것이다. 본 포스팅이 그 복잡한 내용으로 시작한다면 그들과 똑같은 짓을 하는 것일 테니, 학문적으로 복잡한 내용은 글 말미에 언급하도록 하겠다. 

실생활에 존재하는 p-value (One-tailed)

인식은 못하고 있겠지만, 여러분 모두 p-value가 무엇인지 아주 깊이 이해하고 있다. 다음 상황을 보자.

 누구든 저 여행객 A 입장이 되면 야바위꾼이 사기꾼이라고 생각할 것이다. 놀랍게도 이 짧은 이야기에 귀무가설, 대립 가설, p-value에 관한 내용이 모두 담겨있다. 다시 이야기로 돌아가 왜 여행객 A가 야바위꾼이 사기꾼이라고 생각하게 됐는지 그 사고 과정을 낱낱이 살펴보도록 하자. 여행객 A는 아마 다음과 같은 사고 과정을 거쳤을 것이다.

(야바위꾼의 호객행위를 들었을 때)

 - "저 동전은 공정한 동전이라 앞면이 나올 확률이 $1/2$, 뒷면이 나올 확률이 $1/2$일 거야."

 - "10번 던지면 5회쯤은 앞면이, 5회쯤은 뒷면이 나오겠지. 물론 5:5가 아니라 6:4 정도의 결과도 발생할 수는 있겠지. 운이 좋아 뒷면이 6회가 나온다면 내가 20,000원을 따겠구먼."

(동전을 던지고 난 뒤)

 - "어떻게 앞면이 9번이나 나올 수가 있어?! 이 정도의 일이 일어날 확률이 얼마나 된다고?!"

 - "분명 저 동전은 일반 동전이 아닐 거야. 앞 면이 더 잘 나오도록 어떤 조치를 취해뒀을 거야."

이 이야기의 귀무가설, 대립 가설, p-value은 다음과 같다.

귀무가설 : "저 동전은 공정한 동전이라 앞면이 나올 확률이 $1/2$, 뒷면이 나올 확률이 $1/2$일 거야."

대립 가설: "앞 면이 더 잘 나오도록 어떤 조치를 취해뒀을 거야" 

p-value: "이 정도의 일이 일어날 확률"

p-value

위 내용에서 "이 정도의 일이 일어날 확률"은 얼마나 될까?

또한,  "이 정도의 일"이란 무엇일까?

여행객 A는 앞면이 10회 나왔다면 더 화가 났을 것이다. 앞면이 8회 나온 상황에서도 화가 날 수 있지만, "이 정도의 일"까지는 아닌 것이다. 따라서 "이 정도의 일"이란 "앞면이 10회 나오는 일"과 "앞면이 9회 나오는 일"인 것이다. 이런 확률은 다음과 같이 계산된다.

$$\begin{align} P(X=9) &= _{10}C_1  \left( \frac {1} {2} \right)^{10} = \frac {10} {1024} \\ P(X=10) &= \left( \frac {1} {2} \right)^{10} = \frac {1} {1024} \\ \therefore P&= \frac {1+10} {1024} = \frac {11} {1024} \approx 0.0107 \\&& \end{align} $$

즉 여행객 A는 1.07%의 확률을 뚫고 일어난 일이 본인에게 발생했다는 것을 믿을 수 없어 "앞 면이 더 잘 나오도록 어떤 조치를 취해뒀을 거야"라는 생각을 하는 것이 합리적이라고 봤을 것이다.

실생활에 존재하는 p-value (Twe-tailed)

위의 상황은 "one-tailed"의 상황이었다. 이게 뭔지 몰라도 아래 "two-tailed"을 다루는 이야기를 보면 이해가 될 것이다.

이 상황에서도 여행객들은 화가 날 것이다. "이런 일"이 일어날 확률이 얼마 되지 않을 테니 사기를 쳤다고 확신할 것이다. 하지만 앞의 상황과는 "이런 일"의 정의가 조금 바뀌게 된다. 동전을 던지고 난 뒤 여행객들의 생각은 다음과 같을 것이다.

(동전을 던지고 난 뒤)

 - "어떻게 한쪽이 9번이나 나올 수가 있어?! 이 정도의 일이 일어날 확률이 얼마나 된다고?!"

 - "분명 저 동전은 일반 동전이 아닐 거야. 앞 면이든 뒷면이든 한쪽이 더 잘 나오도록 어떤 조치를 취해뒀을 거야."

앞의 상황(one-tailed)에서는 "앞면"이었던 것이 "한쪽"으로 바뀌었다. 누가 이기든 극단적인 양쪽(Two-tailed, two-sided)의 사건이 발생할 확률로 확인하는 것이 양측 검정이다. 

P-value의 학술적인 표현

P-value는 다음과 같이 기술할 수 있다. "귀무가설이 맞다고 할 때 이런 현상이 발생했을 확률"

이 정의는 좀 날 것 그대로니 학술적인 수정을 가하면 다음과 같다. "귀무가설이 맞는 데도 불구하고, 대립 가설을 선택했을 확률" 

보통 관련성이 없는 내용을 귀무가설로 설정하므로 "아무 관련성이 없는 게 학문적 진실인데, 관련성이 있다고 결론 내렸을 확률"이라고도 할 수 있다.

통계 검정으로의 적용

 P-value=0.01이라고 하자.

  (1) 귀무가설: 독립 변수(X)와 종속 변수(Y) 사이에 아무 관련성이 없다는 것이 사실이라고 했을 때 

  (2) p-value: 지금과 같은 현상이 일어날 확률은 0.01이므로 일어나기 어려운 일이다.

  (3) 대립 가설: 따라서 모종의 관련성이 있다고 하자.

아무 관련성이 없다는 말은 분석 방법에 따라 다르게 표현된다.

  (1) 로지스틱 회귀분석에서는 $\beta=0$ 혹은 $OR=1$으로 표현된다.

  (2) 선형 회귀분석에서는 $\beta=0$으로 표현된다.

  (3) Cox 회귀분석에서는 $\beta=0$ 혹은 $HR=1$으로 표현된다.

분석 방법에 따라 귀무가설이 표현되는 방법은 서로 다르다.

통계 검정으로의 적용 - 예외

다음 세 가지 분석은 다른 분석과는 다르게 귀무가설이 채택되기를 바라는 분석이다.

1) 정규성 검정 (Normality test): Shapiro-Wilk test, Kolmogorov-Smirnov test

2) Ordinal logistic regression의 proportional odds assumption 검정인 Score test

3) Cox regression의 proportional hazard assumption 검정인 Schoenfeld residual test

각 통계분석의 귀무가설은 다음과 같다.

 1) 데이터가 정규성을 따른다.

 2) Proportional odds assumption을 만족한다.

 3) Proportional hazards assumption을 만족한다.

즉 p-value가 0.05보다 작은 경우 귀무가설을 기각할 수 있으므로 정규성을 따르지 않는다거나, proportional odds/hazard를 만족하지 않는다고 할 수 있다. 하지만 p-value가 0.05보다 큰 경우 귀무가설을 기각할 수는 없다. 하지만 이 말이 귀무가설이 맞다는 말이 아니므로 정규성을 만족한다든가, proportional odds/hazard를 만족한다고 할 수는 없는 것이다. 이 점에 유의하며 p-value를 해석해야 한다. 따라서 우리가 통계 검정을 할 때에는 귀무가설이 무엇인지, 대립 가설이 무엇인지 항상 생각해야 한다.

[이론] P-value 정복 완료!

작성일: 2022.09.05.

최종 수정일: 2022.09.05.

[이론] 연속성을 수정한 카이 제곱 검정 (Chi-squared test with Yates's correction for continuity)

이 글을 읽기 전에 이전 포스팅을 읽고 오길 권한다.

2022.08.29 - [통계 이론] - [이론] 카이 제곱 검정과 피셔 정확 검정의 관계

범주형 자료를 분석할 때 카이 제곱 검정을 많이 사용하곤 한다. 하지만 카이 제곱 검정에는 치명적인 단점이 존재한다. 이는 이산 분포인 초기하 분포를 연속형 분포인 카이 제곱 분포에 근사하는 과정에서 발생한다. 

이산형 변수 = 소수점이 없는 변수

 이산형 변수는 쉽게 말해 소수점 아래 데이터가 없는 변수다.

 예를 들어 30살인 사람 11명이 모여있을 때 이 사람들의 평균 나이는 어떻게 될까? 모두 30살이므로 평균도 30이라고 할 수 있을까? 도대체 30살이라는 것이 무엇인가?

 30살인 사람은 30년 0일을 산 사람부터 30년 364일 23시간 59분 59초 9999.. 를 산 사람까지를 아우르는 말이다. 일상적으로 내가 $x$살 $y$개월을 살았어도 모두들 $x$살 살았다고 표현한다. 그런데 소수점 이하 숫자를 버린다면 과소평가 (underestimation)하게 된다.  가상의 11명의 나이를 반올림해보기도 하고 버림을 해보기도 하면 다음과 같다.

 반올림을 하면 10명 중 절반 정도는 "억울하게" 더 많은 나이를 갖게 되고 (과대평가), 나머지 절반 정도는 "운이 좋게" 더 어린 나이를 갖게 된다 (과소평가). 그래서 나이의 평균과 반올림 나이의 평균은 거의 일치한다. 하지만 버림을 한 나이는 누구나 "운이 좋게" 어린 나이를 갖게 되므로 (과소평가) 평균은 항상 참값보다 적게 나오게 된다. 

Yates의 의문점

 그렇다면 "Yates"가 제기한 의문은 이것이다.

"버림을 시행한 자료가 원자료를 대표한다고 할 수 있겠는가?"

해결

해결은 두 가지 논리로 설명해 보겠다. 결론은 같으나 함의가 다르니 찬찬히 살펴보길 바란다.

방법 (1) - 분할표는 버림이 된 상태다!

 흡연자이자 폐암 환자인 사람은 30명이었다고 하자. 하지만, 만약 전체 인구가 10,000명이었다면 어땠을까? 폐암 환자이면서 흡연자인 사람은 300명일 수도 있고, 301명, 302,... , 309명 일 수도 있었을 것이다. 그런데 1,000명만 뽑으면서 흡연자이자 폐암 환자인 사람은 30명이 뽑힐 수도 있고, 30.1명이 뽑힐 수도 있고, 30.2, 30.3,..., 30.9명이 뽑힐 수도 있었다. 이 상황에서 깔끔히 버림을 하고 30이라고 적자고 합의를 한 것이다. 

 그러면 폐암 환자이자 흡연자인 "30명"이라는 값은 $30\leq x<31$을 대표하는 값이다. 근사하기 전 분포인 초기하 분포는 이산 분포이므로 30 다음의 숫자는 31이다. 하지만 근사 시킨 카이 제곱 분포는 연속 분포이므로 30 다음의 숫자가 31이 아니다. 30부터 31 사이에는 수많은 숫자가 존재한다. 따라서 여기에는 $30\leq x <31$라는 숫자가 존재하므로 $30\leq x <31$의 숫자를 대표하는 숫자로 "30명"을 쓰는 것은 적절하지 않다. 그 중간값인 "30.5명"을 쓰는 것이 더욱 적절하다. 그리하여 카이 제곱 통계량 계산 식이 다음과 같이 바뀌는 것이다.

$$\chi^{2} _{Yates}=\sum_{ij} \frac{\left( \lvert a_{ij}-e_{ij}\rvert-0.5 \right)^{2}}{e_{ij}}$$

아마 Yates가 의도한 연속성 수정 (Correction for contunuity)이라는 말은 여기에서 기인한 것일 거다.

하지만 결론적으로 연속성 수정한 카이 제곱 검정은 현실적으로 쓰이고 있지 않은데, 아마 다음의 이유 때문일 것이다.

방법 (2) - 분할표는 반올림이 된 상태다!

 버림이 된 상황이라고 가정한 것은 순전히 위 상황을 합리화하기 위한 것이었다. 세상 대부분의 일은 반올림이 시행된다.   예를 들어 소수점 첫째 자리까지 표기할 수 있는 디지털 체중계가 75.3kg을 나타냈다고 하자. 이 사람의 실제 몸무게는 $75.25\leq 몸무게<75.35$ 사이에 있다고 보는 것이 적당할 것이다.

 분할표에 있는 숫자 또한 반올림이 된 상태라는 가정 하에 다시 분할표를 살펴보자. 폐암 환자이자 흡연자인 사람의 수 "30명"이 나타내는 것은 $29.5\leq x<30.5$을 나타내는 숫자다. 그런데 "30명"이라는 숫자로 구한 카이 제곱 검정의 p-value는 30명이 발생하는 것보다 일어나기 어려운 사건들이 일어날 확률을 의미한다. 여기에서 기대 빈도는 60명이므로 60 이하에서는 관찰 빈도가 낮아진다는 것은 일어나기 어려운 상황을 의미한다. 따라서 $29.5\leq x \leg30$이 발생할 확률은 "30명"이라는 숫자로 구한 카이 제곱 검정의 p-value에 포함되어 있다. 하지만 $30<x<30.5$는 포함되지 못한다. 그런데 위에서 구한 식 

$$\chi^{2} _{Yates}=\sum_{ij} \frac{\left( \lvert a_{ij}-e_{ij}\rvert-0.5 \right)^{2}}{e_{ij}}$$

을 쓴다는 것은 다음을 의미한다.

"30이라는 숫자에는 $30<x<30.5$가 발생할 가능성도 내포하고 있는 것이다. 혹시 모르니까 30이 아닌 30.5로 계산하자."

따라서 p-value가 더 커지게 되고 보수적인 계산을 하게 되는 것이다.

그런데 왜 절댓값 기호가 붙었는지 이해하기 위해 다음 내용을 봐보자.

여기에서 2*2 분할표가 위와 같을 때 이 데이터로 구하는 카이 제곱 통계량의 의미를 다시 생각해보자.

 - 폐암 환자이면서 흡연자의 기대 빈도는 60명이다. 

 - 폐암 환자이면서 흡연자인 사람이 30명인 사건이 발생했다.

 - p-value는 "이 정도의 일이 발생할 확률이 얼마인가?"를 의미하는 것이다.

따라서, p-value를 구하기 위해서는 30명이 발생할 확률만 알아서는 안 된다.

 - 30명이 발생하는 것보다 일어나기 어려운 사건들의 확률을 모두 합해야 한다.

 - 기대 빈도인 60명이 발생할 확률이 가장 크고, 60에서 멀어질수록 확률은 낮아진다.

 (1) 30명보다 작은 경우 

   - 29명이 발생할 확률은 더 낮을 것이므로 더해야 한다.

   - 28, 27,... , 1, 0명이 발생할 확률은 더 낮으므로 모두 더해야 한다.

 (2) 30명보다 큰 경우 

   - 30명보다 커지는 상황을 상상해볼 때 60명이 될 때까지는 확률이 증가하다가 그 이후로 감소한다.

   - 하필이면 30명일 확률과 정확히 똑같아지는 시점은 90명일 때다. (계산을 해보면 알 수 있다.)

   - 따라서 90명일 발생할 확률, 91, 92,..., 199, 200명이 발생할 확률을 모두 더해야 한다.

"(1) 30명보다 작은 경우" 30이 아닌 30.5를 쓰는 이유는 위에서 모두 설명이 되었다,

"(2) 30명보다 큰 경우"에 대해 설명을 하겠다.

예를 들어 90이라는 숫자에는 $89.5\leq x<90.5$가 포함되어 있는 것이다. 그런데 $90\leq x<90.5$는 이미 카이 제곱 검정 계산에 포함되어 있다. 하지만 $89.5\leq x<90$은 포함되어 있지 않으므로 "89.5"를 기준으로 계산해야 한다는 것이 Yates의 논리다. 그가 의도했든 의도하지 않았든 말이다. 즉, 기댓값보다 작을 때에는 0.5를 더하고, 기댓값보다 클 때에는 0.5를 빼는 로직을 만들어야 하는데, 이를 수식으로 한 번에 나타내는 방법이 절댓값 기호를 붙이는 것이다.

하지만, 결국 Yates의 연속성 수정은 너무나 보수적이고, 학자들 사이에서도 의견이 분분하여 굳이 쓰이지 않고 있는 분석 방법이다.

[이론] 연속성을 수정한 카이 제곱 검정 (Chi-squared test with Yates's correction for continuity) 정복 완료!

작성일: 2022.08.30.

최종 수정일: 2022.08.30.

[이론] 카이 제곱 검정과 피셔 정확 검정의 관계

 범주형 자료 분석을 할 때 "기대 빈도가 5 미만인 셀이 25% 이상인 경우 카이 제곱 검정을 신뢰할 수 없으며 피셔 정확 검정의 결과를 확인해야 한다."라는 말을 정말 많이 보게 된다. 도대체 카이 제곱 검정과 피셔 정확 검정이 무슨 관계이길래 이렇다는 건지 궁금증이 유발될 것이다. 모든 것을 설명할 수는 없지만, 왜 이런 이야기들이 나오는지 대략적으로 설명하고자 한다.

본 글을 읽기 전에 카이 제곱 검정과 피셔 정확 검정의 이론 내용에 관한 포스트를 읽고 오기를 강력히 권한다.

카이 제곱 검정: 2022.08.16 - [통계 이론] - [이론] 카이 제곱 검정 (Chi-squared test)

피셔 정확 검정: 2022.08.26 - [통계 이론] - [이론] 피셔 정확 검정 (Fisher's exact test)

카이 제곱 검정=근사, 피셔 정확 검정=정확

먼저 이해해야 하는 것은 "왜 카이 제곱 검정을 쓰게 되었는가?"일 것이다. 결론적으로는 정확한 방법인 피셔 정확 검정 시 계산량이 너무 방대하여 그에 근사하는 카이 제곱 분포를 사용했다는 것이다.

 물론 Pearson 경이 어떤 생각으로 카이 제곱 검정을 만들었는지 정확하게 알 수는 없을 뿐 아니라, 그나마 세간에 알려진 이유도 나에게는 그다지 와닿지는 않는다. 피어슨 경의 생각을 읽어보면 아마도 다음과 같을 것이다. (따라서 본 포스팅은 문헌을 리뷰하거나 참고한 것이 아니며 피어슨 입장에서 생각해본 '뇌피셜'이다.)

분포가 무엇인지 대충이라도 상상해보기

흡연과 폐암의 빈도 표가 다음과 같다고 하자.

각 셀 안의 값은 정확히 모르지만 빨간색으로 표시된 총합의 값은 정확히 알려져 있다고 하자. 이때 (A)에 들어갈 수 있는 숫자는 이론적으로 0부터 200 사이의 값이다. 흡연과 폐암에 아무런 관련성이 없을 때 (A)에 들어갈 것으로 생각되는 숫자를 생각해보자. 폐암은 흡연과 관련이 없고, 전체 인구 중 흡연자는 30%를 차지하므로 폐암 환자 200명 중 30%인 60명이 있을 것으로 기대된다. 당연하게도 60명이 있을 확률이 가장 높아야 하고, 60보다 커지거나 작아질수록 그 확률은 감소해야 한다. 따라서 확률 분포 함수는 다음과 같을 것이다.

대립 가설 검정하기 <흡연자일수록 폐암환자일 가능성이 높아지는 관련성이 있다.>  (One-sided, One-Tailed)

실제로 관찰 빈도가 다음과 같았다고 하자.

 아무런 관련성이 없었다면 폐암 환자이면서 흡연자일 것으로 예상되는 사람의 수는 60명이었는데 75명이 관찰되었으므로 흡연자일수록 폐암환자일 가능성이 높아질 것이라고 예상해볼 수 있다. 이것이 통계적으로 유의미한지 확인하는 사고 과정은 다음과 같다.

"흡연자이면서 폐암환자로 예상된 사람은 60명이었는데 75명이나 있네?"

$\rightarrow$ "이 말은 폐암과 흡연이 관련성이 있다는 것 아닐까?"

$\rightarrow$  [가정] "흡연이랑 폐암 간에 아무런 관계가 없다고 가정해보자."

$\rightarrow$ "그런 가정 하에서 흡연자이면서 폐암환자인 사람이 75명 이상일 확률이 얼마인지 구해보자."

계산 방법: 아래 수식을 모두 합한다.

$$75명일 확률 = \frac{\begin{pmatrix} 300 \\  75 \end{pmatrix} \times \begin{pmatrix} 700 \\  125 \end{pmatrix}} {\begin{pmatrix} 1000 \\  200 \end{pmatrix}}$$

$$76명일 확률 = \frac{\begin{pmatrix} 300 \\  76 \end{pmatrix} \times \begin{pmatrix} 700 \\  124 \end{pmatrix}} {\begin{pmatrix} 1000 \\  200 \end{pmatrix}}$$

$$77명일 확률 = \frac{\begin{pmatrix} 300 \\  77 \end{pmatrix} \times \begin{pmatrix} 700 \\  123 \end{pmatrix}} {\begin{pmatrix} 1000 \\  200 \end{pmatrix}}$$

$$\vdots$$

$$200명일 확률 = \frac{\begin{pmatrix} 300 \\  200 \end{pmatrix} \times \begin{pmatrix} 700 \\  0 \end{pmatrix}} {\begin{pmatrix} 1000 \\  200 \end{pmatrix}}$$

위 값들을 모두 더하면 0.00674213이다.

$\rightarrow$ "이 정도의 확률은 현실적에서 일어나기는 어려운 일 아니야?"

$\rightarrow$ "그러면 차라리 흡연자일수록 폐암에 걸릴 확률이 높은 게 진실이고 그런 현실 속에서 발생한 일이라고 보는 게 낫겠다."

위 사고 구조의 비현실성: 방대한 계산량

 위 사고 과정은 틀린 것이 없고 논리적이기만 하다. 하지만 피어슨(Pearson) 경이 활동했을 1900년 경에는 복잡한 계산을 순식간에 해줄 계산기가 없었기 때문에, "75명일 확률", "76명일 확률", "77명일 확률",... , "200명일 확률"을 일일이 계산한다는 것은 굉장히 귀찮고 시간이 많이 걸리지만 그에 비해 그만한 가치는 별로 없는 일이었을 것이다.

 예를 들어 "75명일 확률"만 하더라도 

$$75명일 확률 = \frac{\begin{pmatrix} 300 \\  75 \end{pmatrix} \times \begin{pmatrix} 700 \\  125 \end{pmatrix}} {\begin{pmatrix} 1000 \\  200 \end{pmatrix}}= \frac {\frac {300!} {75! \times 225!}\times \frac {700!} {125! \times 575!}} { \frac {1000!} {200! \times 800!}} $$

와 같은데, 저 계산을 한다는 것은 아주 매우 많이 귀찮은 일일 것이다.

근사하는 분포를 만들자

 그래서 피어슨 경은 초기하 분포 (hypergeometric distribution)를 계산이 이미 되어있는 어떤 분포에 근사 시키고자 하는 욕구가 들었을 것이다. 아마도 피어슨 경은 두 번의 근사를 통해 초기하 분포를 정규분포에 근사 시키고자 했을 것이다.

 1) 초기하 분포 (hypergeometric distribution)에서 이항 분포(binomial distribution)로

 2) 이항 분포 (binomial distribution)에서 정규 분포 (normal distribution)로

 1) 초기하 분포 (hypergeometric distribution)에서 이항 분포(binomial distribution)로

위와 같은 상황에서 다음 조건을 만족하면 초기하 분포는 이항 분포로 근사 시킬 수 있다.

$$K>0.1 \times N$$

 2) 이항 분포 (binomial distribution)에서 정규 분포 (normal distribution)로

어떤 이항 분포의 시행 횟수가 $n$, 발생 확률이 $p$일 때, $np>5, n(1-p)>5$이면 이항 분포는 정규 분포로 근사될 수 있다. 

여기에서 $np$와 $n(1-p)$는 각각 사건이 일어나 기댓값과 일어나지 않을 기댓값을 의미한다. "기대 빈도가 5 미만인 셀이 25% 이상인 경우 카이 제곱 검정을 신뢰할 수 없으며 피셔 정확 검정의 결과를 확인해야 한다."라는 말은 여기에서 기인한 것으로 보인다.

근사하는 정규분포를 지정하자

 정규 분포의 확률 밀도 함수 (probability density function)은 다음과 같다.

$$ {\frac {1}{\sigma {\sqrt {2\pi }}}}e^{-{\frac {1}{2}}\left({\frac {x-\mu }{\sigma }}\right)^{2}}$$

갑자기 수식 이야기를 왜 하느냐? 

 수식을 보면 정규 분포를 정의하는 데에 오직 평균($\mu$)과 표준편차($\sigma$)만이 필요하다는 것을 알 수 있다. 즉, 초기하 분포의 평균과 표준편차 (혹은 분산)를 알아내면 근사하는 정규분포를 지정할 수 있다.

(1) 평균

초기하 분포 (hypergeometric distribution)의 이산 밀도 함수 (discrete density function)은 $$f_{k} (k;K,N,n)= \frac{\begin{pmatrix} K \\  k \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\begin{pmatrix} N \\  n \end{pmatrix}}$$ 이므로 평균($\mu$)은 다음과 같이 구할 수 있다.

$$\begin{align}  \mu &=\sum _{k} {k \times \frac{\begin{pmatrix} K \\  k \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\begin{pmatrix} N \\  n \end{pmatrix}} } \\&=\sum _{k} {K \times \frac{\begin{pmatrix} K-1 \\  k-1 \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\begin{pmatrix} N \\  n \end{pmatrix}} } \\&=\sum _{k} {K \times \frac{\begin{pmatrix} K-1 \\  k-1 \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\frac{N}{n} \times \begin{pmatrix} N-1 \\  n-1 \end{pmatrix}} } \\&= \frac {Kn} {N} \\&&\end{align} $$

(2) 분산

분산은 다음과 같이 계산할 수 있다.

$$ \begin{align} \sigma^{2}&=\sum _{k} {k^{2} \times \frac{\begin{pmatrix} K \\  k \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\begin{pmatrix} N \\  n \end{pmatrix}} }  - \mu^{2} \\ &=\sum _{k} {Kk\times \frac{\begin{pmatrix} K-1 \\  k-1 \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\begin{pmatrix} N \\  n \end{pmatrix}} }  - \left( \frac {Kn} {N} \right)^{2} \\& =K \left( \sum _{k} { \left(k-1\right)\times \frac{\begin{pmatrix} K-1 \\  k-1 \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\begin{pmatrix} N \\  n \end{pmatrix}} }  + \sum _{k} { \times \frac{\begin{pmatrix} K-1 \\  k-1 \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\begin{pmatrix} N \\  n \end{pmatrix}} } \right)  - \left( \frac {Kn} {N} \right)^{2} \\&=K \left( \sum _{k} { \left(K-1\right)\times \frac{\begin{pmatrix} K-2 \\  k-2 \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} { \frac{N \left( N-1 \right)}{ n \left( n-1\right)} \begin{pmatrix} N-2 \\  n-2 \end{pmatrix}} }  + \sum _{k} { \times \frac{\begin{pmatrix} K-1 \\  k-1 \end{pmatrix} \times \begin{pmatrix} N-K \\  n-k \end{pmatrix}} {\frac{N} {n} \begin{pmatrix} N-1 \\  n-1 \end{pmatrix}} } \right)  - \left( \frac {Kn} {N} \right)^{2}\\&=\frac {K(K-1)n(n-1)} {N(N-1)} + \frac {Kn} {N} - \left( \frac {Kn} {N} \right)^{2}\\&= n \cdot \frac {K} {N} \frac {N-K} {N} \frac {N-n} {N-1} \\&& \end{align} $$

다른 방법으로는 factorial moment를 사용하여 $[E(X(X-1))]$을 구하는 방법도 있다. 수식이 조금 더 간단해진다.

따라서 흡연자이면서 폐암 환자인 사람의 수는 평균이 $\frac {Kn} {N} $, 분산이 $n \cdot \frac {K} {N} \frac {N-K} {N} \frac {N-n} {N-1} $인 정규분포를 따른다고 할 수 있다. 

정규 분포가 아닌 카이 제곱 분포를 사용하자

피어슨 경은 이 정규 분포를 바로 사용하는 것이 아니라 카이 제곱 분포를 사용하고자 했다. 이런 2*2 분할표 (contingency table)에서는 정규 분포를 쓰든 카이 제곱 분포를 쓰든 아무 상관이 없겠지만 표가 더 커지면 문제가 발생하기 마련이다.

<2*2 분할표>

<3*3 분할표>

 2*2 분할표에서는 자유도가 1이므로 한 개의 값만 지정하면 되지만, 3*3 분할표만 되어도 자유도가 4이므로 자유롭게 정할 수 있는 값이 4개가 된다. 이런 경우 정규 분포를 바로 사용할 수 없다. 따라서 정규분포를 사용하지만 자유도의 개념이 있는 카이 제곱 분포를 사용하고자 했을 것이다.

 설명을 용이하게 하기 위해 2*2 분할표로 설명을 이어가도록 하겠다. 흡연자이면서 폐암 환자인 사람의 수($X$)를 표준화한 뒤 제곱해주면 카이 제곱 분포를 따른다고 할 수 있다. 평균이 $\frac {Kn} {N} $, 분산이 $n \cdot \frac {K} {N} \frac {N-K} {N} \frac {N-n} {N-1} $이므로 

$$\frac{\left(X-\frac {Kn} {N}\right)^2}{  n \cdot \frac {K} {N} \frac {N-K} {N} \frac {N-n} {N-1} }\sim \chi^2(1) \tag{1}$$

이라고 작성할 수 있다.

 이때 $N$이 굉장히 크면 $N\sim \left(N-1\right)$이므로 $(1)$식은 $$\frac{\left(X-\frac {Kn} {N}\right)^2}{  n \cdot \frac {K \cdot (N-K) \cdot(N-n)} {N^3}}\sim \chi^2(1) \tag{2}$$

에 근사 시킬 수 있다.

수식을 정리해보자

 한편 수식$(2)$에서 쓰인 문자들은 표에서 다음과 같이 나타난다.

수식$(2)$을 어떻게 잘 정리해야 세상에 잘 먹힐지 피어슨 경은 고민이 많았을 것이다. 저 식은 너무 복잡해서 일반적인 연구자가 쓰기엔 복잡할 뿐만 아니라 쓰고 싶지 않게 생겼기 때문이다. 피어슨 경은 뛰어난 직관으로 해결했을 수도 있지만 나에게는 다음과 같이 계산되는 결과가 매력적으로 느껴졌다. 

먼저 기대 빈도를 구하면 다음과 같다.

관찰 빈도에서 기대 빈도를 뺀 뒤 제곱한 값은 놀랍게도 모든 셀에서 같다. (사실 자유도가 1이므로 당연한 현상이긴 하다.)

게다가 이 값은 수식$(2)$인 $\frac{\left(X-\frac {Kn} {N}\right)^2}{  n \cdot \frac {K \cdot (N-K) \cdot(N-n)} {N^3}}$의 분자에 해당하는 내용이다. 따라서 이 수식은 다음과 같이 바꿔보고 싶은 욕구가 차오른다.

$$ \begin{align} \frac{\left(X-\frac {Kn} {N}\right)^2}{  n \cdot \frac {K \cdot (N-K) \cdot(N-n)} {N^3}}&=\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} \frac{\frac{Kn}{N}}{  n \cdot \frac {K \cdot (N-K) \cdot(N-n)} {N^3}}\\&= \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} \left(  1+ \frac{N^2}{(N-K)(N-n)}-1\right)\\&= \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} + \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} \frac{KN+Nn-Kn}{(N-K)(N-n)}\\&=\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} + \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{K(N-n)}{N}} \frac{KN+Nn-Kn}{(N-K)n}\\&=\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} + \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{K(N-n)}{N}} \left( 1 + \frac{KN} {(N-K)n} \right)\\&=\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} + \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{K(N-n)}{N}}  +\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{K(N-n)}{N}}  \frac{KN} {(N-K)n} \\&=\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} + \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{K(N-n)}{N}}  +\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{(N-K)n}{N}}  \frac{N} {(N-n)} \\&= \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} + \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{K(N-n)}{N}}  +\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{(N-K)n}{N}} \left( 1+  \frac{n} {(N-n)} \right)\\&=\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{Kn}{N}} + \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{K(N-n)}{N}}  +\frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{(N-K)n}{N}} +  \frac{\left(X-\frac {Kn} {N}\right)^2}{\frac{(N-K)(N-n)}{N}} \\&= \frac  {\left(a_{11}-e_{11}\right)^2} {e_{11}} +\frac  {\left(a_{21}-e_{21}\right)^2} {e_{21}} +\frac  {\left(a_{12}-e_{12}\right)^2} {e_{12}}+\frac  {\left(a_{22}-e_{22}\right)^2} {e_{22}} \\&= \sum_{i, j} \frac  {\left(a_{ij}-e_{ij}\right)^2} {e_{ij}} \tag{3} \\&& \end{align}$$

정리된 수식의 위엄

 이제 정확한 초기하 분포의 확률을 구하기 위해 $1000!$같은 무식한 계산을 하지 않아도 된다. $2*2$ 분할표라면 관찰 빈도와 기대 빈도로 계산한 $(3)$ 식의 값이 3.84를 넘기만 하다면 분포에 통계적으로 유의한 차이가 있다고 양측 검정 (Two-tailed or two-sided)을 한 셈이니 말이다. (3.84는 자유도가 1인 카이 제곱 분포의 누적 확률이 0.95가 되는 지점이다.) $m*n$ 분할표일 때에도 식 $(3)$과 비슷하게 계산을 먼저 하고 그 값이 카이 제곱 분포표에서 자유도 $(m-1)\times(n-1)$일 때 $\alpha=0.05$인 지점보다 큰지만 확인하면 된다.

[이론] 카이 제곱 검정과 피셔 정확 검정의 관계 정복 완료!

작성일: 2022.08.29.

최종 수정일: 2022.08.29.