11 일반화 선형모형

주의

이 장에서는 벡터와 스칼라의 표기를 구분없이 모두 보통체로 사용한다.

11.1 일반화 선형모형

선형모형의 의미는 모형에서 고려하는 설명변수가 변할 때 반응값의 평균이 변하는 관계가 선형이라는 것이다. 즉, 반응값의 평균을 설명하는 회귀식이 회귀계수에 대하여 선형이라는 의미이다.

\[ E(y|x_1,x_2,\dots,x_p) = \beta_ 1 x_1 + \dots + \beta_p x_p \equiv \eta \tag{11.1}\]

참고로 식 11.1 의 오른쪽에 나타나는 식을 선형예측식(linear predictor, \(\eta\))라고 부른다.

이러한 평균의 선형성의 가정이 적절하지 않은 경우가 있다. 예를 들어 공학이나 생물학에서 사용되는 비선형 회귀모형(nonlinear regression model)처럼 반응변수 평균의 변화가 설명변수들의 복잡한 비선형 관계(예를 들어 미분방정식의 관계)로 나타나는 경우로 흔히 나타난다. 이러한 비선형 회귀모형은 섹션 10.4 장에서 다루었다.

반응변수가 가질 수 있는 평균값의 범위에 제한이 있을 수 있다. 예를 들어 베르누이 분포의 경우 평균이 성공확률이기 때문에 0과 1사이에 있으며 포아송 분포의 경우 반응값은 음의 값을 가질 수 없다. 따라서 식 11.1 의 선형예측식 \(\eta\)와 반응값의 평균 \(E(y| x)\)의 관계를 선형모형 식 11.1 처럼 정의할 수 없다.

이렇게 반응변수의 평균과 선형에측식의 범위가 일치하지 않는 경우 임의의 단조증가 함수 \(g\)를 사용하여 그 범위를 일치하게 만들어 줄 수 있다. 예를 들어 베르누이 분포의 경우 표준 정규분포의 누적분포함수 \(\Phi\)를 사용하여 확률의 범위와 선형에측식의 범위를 맞추어 줄 수 있다. 이러한 회귀모형을 프로빗(probit)모형이라고 부른다.

\[ \Phi^{-1} [p(y| x)] = \beta_ 1 x_1 + \dots + \beta_p x_p =\eta \tag{11.2}\]

이제부터 정규분포하에서 정의되는 선형모형을 다른 분포들로 확장한 모형인 일반화 선형모형(Generalized Linear Model; GLM)을 살펴보기로 하자.

11.1.1 지수군 분포와 일반화 선형모형

지수군 분포에 대한 자세한 내용은 부록 H 에 있으니 참고하자.

일변량 확률변수 지수군(exponential family) 분포를 따른다고 가정하자.

\[ f (y | \theta, \phi ) = \exp \left \{ \frac{y \theta-b (\theta) }{a(\phi) } + c(y,\phi) \right \} \]

충분통계량이 관측값 \(y\) 이고 1차원의 기본형 모수 \(\theta\) 로 정의된 분포임을 유의하자.

확률변수 \(y\)의 평균을 \(\mu = E(y)\) 이라고 하고 독립변수 벡터 \(x\)와 회귀계수 벡터 \(\beta\)로 구성된 선형예측식을 \(\eta\)라고 하자.

\[ \eta = { x}^t \beta \tag{11.3}\]

일반화 선형모형은 분포의 특성에 따라 주어진 단조증가함수 \(g\) 를 이용하여 \(y\) 의 평균과 선형예측식의 관계를 설정하는 모형이다. 이러한 함수 \(g\) 를 연결함수(link function)라고 부른다.

\[ g(\mu) = g(E[y| x]) = { x}^t \beta \tag{11.4}\]

반응값의 분포가 주어진 경우 연결함수는 평균의 범위와 선형예측식의 범위를 연속적으로 1-1 대응하게 해주는 함수이며 사용할 수 있는 가능한 함수는 무한히 많다. 예를 들어 베르누이 분포의 경우 위에서 정의된 프로빗 모형 식 11.2 에서 \(\Phi^{-1}\) 같이 (0,1) 에서 실수전체 집합으로 단조증가하는 함수는 모두 연결함수로 고려할 수 있다.

지수군 분포에서 모수 \(\theta\)를 기본형모수(canonical parameter)라고 부르며 일반적으로 \(\theta\) 는 평균 \(\mu\)의 비선형 함수로 나타난다. 만약 다음과 같은 관계를 나타내는 연결함수 \(g\) 가 있다면 그 함수를 기본형 연결함수(canonical link function) 또는 자연결합함수(natural link function)라고 부른다.

\[ \theta = \eta \tag{11.5}\]

예를 들어 이항분포의 확률밀도함수를 지수군 분포의 형태로 표현했을 때 그 형태를 보면 다음과 같은 관계를 알 수 있다.

\[ \theta = \log \frac{\mu}{1-\mu} = \log \frac{p}{1-p} \]

따라서 식 11.5 을 만족하는 연결함수는 다음과 같고

\[ \log \frac{p}{1-p} = \eta = { x}^t \beta \tag{11.6}\]

이를 로짓 연결함수(logit link function)이라고 부르며 이는 이항분포의 기본형 연결함수이다.

노트

만약 \(y\)의 분포가 정규분포이며 연결함수 \(g\)가 \(g(\mu) = \mu\)이면 선형회귀모형이 된다. \[ E[y| x] = { x}^t \beta \]

11.2 일반화 선형모형의 가능도함수

하나의 표본 \(y\)에 대하여 기본형 모수 \(\theta\) 하나인 로그가능도함수(log likelihood function) \(\ell\)은 다음과 같이 정의된다.

\[ \ell = \log f(y) = \frac{y\theta-b(\theta)}{a(\phi) } + \log c(y, \phi) \]

표본 \(y_1,y_2,\dots,y_n\) 가 각각 설명변수 벡터\({ x}_1, { x}_2, \dots,{ x}_n\)에서 독립적으로 얻어졌다면 로그가능도함수 \(\ell_n\) 은 다음과 같다.

\[ \ell_n = \log \prod_{i=1}^n f(y_i) = \sum_{i=1}^n \frac{y_i \theta_i-b(\theta_i)}{a(\phi_i) } + \sum_{i=1}^n \log c(y_i, \phi) \tag{11.7}\]

여기서 \(\theta_i\)는 \(i\) 번째 관측값에 대한 모수로서 첨자 \(i\) 를 붙이는 이유는 관측치의 기대값 \(\mu_i = E(y_i | { x}_i)\)가 독립변수의 값 \({ x}_i\)에 따라 다를 수 있고 \(\theta_i\)는 평균 \(\mu_i\)의 함수이기 떄문이다.

이제 식 11.4 과 같이 설명변수와 반응변수 평균과의 관계가 연결함수 \(g\)로 정의되었다고 하자.

\[ g(\mu_i) = g(E[y_i| { x}_i ]) = { x}_i^t \beta \equiv \eta_i, \quad i=1,2,\dots,n \tag{11.8}\]

보기 11.1 (이항분포) 주어진 예측변수 \(\pmb x_i=(x_{1i},\dots,x_{pi})^t\)에서 실행횟수가 \(m_i\)인 이항분포\(B(m_i, p(x_i))\)를 생각하자. \(m_i\)의 시행 중에 성공의 횟수가 \(y_i\)라고 하면 \(y_i\)의 평균과 분산은 다음과 같다.

\[ E(y_i | x_i ) = m_i p(x_i), \quad \quad Var(y_i | x_i) = m_i p(x_i) (1-p(x_i)) , \quad i=1,2,\dots,n \]

이항분포를 위한 로지스틱 회귀방정식은 선형예측식과 성공의 확률의 관계를 다음과 같이 정한다.

\[ \log \left [ \frac{p(x_i)}{1-p(x_i)} \right ] = {\pmb x}^t_i {\pmb \beta} = \beta_0 + \beta_1 x_{1i} + \dots \beta_p x_{pi} \tag{11.9}\]

서로 독립인 관측값 \((y_1,y_2,\dots,y_n)\)의 가능도함수(likelihood function) \(L\)은 이항분포들의 결합확률밀도함수와 같고 아래와 같이 주어지며

\[ L = \prod_{i=1}^n f(y_i|p(x_i)) = \prod_{i=1}^n \left [ {{m_i}\choose{y_i}} \right] {p(x_i)}^y_i {(1-p(x_i))}^{m-y_i} \]

로그가능도함수(log likelihood function) \(l\) 은 다음과 같이 나타낼 수 있다.

\[ \begin{aligned} l & = \log L = \sum_i \log {{m_i}\choose{y_i}} + \sum_i y_i \log p(x_i) + \sum_i (m_i -y_i) \log (1-p(x_i)) \notag \\ & = c(\pmb y,\pmb m) + \sum_i y_i \log \left [ \frac{p(x_i)}{1-p(x_i)} \right ] + \sum_i m_i \log (1-p(x_i)) \end{aligned} \tag{11.10}\]

위의 로그가능도함수에서 볼 수 있듯이 충분통계량인 성공의 횟수 \(y_i\)와 곱으로 나타내어진 함수가 로짓함수이며 이렇게 가능도함수에서 얻어진 결합함수를 기본형 연결함수라고 한다.

노트

Nelder 와/과 Wedderburn (1972) 에서 연결 함수(link function)의 개념이 제시할 때 다음과 같은 작업 변량(working variate) \(z_i\)를 이용하여 선형모형을 일반화하고자 하였다. 작업 변량 \(z_i\) 은 다음과 같이 정의된다.

\[ \begin{aligned} z_i & = g(y_i) \\ & \simeq g(\mu_i ) + g_{\mu}(\mu_i)(y_i - \mu_i) \\ & = { x}_i^t \beta + g_{\mu}(\mu_i)(y_i - \mu_i) \\ & \simeq { x}_i^t \beta + r_i \end{aligned} \]

위의 식에서 \(g_{\mu}\) 은 연결함수 \(g\) 를 미분한 함수이다. 위의 식에서 오른쪽 식의 두번째 항의 기대값이 0 이므로

\[ E(y_i - \mu_i) =0 \]

\(y_i - \mu_i\) 를 오차항과 같이 생각하면 위의 모형을 오차항의 분산이 다른 선형모형으로 생각할 수 있다. 지수군 분포의 성질을 이용하면 작업 변량 \(z_i\)의 분산은 다음과 같다.

\[ Var(z_i) = Var(r_i) = [g_\mu(\mu_i)]^2 Var(y_i) = [g_\mu(\mu_i)]^2 [ a(\phi_i)v(\mu_i)] \tag{11.11}\]

이러한 가정하에서 작업 변량 \(z_i\)를 반응변수로 놓고 분산의 역수를 가중치로하는 가중 선형모형(wighted linear regression) 을 최소제곱법으로 적합하고 계수의 값이 수렴할 때까지 반복적으로 수행하는하는 계산법을 제공하였다. 이러한 방법을 반복가중최소최곱법(iterative weighted least square method; IWLS)라고 부른다.

(Searle 와/과 McCulloch 2001 의 136 쪽 참조)

11.3 최대가능도추정

이제 회귀계수 \(\beta\)를 최대가능도추정법(Maximum Likelihood Estimation)으로 구하기 위하여 로그가능도함수 식 11.7 를 회귀계수 벡터 \(\beta\)로 미분한 가능도함수 방정식을 고려하자.

\[ \pardifftwo{ \ell_n }{ \beta} = 0 \tag{11.12}\]

중요

여기서 일반화 선형모형에서 나타나는 모수들 \(\beta\), \(\mu_i\), \(\theta_i\)의 관계를 살펴보자.

회귀계수 벡터 \(\beta\)는 설명변수 벡터 \({ x}_i\)와 내적 형태로 연결되어 있으며 이를 선형 예측식이라고 한다.

\[ \eta_i = { x}^t_i \beta \]

선형 예측식 \(\eta_i\)는 관측값의 평균 \(\mu_i\)와 연결함수 \(g\)로 연결되어 있다.

\[ g(\mu_i) = \eta_i \]

관측값의 평균 \(\mu_i\)는 기본형 모수 \(\theta_i\)와 함수 \(b\)로 연결되어 있다.

\[ b'(\theta_i) = \mu_i \]

일반화 선형모형에서 최종적으로 추정해야 하는 모수는 회귀계수 벡터 \(\beta\) 이며 모수 \(\beta\), \(\mu_i\), \(\theta_i\) 다음과 연결되어 있음을 알 수 있다.

\[ \beta \underset{ g}{\longrightarrow} \mu_i \underset{ b }{\longrightarrow} \theta_i \tag{11.13}\]

최대가능도 추정량을 구하는 방정식을 유도할 때 다음과 같은 일반화 선형모형의 지수군 분포에서 나타나는 미분공식이 적용된다.

선형 예측식의 미분:

\[ \pardifftwo{ \eta }{ \beta} = \pardifftwo{ { x}^t \beta }{ \beta } = x \]

연결함수의 역함수에 대한 미분: \(g(\mu) = \eta\) 의 관계를 이용하면

\[ \pardifftwo{\mu}{\eta} = \pardifftwo{\mu}{g(\mu)} = \left [ \pardifftwo{g(\mu)}{\mu} \right ]^{-1} = [g'(\mu)]^{-1} = g^{-1}_{\mu} (\mu) \]

여기서 \(g_{\mu}(\mu) = g'(\mu)\)로서 연결함수의 미분을 나타내는 기호이다.

평균과 기본형모수의 미분, 분산함수: 평균과 분산의 관계식을 이용하면

\[ \frac{\partial \theta }{\partial \mu } = \left [ \frac{\partial \mu }{\partial \theta } \right ]^{-1} = \left [ b''(\theta) \right ]^{-1} = \frac{1}{v(\mu)} \]

11.3.1 가능도 방정식의 유도: 첫 번째 방법

이제 가능도 함수 식 11.7 의 형태를 이용하여 방정식 식 11.12 를 유도해 보자.

\[ \begin{aligned} 0 & =\pardifftwo{ \ell_n}{ \beta }\\ &= \sum_{i=1}^n \left [ \frac{1}{a(\psi_i)} \right ] \left [ y_i \pardifftwo{ \theta_i}{ \beta } -\pardifftwo{ b(\theta_i)}{ \beta } \right ] \\ &= \sum_{i=1}^n \left [ \frac{1}{a(\psi_i)} \right ] \left [ y_i \pardifftwo{ \theta_i}{ \beta } -\pardifftwo{ \theta_i }{ \beta } \pardifftwo{ b(\theta_i)}{ \theta_i } \right ] \\ &= \sum_{i=1}^n \left [ \frac{1}{a(\psi_i)} \right ] \left [ \pardifftwo{ \theta_i }{ \beta } (y_i - \mu_i) \right ] \\ &= \sum_{i=1}^n \left [ \frac{1}{a(\psi_i)} \right ] \left [ \pardifftwo{ \eta_i }{ \beta } \pardifftwo{ \mu_i }{ \eta_i} \pardifftwo{ \theta_i }{ \mu_i}(y_i - \mu_i) \right ] \\ &= \sum_{i=1}^n \left [ \frac{1}{a(\psi_i)} \right ] \left [ { x}_i \frac{ (y_i - \mu_i) }{ v(\mu_i) g_\mu(\mu_i) } \right ] \\ &= \sum_{i=1}^n \left [ x_i w_i g_\mu(\mu_i) (y_i - \mu_i) \right ] \\ & = \begin{bmatrix} \sum_{i=1}^n x_{i1} w_i g_{\mu}(\mu_i) (y_i - \mu_i) \\ \sum_{i=1}^n x_{i2} w_i g_{\mu}(\mu_i) (y_i - \mu_i) \\ \vdots \\ \sum_{i=1}^n x_{ip} w_i g_{\mu}(\mu_i) (y_i - \mu_i) \end{bmatrix} \end{aligned} \tag{11.14}\]

여기서 가중치 \(w_i\)는 다음과 정의한다. 가중치 \(w_i\)는 앞에서 설명한 작업 변량의 분산의 역수와 동일하다. 식 11.11 을 참조하자.

\[ w_i \equiv \frac{1}{ g^2_\mu(\mu_i) a(\phi_i) v(\mu_i) } \]

11.3.2 가능도 방정식의 유도: 두 번째 방법

위의 방정식 식 11.12 에 미분의 연쇄법칙(chain rule)을 적용하면 다음과 같은 방정식을 얻는다.

\[ \pardifftwo{ \ell_n }{ \beta} = \pardifftwo{ \eta }{ \beta } \pardifftwo{ \mu }{ \eta } \pardifftwo{ \theta }{ \mu } \pardifftwo{ \ell_n }{ \theta } =0 \tag{11.15}\]

위의 식에서 \(\eta\), \(\mu\), \(\theta\)는 다음과 같이 \(n\)개의 대응되는 원소로 이루어진 벡터이다.

\[ \eta = \begin{bmatrix} \eta_1 \\ \eta_2 \\ \vdots \\ \eta_n \end{bmatrix}, \quad \mu = \begin{bmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_n \end{bmatrix}, \quad \theta = \begin{bmatrix} \theta_1 \\ \theta_2 \\ \vdots \\ \theta_n \end{bmatrix} \]

이제 식 11.15 에서 나타난 도함수들을 각각 구해보자. 먼저

\[ \begin{aligned} \pardifftwo{ \eta }{ \beta } & = \left [ \pardifftwo{ \eta_1 }{ \beta } ~~ \pardifftwo{ \eta_2 }{ \beta } ~~ \dots ~~ \pardifftwo{ \eta_n }{ \beta } \right ] \\ & = \left [ \pardifftwo{ { x}^t_1 \beta }{ \beta } ~~ \pardifftwo{ { x}^t_2 \beta }{ \beta } ~~ \dots ~~ \pardifftwo{ { x}^t_n \beta }{ \beta } \right ] \\ & = [ { x}_1 ~~ { x}_2 ~~ \dots ~~ { x}_n ] \\ & = { X}^t \end{aligned} \]

위의 식에서 행렬 \(X\)는 \(n \times p\) 계획행렬이다. 또한

\[ \begin{aligned} \pardifftwo{ \mu }{ \eta } & = \left [ \pardifftwo{ \mu_1 }{ \eta } ~~ \pardifftwo{ \mu_2 }{ \eta } ~~ \dots ~~ \pardifftwo{ \mu_n }{ \eta } \right ] \\ & = \left [ \pardifftwo{ g^{-1}(\eta_1) }{ \eta } ~~ \pardifftwo{ g^{-1}(\eta_2) }{ \eta } ~~ \dots ~~ \pardifftwo{ g^{-1}(\eta_2) }{ \eta } \right ] \\ & = \begin{bmatrix} \frac{1}{g'(\mu_1)} & 0 & \dots & 0 \\ 0 & \frac{1}{g'(\mu_2)} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \frac{1}{g'(\mu_n)} \end{bmatrix} \\ & = \begin{bmatrix} g^{-1}_{\mu}(\mu_1) & & & \\ & g^{-1}_{\mu}(\mu_2) & & \\ & & \ddots & \\ & & & g^{-1}_{\mu}(\mu_n) \end{bmatrix} \end{aligned} \]

위에서 \(\pardifftwo{ \mu }{ \eta }\)는 \(n\)-차원 대각행렬이며 \(g_\mu(\mu) = g'(\mu)\)로서 연결함수 \(g\)를 1차 미분한 함수이다.

또한 다음과 같은 결과를 얻는다.

\[ \begin{aligned} \pardifftwo{ \theta }{ \mu } & = \left [ \pardifftwo{ \theta_1 }{ \mu } ~~ \pardifftwo{ \theta_2 }{ \mu } ~~ \dots ~~ \pardifftwo{ \theta_n }{ \mu } \right ] \\ & = \begin{bmatrix} \frac{1}{b''(\theta_1)} & 0 & \dots & 0 \\ 0 & \frac{1}{b''(\theta_2)} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \frac{1}{b''(\theta_n)} \end{bmatrix} \\ & = \begin{bmatrix} v^{-1}(\mu_1) & & & \\ & v^{-1}(\mu_2) & & \\ & & \ddots & \\ & & & v^{-1}(\mu_n) \end{bmatrix} \end{aligned} \]

마지막으로 식 11.7 를 이용하면 로그가능도함수 \(\ell_n\) 을 \(\theta\)로 미분한 \(n\)-차원 벡터는 다음과 같이 얻어진다.

\[ \pardifftwo{ \ell_n }{ \theta} = \begin{bmatrix} \frac{y_1 - b'(\theta_1)}{a(\phi_1)} \\ \frac{y_2 - b'(\theta_2)}{a(\phi_2)} \\ \vdots \\ \frac{y_n - b'(\theta_n)}{a(\phi_n)} \end{bmatrix} = \begin{bmatrix} \frac{y_1 - \mu_1}{a(\phi_1)} \\ \frac{y_2 - \mu_2}{a(\phi_2)} \\ \vdots \\ \frac{y_n - \mu_n}{a(\phi_n)} \end{bmatrix} \]

이제 가능도추정을 위한 방정식 식 11.15 을 위에서 유도한 도함수 벡터와 행렬을 이용하여 다시 쓰면 다음과 같다.

\[ \begin{aligned} 0 & = \pardifftwo{ \ell_n }{ \beta} = \pardifftwo{ \eta }{ \beta } \pardifftwo{ \mu }{ \eta } \pardifftwo{ \theta }{ \mu } \pardifftwo{ \ell }{ \theta } \\ & = { X}^t \begin{bmatrix} \frac{1}{g_{\mu}(\mu_1)} & & & \\ & \frac{1}{g_{\mu}(\mu_1)} & & \\ & & \ddots & \\ & & & \frac{1}{g_{\mu}(\mu_n)} \end{bmatrix} \begin{bmatrix} \frac{1}{v(\mu_1)} & & & \\ & \frac{1}{v(\mu_2)} & & \\ & & \ddots & \\ & & & \frac{1}{v(\mu_n)} \end{bmatrix} \begin{bmatrix} \frac{y_1 - \mu_1}{a(\phi_1)} \\ \frac{y_2 - \mu_2}{a(\phi_2)} \\ \vdots \\ \frac{y_n - \mu_n}{a(\phi_n)} \end{bmatrix} \\ &= { X}^t \begin{bmatrix} \frac{1}{g^2_{\mu}(\mu_1) a(\phi_1) v(\mu_1)} & & & \\ & \frac{1}{g^2_{\mu}(\mu_1)a(\phi_2)v(\mu_2)} & & \\ & & \ddots & \\ & & & \frac{1}{g^2_{\mu}a(\phi_n)(\mu_n)v(\mu_n)} \end{bmatrix} \\ & \quad \quad \times \begin{bmatrix} g_{\mu}(\mu_1) & & & \\ & g_{\mu}(\mu_2) & & \\ & & \ddots & \\ & & & g_{\mu}(\mu_n) \end{bmatrix} \begin{bmatrix} y_1 - \mu_1 \\ y_2 - \mu_2 \\ \vdots \\ y_n - \mu_n \end{bmatrix} \\ & = { X}^t W \Delta ( y - \mu) \end{aligned} \tag{11.16}\]

위의 식에서 가중값 대각행렬 \(W\), 연결함수 미분값 대각행렬 \(\Delta\), 관측값 벡터 \(y\), 평균 벡터 \(\mu\)는 다음과 같이 정의된다.

\[ \begin{aligned} W & = \begin{bmatrix} w_1 & & & \\ & w_2 & & \\ & & \ddots & \\ & & & w_n \end{bmatrix} = \begin{bmatrix} \frac{1}{g^2_{\mu}(\mu_1) a(\phi_1) v(\mu_1)} & & & \\ & \frac{1}{g^2_{\mu}(\mu_1)a(\phi_2)v(\mu_2)} & & \\ & & \ddots & \\ & & & \frac{1}{g^2_{\mu}a(\phi_n)(\mu_n)v(\mu_n)} \end{bmatrix} \\ \Delta & = \begin{bmatrix} \delta_1 & & & \\ & \delta_2 & & \\ & & \ddots & \\ & & & \delta_n \end{bmatrix} = \begin{bmatrix} g_{\mu}(\mu_1) & & & \\ & g_{\mu}(\mu_2) & & \\ & & \ddots & \\ & & & g_{\mu}(\mu_n) \end{bmatrix} \\ y & = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, \quad \mu = \begin{bmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_n \end{bmatrix} \end{aligned} \]

11.4 최대가능도추정량의 계산

이제 회귀계수 \(\beta\)를 최대가능도추정법(Maximum Likelihood Estimation)으로 가능도 방정식을 식 11.16 를 이용하면 다음과 같은 행렬 방정식으로 표시된다.

\[ { X}^t W \Delta y = { X}^t W \Delta \mu \tag{11.17}\]

위의 방정식은 일반적으로 회귀계수 벡터 \(\beta\)에 대하여 선형방정식이 아니므로 최소제곱법과 같이 최대가능도 추정량을 직접 구할 수 없다.

노트

정규분포 가정 하에서 선형회귀 모형에서는 식 11.17 이 최소제곱법의 방정식 \({ X}^t y = { X}^t { X} \beta\)로 유도되고 직접적으로 구할 수 있다.
많은 경우 스케일 모수 \(a(\phi_i)\)는 관측값 \(y_i\)에 따라 변하지 않고 상수인 경우가 흔하다. 즉 \(a(\phi_i) \equiv a(\phi)\). 이러한 경우 가능도 방정식 식 11.14 또는 식 11.16 에서 스케일 모수 \(a(\phi_i)\)를 1로 놓고 방정식을 푼다.

최대가능도추정량을 실제 계산하기 위하여 로그 가능도 함수의 2차 도함수(헤시안) 행렬을 구해보자. 식 11.14 에서 얻은 1차 도함수를 한번 더 미분하면 다음과 같은 결과를 얻는다.

\[ \begin{aligned} \pardiffdd{\ell_n}{ \beta}{ \beta^t} & = \sum \pardifftwo{}{ \beta} \left [ { x}_i w_i g_\mu(\mu_i) (y_i - \mu_i) \right ] \\ & = \sum \pardifftwo{}{ \beta} \left [ { x}_i c_i (y_i - \mu_i) \right ] \quad [ c_i \equiv w_i g_\mu(\mu_i) ] \\ & = \sum \left [ \pardifftwo{ c_i (y_i - \mu_i)}{ \beta} { x}^t_i \right ] \\ & = \sum \left [ \pardifftwo{ c_i }{ \beta} (y_i - \mu_i)+ \pardifftwo{ (y_i - \mu_i)}{ \beta} c_i \right ] { x}^t_i \\ & = \sum \left [ \pardifftwo{ c_i }{ \beta} (y_i - \mu_i) - \pardifftwo{ \eta_i}{ \beta} \pardifftwo{ \mu_i}{\eta_i} c_i \right ] { x}^t_i \\ &= \sum \left [ \pardifftwo{ c_i }{ \beta} (y_i - \mu_i)- { x}_i [g_\mu(\mu)]^{-1} c_i \right ] { x}^t_i \quad [ c_i [g_\mu(\mu)]^{-1} = w_i] \\ &= \sum \left [ \pardifftwo{ c_i }{\beta} (y_i - \mu_i) \right ] { x}^t_i - \sum { x}_i w_i { x}^t_i \\ &= \sum \left [ \pardifftwo{ c_i }{\beta} (y_i - \mu_i) \right ] { x}^t_i - { X}^t W { X} \end{aligned} \]

그러므로 피셔정보 \(I( \beta)\)는 다음과 같이 얻어진다.

\[ \begin{aligned} I( \beta) & = - E \left [ \pardiffdd{\ell_n}{ \beta}{ \beta^t} \right ] \\ & = E \left [ - \sum \left [ \pardifftwo{ c_i }{\beta} (y_i - \mu_i) \right ] { x}^t_i + { X}^t W { X} \right ] \\ & = 0+ { X}^t W { X} \end{aligned} \tag{11.18}\]

또는 식 11.16 에서 얻은 1차 도함수 방정식을 를 한번 더 미분하여 기대값을 취하면 식 11.18 와 동일한 결과를 얻는다.

\[ \begin{aligned} \pardiffdd{\ell_n}{ \beta}{ \beta^t} & = \pardifftwo{}{ \beta} \left [ { X}^t W \Delta ( y - \mu) \right ] \\ & = \left \{ \pardifftwo{}{ \beta} \left [ { X}^t W \Delta \right ] \right \} ( y - \mu) + \left \{ \pardifftwo{}{ \beta} \left [ ( y - \mu)^t \right ] \right \} \Delta W { X} \\ & = \left \{ \pardifftwo{}{ \beta} \left [ { X}^t W \Delta \right ] \right \} ( y - \mu) -\left [ \pardifftwo{ \mu}{ \beta} \right ] \Delta W { X} \\ & = \left \{ \pardifftwo{}{ \beta} \left [ { X}^t W \Delta \right ] \right \} ( y - \mu) -\left [ \pardifftwo{ \eta}{ \beta} \pardifftwo{ \mu}{ \eta} \right ] \Delta W { X} \\ & = \left \{ \pardifftwo{}{ \beta} \left [ { X}^t W \Delta \right ] \right \} ( y - \mu) -\left [ { X}^t { \Delta}^{-1} \right ] \Delta W { X} \\ & = \left \{ \pardifftwo{}{ \beta} \left [ { X}^t W \Delta \right ] \right \} ( y - \mu) -{ X}^t W { X} \\ \end{aligned} \]

최대가능도추정량 \(\hat { \beta}\)는 가능도 방정식 식 11.17 을 직접 풀어서 계산할 수 있지만 대부분의 경우 직접해(explicit solution)를 구하는 것이 불가능 하다. 따라서 보통의 경우 선형화된 작업변량에 반복가중최소제곱법(iterative weighted least square; IWLS)을 적용하여 최대가능도 추정량을 구하며 IWLS로 구하는 해는 가능도 방정식 식 11.17 의 해와 동일하다.

주어진 분포에서 기본 연결함수를 \(g\) 라고 하고 관측값 \(y\)를 변환한 작업변량 \(z=g(y)\)의 테일러 전개를 다음과 같이 고려해 보자.

\[ z \equiv g(y) \cong g(\mu) + g_\mu (\mu)(y-\mu) \]

작업 변량 \(z\) 의 분산은 식 11.11 와 같이 다음으로 주어진다.

\[ var(z) = v(\mu) g^2_\mu (\mu) \equiv w^{-1} \]

회귀계수 벡터의 초기값을 \({ \beta}_0\)라고 하자. 그러면 작업 변량의 초기값 \(z_0\)는 \({ \beta}_0\) 로 계산된 \(\mu_0\)를 이용하여 다음과 같이 구할 수 있다.

\[ z_0 = g(\mu_0) + g_\mu (\mu_0) (y-\mu_0) = \eta_0 + g_\mu (\mu_0) (y-\mu_0) \]

IWLS 추정량 \(\hat { \beta}\)는 \(z_0\) 를 설명변수 벡터 \(x\)로 선형회귀분석을 적합할 때 가중치를 \(w_0\) 로 이용하는 가중최소제곱법으로 반복적으로 적용하여 개선할 수 있다.

실제로 IWLS 추정량은 피셔정보를 이용한 스코어 방법(Fisher scoring method)로 구한 최대가능도 추정량과 동일함을 보일 수 있다. 일단 회귀계수 벡터의 초기값 \(\hat { \beta}^0\) 으로 계산된 피셔정보 행렬을 \(A\)로 아래와 같이 정의하자.

\[ A = -E \left [ \pardiffdd{\ell_n}{ \beta}{ { \beta}^t} \right ]_{ \beta= \hat { \beta}^0} \]

새로운 추정량 \(\hat { \beta}^1\) 가 이전의 추정량 \(\hat { \beta}^0\)에서 다음과 같은 축차식으로 계산되는 방법이 피셔 스코링 방법이다.

\[ 0 = \pardifftwo{\ell_n }{ \beta} |_{ \beta=\hat { \beta}^0}- A(\hat { \beta}^1 - \hat { \beta}^0) \quad \Leftrightarrow \quad A(\hat { \beta}^1 - \hat { \beta}^0) = \pardifftwo{\ell_n }{ \beta} |_{ \beta=\hat { \beta}^0} \]

위의 피셔의 스코링 방법을 더 정리하면 다음과 같이 유도할 수 있다.

\[ \begin{aligned} A( \hat { \beta}^1 - \hat { \beta}^0) & = \pardifftwo{ \ell_n }{ \beta} |_{ \beta=\hat { \beta}^0} \\ \Leftrightarrow { X}^t { W}_0 X (\hat { \beta}^1 - \hat { \beta}^0) & = { X}^t { W}_0 { \Delta}_0 ( y-{ \mu}_0) \\ \Leftrightarrow { X}^t { W}_0 X \hat { \beta}^1 & = { X}^t { W}_0 [ X \hat { \beta}^0+ { \Delta}_0 ( y-{ \mu}_0)] \\ \Leftrightarrow { X}^t { W}_0 X \hat { \beta}^1 & = { X}^t { W}_0 [ { \eta}_0+ { \Delta}_0 ( y-{ \mu}_0)] \\ \Leftrightarrow { X}^t { W}_0 X \hat { \beta}^1 & = { X}^t { W}_0 { z}_0 \\ \end{aligned} \]

위의 방정식은 \(z_0 = \eta_0 + g_\mu (\mu_0) (y-\mu_0)\)를 가중치 \(w_0\)를 사용하여 얻은 가중최소제곱법에서 나온 방정식임을 알 수 있다. 따라서 최대가능도 추정량을 구하는 피셔의 스코링 방법은 앞에서 알아본 반복가중최소제곱법과 동일하다.

11.5 편차

선형모형에서 잔차제곱합(residual sum of square; SSE)에 대한 의미를 살펴보고 이를 일반화 선형모형에 확장하는 개념인 편차(deviance)의 정의를 알아보자.

먼저 다음과 같은 선형회귀식을 고려한다.

\[ y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \dots \beta_p x_{pi} + e_i , \quad i=1,2,\dots,n \]

여기서 오차항 \(e_i\)를 서로 독립이며 평균이 0이고 분산이 \(\sigma^2\)인 정규분포를 따른다고 가정하고 (\(\sigma^2\)는 알고있다고 가정하자) 각 관측변수의 평균을 다음과 같이 \(\mu_i\)로 하자.

\[ \mu_i =E(y_i|x_i)= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \dots +\beta_p x_{pi} \]

서로 독립인 관측변수 \(y_i\)의 분포는 정규분포를 따르므로

\[ y_i \sim N(\mu_i,\sigma^2) \]

관측치 \(\pmb y=(y_1,y_2,\dots,y_n)^t\)의 로그가능도함수는 다음과 같이 나타낼 수 있다.

\[ l(\pmb \mu|\pmb y) = C -\frac{n}{2} \log \sigma^2 -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu_i)^2 \]

예측변수 \(x_1,x_2,\dots,x_p\)를 고려한 선형회귀모형에서 각 반응변수 평균의 예측에 선형예측식을 사용하며 모형에서 모수의 개수는 \(p+1\)로 보통 모수의 수가 관측값의 수 \(n\)보다 작다.

\[ \hat \mu_i = \hat \beta_0 + \hat \beta_1 x_{1i} + \hat \beta_2 x_{2i} + \dots +\hat \beta_p x_{pi} \equiv \hat y_i \]

이 때 선형회귀모형의 로그가능도함수의 최대값는 다음과 같다.

\[ \begin{aligned} l(\hat{\pmb \mu}|\pmb y) & = l_{regession}(\hat {\pmb \beta}|y) \\ & = C -\frac{n}{2} \log \sigma^2 -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \hat \mu_i)^2 \\ &= C -\frac{n}{2} \log \sigma^2 -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \hat y_i)^2 \\ &= C -\frac{n}{2} \log \sigma^2 -\frac{1}{2\sigma^2} SSE \end{aligned} \]

이제 위의 선형회귀모형에서 예측변수 \(x_1,x_2,\dots,x_p\)를 고려하지 않는 포화 모형을 생각해보자.

\[ y_i = \beta_{0i} +e_i \quad \text{or} \quad E(y_i) = \beta_{0i} \]

이러한 포화 모형은 \(n\)개의 반응변수의 평균을 \(n\)개의 모수를 가진 모형으로 추정하는 것으로 위와 같은 모형을 포화 모형(saturated model)이라고 한다. 포화 모형에서 모수 \(\beta_{0i}\)의 최소제곱 추정량(또는 최대가능도 추정량)는 관측값 \(y_i\)임을 쉽게 알수 있다.

\[ \min_{\beta_{0i}} \sum_{i=1}^n (y_i -\beta_{0i})^2 \quad \Rightarrow \hat \beta_{0i} = y_i, \quad i=1,2,\dots,n \]

포화 모형의 의미는 우리가 생각할 수 있는 모형 중에 가장 큰 모형으로 포화모형보다 큰 모형을 생각할 수 없다. 위에서 언급한 바와 같이 \(n\)개의 관측값에 대하여 모수의 수가 \(n\)개보다 큰 모형을 생각하면 유일한 모수의 추정이 불가능하다.

선형회귀모형에서 포화모형은 \(\hat \beta_{0i}=y_i\)이며 로그가능도함수의 최대값은 \(l(\pmb y | \pmb y)\)로 표시하며 다음과 같다.

\[ \begin{aligned} l(\pmb y | \pmb y) & = l_{saturated}(\hat {\pmb \beta_{0}} | y) \\ & = C -\frac{n}{2} \log \sigma^2 -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \hat \beta_{0i})^2 \\ &= C -\frac{n}{2} \log \sigma^2 -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - y_i)^2 \\ &= C -\frac{n}{2} \log \sigma^2 + 0 \end{aligned} \]

포화모형은 설정할 수 있는 최대의 모수를 가진 가장 큰 모형이므로 우리가 생각할 수 있는 모형 중에서 관측값을 예측하는 예측력은 가장 좋다는 것을 알 수 있다(하지만 과적합모형이다). 따라서 예측변수 \(\pmb x\)들을 사용하는 선형회귀모형의 예측력이 포화모형이 가지는 예측력에 가까우면 좋은 모형이라고 생각할 수 있다. 반응변수의 평균을 예측하는 예측력은 로그가능도함수의 크기로서 나타낼 수 있다. 포화모형과 선형회귀모형의 로그가능도함수를 비교하면 포화모형의 로그가능도함수가 크다는 것을 알수 있고 (why?) 두 로그가능도함수의 의 차이를 비교하면 다음과 같다.

\[ l(\pmb y | \pmb y) - l(\hat {\pmb \mu} | \pmb y) = l_{saturated}(\hat {\pmb \beta_{0}} | y) -l_{regession}(\hat {\pmb \beta}|y) = \frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \hat y_i)^2 = \frac{1}{2\sigma^2} SSE \]

따라서 포화모형과 로그가능도함수의의 차이가 작다는 것은 선형회귀모형의 잔차제곱합(SSE)이 작다는 것을 의미한다. 보통 잔차제곱합이 작으면 선형회귀모형의 예측력이 좋은 모형이며 이는 선형회귀모형의 가능도함수가 포화모형의 가능도함수에 가깝다는 의미이다.

이렇게 모형의 예측능력을 평가하는 측도로서 편차(deviance)를 포화모형과 고려한 회귀모형의 로그가능도함수의 차이에 2를 곱한 양으로 정의한다. 따라서 편차는 작을 수록 좋다.

\[ deviance \equiv D(\pmb y;\hat {\pmb \mu}) = 2 \left [ l(\pmb y | \pmb y) - l(\hat {\pmb \mu} | \pmb y) \right ] \tag{11.19}\]

정규분포인 경우 편차는 다음과 같이 주어진다.

\[ D(\pmb y;\hat {\pmb \mu}) = 2 \left [ l(\pmb y | \pmb y) - l(\hat {\pmb \mu} | \pmb y) \right ] = \frac{1}{\sigma^2} SSE \]

이제 이항분포들에서 나온 관측값에 대한 포화모형을 생각해 보자.

\[ y_i \sim B(m_i, p_i(x_i)), \quad i=1,2,\dots,n \]

위의 모형에서 포화모형은 어떤 모형일까? 포화모형은 \(n\)개의 관측변수의 평균, 여기서 \(E(y_i/m_i) = p(\pmb x_i)\)를 \(n\)개의 관측치 \(y_i\)를 이용하여 추정한 모형으로서 각 성공확률은 해당하는 관측된 성공의 비율에 의해 추정된다. 즉,

\[ \hat p(x_i) = \frac{y_i}{m_i} \]

이러한 경우의 로그가능도함수의 값은 다음과 같이 주어진다.

\[ \begin{aligned} l(\pmb y | \pmb y) & =l_{saturated} \\ & = \sum_i \log {{m_i}\choose{y_i}} + \sum_i y_i \log \hat p(x_i) + \sum_i (m_i -y_i) \log (1- \hat p(x_i)) \\ & = \sum_i \log {{m_i}\choose{y_i}} + \sum_i y_i \log \frac{y_i}{m_i} + \sum_i (m_i -y_i) \log (1- \frac{y_i}{m_i}) \end{aligned} \]

따라서 위에서 주어진 포화함수의 로그가능도함수에서 로지스틱회귀식의 로그가능도함수 식 11.10 를 빼고 2를 곱해서 편차를 정의할 수 있다.

\[ \begin{aligned} D(\pmb y;\hat {\pmb \mu}) &= 2 \left [ l(\pmb y | \pmb y) - l(\hat {\pmb \mu} | \pmb y) \right ] \\ &= 2 ( l_{saturated}-l_{regession}) \\ & = 2 \left [\sum_i y_i \log \frac{y_i}{m_i} + \sum_i (m_i -y_i) \log (1- \frac{y_i}{m_i}) -\sum_i y_i \log \hat p(x_i) - \sum_i (m_i -y_i) \log (1-\hat p(x_i)) \right ]\\ & = 2\left [\sum_i y_i \log \frac{y_i}{m_i \hat p(x_i) } + \sum_i (m_i -y_i) \log \frac{1-y_i/m_i}{1-\hat p(x_i)}\right ] \\ & = 2\left [\sum_i y_i \log \frac{y_i}{m_i \hat p(x_i) } + \sum_i (m_i -y_i) \log \frac{m_i-y_i}{m_i- m_i \hat p(x_i)} \right ] \\ & = 2\left [\sum_i y_i \log \frac{y_i}{\hat y_i } + \sum_i (m_i -y_i) \log \frac{m_i-y_i}{m_i- \hat y_i } \right ]\\ \end{aligned} \]

위에서 \(\hat y_i = m_i \hat p(x_i)\)으로 로지스틱 회귀에서 성공의 횟수의 평균에 대한 예측값이다.

위의 논의에서 알 수 있듯이 로지스틱 회귀에서의 편차는 선형회귀 분석에서 잔차 제곱합 SSE의 의미로 해석할 수 있으며 작을 수록 모형의 예측력이 좋다는 것을 알 수 있다.

또한 편차는 표본의 개수가 충분히 크면 자유도가 \(n-p\) 인 카이제곱분포를 따른다. 여기서 \(p\)는 회귀계수 벡터 \(\pmb \beta\)의 크기이다.

정규분포와 이항분포의 편차를 비교하면 정규분포의 편차에는 산포를 나타내는 모수 \(\sigma^2\) 이 포함되어 있지만 이항분포의 편차에는 다른 모수가 나타나지 않는다. 식 11.19 에서 주어진 편차를 척도 모수(scaled parameter) 또는 산포 모수(dispersion parameter) \(\phi\) 를 곱해준 값을 척도화 편차(scaled deviance) \(D^*(\pmb y;\hat {\pmb \mu})\) 라고 부른다.

\[ D^*(\pmb y;\hat {\pmb \mu}) = \phi D(\pmb y;\hat {\pmb \mu}) \tag{11.20}\]

정규분포에서 산포 모수가 분산 \(\phi=\sigma^2\) 이므로 척도화 편차는 잔차제곱합 \(D^*(\pmb y;\hat {\pmb \mu}) = SSE\) 가 되며 이항분포에서는 편차와 척도화 편차가 같다.

이제 선형모형에서 고려하는 다음과 같은 가설 검정을 고려해 보자.

\[ H_0 : \text{ reduced model} \quad vs. \quad H_1: \text{full model} \tag{11.21}\]

예를 들어 다음과 같은 가설검정을 자주 고려하게 된다.

\[ H_0: \beta_{1} = \beta_2 =\dots =\beta_q =0 \quad vs. H_1: \text{ not } H_0 \]

만약 축소모형(reduced moldel) 에 대한 편차를 \(D_0\) 라고 하고 큰 모형(full model)에 대한 척도화 편차를 \(D_1\) 라고 하면 귀무가설이 참인 경우 두 편차의 차이 \(D_0 - D_1\) 은 근사적으로 자유도가 \(d\) 인 카이제곱분포를 따른다. 여기서 자유도 \(d\) 는 두 모형의 회귀계수의 갯수 차이이다. 이러한 두 편차의 차이의 점근적 분포 가능도비검정 이론에 의하여 유도할 수 있다.

11.6 이항변수에 대한 회귀모형

11.6.1 단순 로지스틱 회귀모형

이제 반응변수의 값이 연속형 변수가 아니라 두 개의 가능한 결과만을 가지는 이항변수라면 선형 회귀식은 적절하지 못하다. 왜냐하면 반응변수의 기대값이 0과 1사이의 확률로 나타나기 때문이다.

\[ E(y|x) = 1\cdot P(y=1|x) + 0 \cdot P(y=0|x) = P(y=1|x) \]

따라서 반응변수의 기대값의 범위와 예측변수가 있는 선형예측식(linear predictor) \(\beta_0 + \beta_1 x\)의 범위가 일치하지 않아서 선형회귀식을 그대로 사용할 수 없다.

위의 문제를 해결하기 위한 방법중의 하나는 다음과 같은 함수 \(m\)를 생각하여 변환된 선형예측식의 범위를 \([0,1]\)로 만드는 것이다.

\[ m:\Re \rightarrow [0,1] \quad \text{and } g(x) \text{ is monotone function}. \]

따라서 다음과 같은 이항변수를 반응변수로 하는 새로운 회귀식을 만들 수있다.

\[ E(y|x) = m(\beta_0 + \beta_1 x) \tag{11.22}\]

주로 쓰이는 변환함수로 다음과 같은 로지스틱 함수(logistic function)가 있다.

\[ m(x) = \frac{ \exp(\beta_0 + \beta_1 x) }{ 1+ \exp(\beta_0 + \beta_1 x) } \tag{11.23}\]

반응변수가 베르누이분포를 따를 때 위의 로지스틱홤수를 사용하는 회귀식을 로지스틱 회귀식이라고 한다.

\[ P(y=1|x) = \frac{ \exp(\beta_0 + \beta_1 x) }{ 1+ \exp(\beta_0 + \beta_1 x) } = \{ 1+ \exp[-(\beta_0 + \beta_1 x)] \}^{-1} \tag{11.24}\]

위의 로지스틱 회귀식을 다시 역으로 정리하면 다음과 같은 식을 얻을 수 있다.

\[ \log \left [ \frac{P(y=1|x)}{1-P(y=1|x)} \right ] = \log \frac{p(x)}{1-p(x)}=\beta_0 + \beta_1 x \tag{11.25}\]

식 11.25 에서 나타난 함수 \(g(p)=\log[p/(1-p)]\)를 로짓함수(logit function)이라고 부르며 이는 로지스틱 함수의 역함수로서 0과 1 사이의 값을 가지는 확률을 실수 전체로 변환하는 함수로서 선형예측식의 범위와 일치하게 한다.

이렇게 관측값의 평균 (베르누이분포에서는 성공확률)과 선형예측식의 관계를 설정하는 함수는 식 11.4 에서 정의한 결합함수(link function)라고 한다. 아래 로짓함수는 특별히 식 식 11.5 에 정의된 것과 같이 기본형 연결함수라고 부른다.

\[ g[E(y|x)] = g(p(x)) = \log \frac{p(x)}{1-p(x)}=\beta_0 + \beta_1 x \]

다른 종류의 결합함수도 생각할 수 있다. 식 11.2 에 나타난 프로빗 함수 \(\Phi(x)=P(Z \le x)\) 또한 결합함수로 생각할 수 있다.

11.6.2 오즈비

일반적인 선형 회귀분석의 모형에서 기울기 계수 \(\beta_1\)은 기울기로서 예측변수 \(x\)의 단위가 1 증가할 때 반응변수의 평균이 \(\beta_1\)만큼 증가하는 것으로 해석할 수 있다. 하지만 로지스틱 회귀모형 식 11.24 에서는 이러한 해석을 할 수 없다.

로지스틱 회귀모형에서 기울기 \(\beta_1\)의 의미를 알아보기 위하여 다음과 같은 주어진 확률 \(p\)에 대한 몇 가지 함수를 알아야 한다.

먼저 하나의 확률 \(p\) 에 대한 오드(odd)는 다음과 같이 정의된다.

\[ \text{odd} = \frac{p}{1-p} \]

예로 성공의 확률이 \(1/3\)일 떄 오드는 \(1/2\)가 되며 이는 평균적으로 세 번의 시행할 때 한 번 성공하고 두 번 실패한다는 의미이다. 반대로 성공의 확률이 \(2/3\)일 떄 오드는 \(2=2/1\)가 되며 이는 평균적으로 세 번의 시행할 때 두 번 성공하고 한 번 실패한다는 의미이다. 성공의 확률이 \(1/2\)일 떄 오드는 \(1\)이 된다.

두 개의 확률 \(p_1\)과 \(p_2\) 에 대한 오즈비(odds ratio) 는 다음과 같이 정의된다.

\[ \text{odds ratio} = \frac{p_1/(1-p_1)}{p_2/(1-p_2)} \]

오즈비는 두 개의 성공 확률 \(p_1\)과 \(p_2\)를 비교할 때 쓰는 양이다. 두 개의 오드를 비율로서 비교하는 양이며 오즈비가 1일 경우에 두 확률은 같다.

이제 단순 로지스틱 회귀식 식 11.25 을 생각하고 예측변수 \(x\)를 0과 1의 값을 가지는 이항변수로 가정한다. \(x=1\)인 경우는

\[ \frac{P(y=1|x=1)}{1-P(y=1|x=1)} = \exp( \beta_0 + \beta_1) \]

이며 \(x=0\)인 경우는

\[ \frac{P(y=1|x=0)}{1-P(y=1|x=0)} = \exp(\beta_0) \]

위에서 주어진 두 개의 오드, 즉 \(x=1\)인 경우와 \(x=0\)인 경우의 두 오드의 비를 구하면 다음과 같다.

\[ \frac{ \frac{P(y=1|x=1)}{1-P(y=1|x=1)} } {\frac{P(y=1|x=0)}{1-P(y=1|x=0)}} = \exp (\beta_1) \]

이는 다시 쓰면

\[ \frac{P(y=1|x=1)}{1-P(y=1|x=1)} = \exp (\beta_1) \frac{P(y=1|x=0)}{1-P(y=1|x=0)} \]

위의 식에서 볼 때 예측변수 \(x\)가 1 의 값을 가질 때 반응 변수의 오드가 예측변수 \(x\)가 0일 경우의 오드의 \(\exp (\beta_1)\)배로 변하는 것을 알 수 있다. 따라서 \(\exp (\beta_1)\)는 반응변수의 오드의 증가량으로 볼 수 있다. 이는 두 성공확률의 오즈 비가 \(\exp (\beta_1)\)을 말한다. 위의 식에 로그를 취하면 다음과 같은 관계를 얻는다.

\[ \log \left [ \frac{P(y=1|x=1)}{1-P(y=1|x=1)} / \frac{P(y=1|x=0)}{1-P(y=1|x=0)} \right ] = \beta_1 \]

즉 오즈 비의 로그값이 단순 로지스틱 회귀식에서 기울기 \(\beta_1\)으로 나타난다.

간단한 예제를 통하여 오즈비와 로지스틱 회귀의 기울기의 관계를 명확히 해보자. 100명의 사람들을 55세 이상의 사람(\(x=1\))과 55세 미만의 사람(\(x=0\))의 그룹으로 나누었을 떄 각 그룹에서 만성심장질환(CHD)가 있는 사람(\(y=1\))과 없는 사람(\(y=0\))의 수가 다음 표에 주어져있다.

나이와 만성심장질환의 관계
	나이 \(\ge 55\) (\(x=1\))	나이 \(< 55\) (\(x=0\))	합
CHD 있음 (\(y=1\))	21	22	43
CHD 없음 (\(y=0\))	6	51	57
합	27	73	100

여기서 나이에 대한 CHD유무의 오즈비는 다음과 같이 계산된다.

\[ \text{Odds Ratio } = \frac{ \tfrac{21/27}{6/27}}{ \tfrac{22/73}{51/73}} = \frac{ \tfrac{21}{6}}{ \tfrac{22}{51}} = 8.11 \]

위의 표 \(\ref{twotwotable}\)의 자료를 이용하여 로지스틱회귀를 적합시키면 결과가 아래와 같고 회귀계수 \(\beta_1\)의 추정값은 오즈비의 로그값임을 알 수 있다.

\[ \hat \beta_1 = \log (8.11) = 2.094 \]

yes <- c(21,22) 
no <- c(6,51) 
x <- c (1,0) 
m1 <- glm(cbind(yes,no) ~ x,family=binomial()) 
summary(m1)


Call:
glm(formula = cbind(yes, no) ~ x, family = binomial())

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.8408     0.2551  -3.296  0.00098 ***
x             2.0935     0.5285   3.961 7.46e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.8704e+01  on 1  degrees of freedom
Residual deviance: 1.4211e-14  on 0  degrees of freedom
AIC: 11.987

Number of Fisher Scoring iterations: 3

11.6.3 로지스틱 회귀식 추정 및 예측

주어진 예측변수 \({\pmb x}_i=(x_{1i},\dots,x_{pi})^t\)에서 실행횟수가 \(m_i\)인 이항분포\(B(m_i, p(x_i))\)를 생각하자. \(m_i\)의 시행 중에 성공의 횟수가 \(y_i\)라고 하면 \(y_i\)의 평균과 분산은 다음과 같다.

\[ E(y_i | x_i ) = m_i p(x_i), \quad \quad Var(y_i | x_i) = m_i p(x_i) (1-p(x_i)) , \quad i=1,2,\dots,n \]

이항분포를 위한 로지스틱 회귀방정식은 선형예측식과 성공의 확률의 관계를 다음과 같이 정한다.

\[ \log \left [ \frac{p(x_i)}{1-p(x_i)} \right ] = {\pmb x}^t_i {\pmb \beta} = \beta_0 + \beta_1 x_{1i} + \dots \beta_p x_{pi} \tag{11.26}\]

서로 독립인 관측값 \((y_1,y_2,\dots,y_n)\)의 가능도함수 \(L\)은 이항분포들의 결합확률밀도함수와 같고 아래와 같이 주어지며

\[ L = \prod_{i=1}^n f(y_i|p(x_i)) = \prod_{i=1}^n \left [ {{m_i}\choose{y_i}} \right] {p(x_i)}^y_i {(1-p(x_i))}^{m-y_i} \]

로그가능도함수 \(l\) 은 다음과 같이 나타낼 수 있다.

로지스틱 회귀식이 추정된 후에 새로운 예측변수 \(x=x^*\)에 대하여 성공의 확률을 예측하고 싶다면 다음과 같은 식을 써서 예측할 수 있다.

\[ P(y=1|x=x^*) = \frac{ \exp(\hat \beta_0 + \hat \beta_1 x_1^* + \dots +\hat \beta_p x_p^* ) } { 1+ \exp(\hat \beta_0 + \hat \beta_1 x_1^* + \dots +\hat \beta_p x_p^*) } = \{ 1+ \exp[-(\hat \beta_0 + \hat \beta_1 x_1^* + \dots+ \hat \beta_p x_p^*)] \}^{-1} \]

11.7 발생 횟수에 대한 모형

11.7.1 포아송 회귀모형

반응변수 \(y\)가 어떤 사건이 일어난 횟수(count)라면 주로 포아송분포를 확률 모형으로 사용한다.

\[ P(Y=y) = f(y|\mu)= \frac{ e^{-\mu} \mu^y }{y!}, \quad y=0,1,2,\dots \tag{11.28}\]

이러한 포아송 분포에서 나온 반응변수에 대하여 예측변수 \(x\)의 영향에 대한 회귀분석을 포아송 회귀식이라고 한다. 포아송 분포의 평균 \(\mu\)는 양의 실수이고 선형예측식 \(\eta= {\pmb x}^t \pmb beta\)의 범위는 실수이기 때문에 로그함수를 결합함수(link function)으로 이용하여 회귀식을 세운다.

\[ \log E(y|{\pmb x}_i) =\log \mu({\pmb x}_i) = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \tag{11.29}\]

사실 포아송분포의 로그 가능도함수에서 로그함수가 기본형 결합함수임을 쉽게 알 수 있다. 즉, \(\pmb y=(y_1,y_2,\dots,y_n)^t\)를 서로 독립이고 평균이 \(\mu_i = \mu({\pmb x}_i)\)인 포아송 확률변수라고 한다면 로그가능도함수는 다음과 같다.

\[ \begin{aligned} l &= \log \prod_{i=1}^n f(y_i|\mu_i) \notag \\ &= \sum_{i=1}^n [y_i \log \mu_i - \mu_i - \log y_i!] \end{aligned} \]

위에서 볼 수 있듯이 충분통계량 \(y_i\)에 대응하는 모수에 대한 항은 \(\log \mu_i\)로서 이는 로그함수가 기본형 결합함수임을 나타낸다.

포아송 회귀분석는 다음과 같은 몇 가지 특징을 가지고 있다.

만약에 어떤 사건이 일어난 횟수가 몇 가지 가능한 수중 하나라면 (예: \(y \le M\)) 포아송분포를 이항분포의 근사(approximation)로 생각할 수 있다. 즉 \(n\)이 크고 성공확률 \(p\)가 작으면 이항분포는 평균이 \(\mu=np\)인 포아송 분포와 매우 가깝기 때문에 가능한 횟수가 제한되었다 하더라도 포아송 회귀식을 이용할 수 있다.
사건의 일어난 횟수가 주어진 시간의 길이에 비례하고 다른 사건과 독립이면 포아송 분포를 따른다. 또한 포아송 분포는 두 개의 사건이 일어날 때 시간 간격이 지수분포(exponential distribution)을 따른다면 주어진 시간 간격동안 일어난 사건의 횟수는 포아송 분포를 따른다.
\(y_i\)가 서로 독립이고 평균이 \(\mu_i\)인 포아송 분포를 따른다면 합 \(\sum_i y_i\)는 평균이 \(\sum_i \mu_i\)인 포아송분포를 따른다

포아송 회귀분석에서 식 11.19 에서 정의된 편차 \(D(\pmb y; \hat{\pmb \mu})\)를 구해보기 위하여 포화모형을 생각해보자. 각 관측값의 평균 \(\mu_i\)를 자신의 관측값 \(y_i\)로 추정하는 것이 포화모형이다. 따라서 포화모형의 로그가능도함수는 다음과 같이 주어지고

\[ l(\pmb y| \pmb y) = \sum_{i=1}^n [y_i \log y_i - y_i - \log y_i!] \]

포아송 회귀분석의 로그가능도함수에 적용하면 \(D(\pmb y; \hat{\pmb \mu})\)를 얻을 수 있다.

\[ \begin{aligned} D(\pmb y; \hat{\pmb \mu}) & = 2 \left [ l(\pmb y | \pmb y) - l(\hat {\pmb \mu} | \pmb y) \right ] \\ & = 2 \left \{ \sum_{i=1}^n [y_i \log y_i - y_i - \log y_i!] - \sum_{i=1}^n [y_i \log \hat \mu_i - \hat mu_i - \log y_i!] \right \} \\ &= 2 \sum_{i=1}^n [ y_i \log (y_i/\hat \mu_i) - (y_i - \hat \mu_i) ] \end{aligned} \]

보기 11.2 (포아송 회귀분석) Galapagos 군도에 있는 30개의 섬에서 사는 거북이의 개체 수를 반응변수 \(y\)로하고 5개의 지리적 변수를 예측변수로 하는 포아송 회귀식을 적합하려고 한다 (자료 출처는 Faraway (2016)).

library(faraway)
data(gala)
gala <- gala[,-2]
modp <- glm(Species ~ .,family=poisson, gala)
summary(modp)


Call:
glm(formula = Species ~ ., family = poisson, data = gala)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.155e+00  5.175e-02  60.963  < 2e-16 ***
Area        -5.799e-04  2.627e-05 -22.074  < 2e-16 ***
Elevation    3.541e-03  8.741e-05  40.507  < 2e-16 ***
Nearest      8.826e-03  1.821e-03   4.846 1.26e-06 ***
Scruz       -5.709e-03  6.256e-04  -9.126  < 2e-16 ***
Adjacent    -6.630e-04  2.933e-05 -22.608  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3510.73  on 29  degrees of freedom
Residual deviance:  716.85  on 24  degrees of freedom
AIC: 889.68

Number of Fisher Scoring iterations: 5

11.7.2 비율 모형

어떤 사건이 일어날 횟수(count)는 집단이나 시간의 크기(group size)에 의존할 수 있다. 예를 들어 각 도시의 1년 범죄 발생 횟수는 그 도시의 인구수나 크기에 비례하게 된다. 이러한 모형은 이항분포를 이용하여 분석할 수 도 있지만 사건의 발생 확률이 매우 작고 집단의 크기가 크면 포아송 근사를 통한 분석도 가능하다. 또한 어떤 경우에는 집단의 크기에 대한 정보가 부족할 수 있다.

이러한 비율에 대한 모형(Rate Models)을 나타내면 아래와 같고

\[ \log \frac {\text{ count } } { \text{ group size } } = {\pmb x}^t {\pmb \beta} \]

이는 다시 발생횟수에 대한 포아송 회귀모형의 형태로 나타내면 다음과 같이 쓸 수 있다. \[ \log \text{ count } = (1)(\log \text{ group size }) + {\pmb x}^t {\pmb \beta} \]

따라서 발생횟수에 대한 포아송 회귀분석을 적합할 때 집단의 크기를 안다면 그 로그 변환값을 회귀식에 포함하여 적합할 수 있다. 위의 식에서 알 수 있듯이 크기의 로그 변환변수는 회귀계수를 강제로 1 로 놓는 제약을 둘 수 있다. 이러한 변수를 오프셋(offset)이라고 한다.

보기 11.3 (비율 모형) 세포(cells)에 감마 방사능을 쏘였을 떄 비정상성(ca)를 나타내는 횟수에 대하여 비율 모형을 적합시켰다. 독립변수는 방사능의 양(doseamt)와 비율(doserate)이다. 여기서 세포의 수(cells)를 오프셋(offset) 변수로 사용한다 (자료 출처는 Faraway (2016)).

data(dicentric)
round(xtabs(ca/cells ~ doseamt+doserate, dicentric),2)

       doserate
doseamt  0.1 0.25  0.5    1  1.5    2  2.5    3    4
    1   0.05 0.05 0.07 0.07 0.06 0.07 0.07 0.07 0.07
    2.5 0.16 0.28 0.29 0.32 0.38 0.41 0.41 0.37 0.44
    5   0.48 0.82 0.90 0.88 1.23 1.32 1.34 1.24 1.43

dicentric$dosef <- factor(dicentric$doseamt)
rmod <- glm(ca ~ offset(log(cells))+log(doserate)*dosef, family=poisson,dicentric)
summary(rmod)


Call:
glm(formula = ca ~ offset(log(cells)) + log(doserate) * dosef, 
    family = poisson, data = dicentric)

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)            -2.74671    0.03426 -80.165  < 2e-16 ***
log(doserate)           0.07178    0.03518   2.041 0.041299 *  
dosef2.5                1.62542    0.04946  32.863  < 2e-16 ***
dosef5                  2.76109    0.04349  63.491  < 2e-16 ***
log(doserate):dosef2.5  0.16122    0.04830   3.338 0.000844 ***
log(doserate):dosef5    0.19350    0.04243   4.561  5.1e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 4753.00  on 26  degrees of freedom
Residual deviance:   21.75  on 21  degrees of freedom
AIC: 209.16

Number of Fisher Scoring iterations: 4

11.7.3 음이항 분포

베르누이 독립시행에서 \(k\)번째의 성공까지의 시행회수 \(z\)는 음이항 분포(negative bionomial)을 따른다. 음이항분포는 포아송 분포에서 모수가 감마를 따를 때 근사분포로 사용될 수 있다.

\[ P(z) = {{z-1}\choose {k-1}} p^k (1-p)^{z-k},\quad z=k,k+1,\dots \tag{11.30}\]

위의 분포에서 확률 변수와 모수를 다시 아래와 같이 정의하면

\[ y=z-k, \quad p= \frac{1}{1+\alpha} \]

\(y\)의 확률분포는 다음과 같고

\[ P(y) = {{y+k-1}\choose {k-1}} \frac{\alpha^y}{(1+\alpha)^{y+k}},\quad y=0,1,2,\dots \]

따라서 \(y\)의 평균과 분산은 다음과 같이 주어진다.

\[ E(y) = \mu =k\alpha, \quad Var(y) = k\alpha + k\alpha^2= \mu + \mu^2/k \]

또한 로그가능도함수는 다음과 같이 주어지고

\[ l= \sum_{i=1}^n \left ( y_i \log \frac{\alpha}{1+\alpha} -k \log (1+\alpha) + \sum_{j=0}^{y_i-1} \log (j+k) -\log y_i! \right ) \]

결합함수는 다음과 같다.

\[ \log \frac{\alpha}{1+\alpha} = \log \frac{\mu}{\mu+k} = \eta={\pmb x}^t {\pmb \beta} \]

보통의 경우 \(k\)는 고정된 상수로 생각할 수도 있고 또는 모수로 보고 추정할 수 도 있다.

11.7.4 영과잉모형

어떤 사건의 발생 횟수에 대한 자료를 수집할 때 영(0, zero)이 비정상적으로 많이 나타나는 경우가 있다. 만약 발생횟수의 분포를 포아송분포 식 11.28 으로 가정하면 0이 관측될 확률은 크지 않다.

\[ P(y=0) = e^{-\mu} \]

자료에서 0의 발생 빈도가 비정상적으로 많은 자료를 영과잉 자료(zero inflated data)라고 하며 이러한 자료에 포아송 분포를 그대로 적용하면 회귀 계수의 추정량에 편이(bias)가 발생할 수 있는 여러 가지 문제가 생긴다.

발생 횟수에 0이 많은 이유는 매우 다양하다. 0이 많이 발생하는 대표적인 이유를 살펴보자.

외부 요인에 의하여 사건의 발생이 제약을 받는 경우
발생은 했는데 관측이 안된 경우
원래 0이 많은 경우

이렇게 영과잉 자료를 분석할 수 있는 대표적인 모형은 영과잉 포아송 모형(zero inflated poission model; ZIP)이다. 확률변수 \(y_i\)를 사건의 발생 회수라고 하면 ZIP 모형에서 0이 관측될 확률을 다음과 같이 나타낼 수 있다.

\[ P(y=0) = P(\text{ Extra zeros }) + [1-P(\text{ Extra zeros })] P(\text{ count process gives a zero }) \]

즉 0이 관측될 확률은 별도로 나타난 0이 관찰 될 확률과 원래 확률 과정에서 0이 관찰 될 확률의 조합(mixture)으로 나타난다. 이제 \(i\)번째 관측에서 별도로 나타난 0이 관찰될 확률을 \(\pi_i\)라 하면

\[ P(y_i=0) = \pi_i + (1-\pi_i) P(\text{ count process gives a zero }) \]

더 나아가 확률 과정이 평균이 \(\mu_i\)인 포아송 분포를 따른다고 가정하면

\[ \begin{aligned} P(y_i| y_i=0) & = \pi_i + (1-\pi_i) e^{-\mu_i} \\ P(y_i| y_i>0) & = (1-\pi_i)\frac{ e^{-\mu_i} \mu_i^{y_i} }{y_i!} \end{aligned} \]

위의 분포에서 \(y\)의 평균과 분산을 구해보면 다음과 같이 주어진다.

\[ \begin{aligned} E(y_i) & = (1-\pi_i)\mu_i \\ Var(y_i) & = (1-\pi_i)\mu_i + (1-\pi_i)\pi_i \mu_i^2 \end{aligned} \]

위의 식에서 볼 수 있듯이 영과잉 포아송 모형은 과포화(overdispersion)을 보인다 (\(Var(y_i) > E(y_i)\)). 영과잉 포아송 모형에 대한 회귀분석은 별도로 나타난 0이 관측될 확률 \(\pi_i\)에 대한 로지스틱회귀와 발생회수에 대한 포아송회귀의 결합으로 분석할 수 있다.

\[ \begin{aligned} \log \frac{\pi_i} {1-\pi_i} & = {\pmb x}_b^t {\pmb \beta}_b \\ \log \mu_i & = {\pmb x}_p^t {\pmb \beta}_p \end{aligned} \]