이번 포스팅에서는 일반화선형모형(Generalized Linear Model : GLM)에 대한 우도방정식(Likelihood Equation)과 반응(Response) 변수의 확률밀도함수가 Exponential Dispersion Family일때 일반화 선형 모형을 적합하는 과정에 대해서 정리하겠다. Exponential Dispersion Family에 대한 설명은 여기를 참고하기 바란다. 이번 포스팅을 이해하기 위해서 이전 포스팅을 복습하고 오면 좋다.
[일반화 선형 모형(Generalized Linear Model)] 2. Exponential Dispersion Family에 대하여
[일반화 선형 모형(Generalized Linear Model)] 1. 일반화 선형 모형 소개
목차
1. 우도 방정식(Likelihood Equation)
1. 우도 방정식(Likelihood Equation)
지금 우리에게 데이터 $(\tilde{x}_i, y_i), i=1,2,\cdots,n$가 있다. 여기서 $\tilde{x}_i = (x_{1i}, x_{2i}, \cdots, x_{pi})^t$이다. 반응 변수 $y_i$의 확률밀도함수가 Exponential Dispersion Family이므로 $y_i$의 확률밀도함수는 다음과 같다.
$$f(y_i ; \theta_i, \phi) = \exp \{[y_i\theta_i-b(\theta_i)]/a(\phi)+c(y_i,\phi)\}\tag{1}$$
여기서 $\theta_i$는 자연모수(natural parameter), $\phi$는 분산모수(dispersion parameter)이다.
그리고 연결함수를 $g$, $E(y_i)=\mu_i$, $\eta_i = \sum_{j=1}^p\beta_jx_{ij} = g(\mu_i), i=1, \cdots, n$ 라고 하자. 연결함수 중에서 $\eta_i = \theta_i$를 만족하는 연결함수를 canonical 연결함수라 한다. 즉, $ \theta_i= g(\mu_i)$ 이면 $g$를 canonical 연결함수라고 한다.
우도방정식(Likelihood Equation)이란 로그 우도함수의 1차 편미분값을 0으로 놓았을 때 만들어지는 방정식을 말하며 우도방정식을 만족하는 파라미터가 최대우도추정량(Maximum Likelihood Estimator)이다.
먼저 $L$을 Exponential Dispersion Family의 로그 우도함수라고 한다면 $L=\sum_{i=1}^nL_i, L_i=\log f(y_i ; \theta_i, \phi)$이고 식 (1)을 이용하면 다음과 같다.
$$L=\sum_{i=1}^nL_i=\sum_{i=1}^n\frac{y_i\theta_i-b(\theta_i)}{a(\phi)}+\sum_{i=1}^nc(y_i,\phi)$$
일단 목표는 우도방정식을 만족하는 $\beta_j\;(j=1,\cdots,p)$를 찾는 것이다. 따라서 우도방정식은 다음과 같다.
$$\frac{\partial L}{\partial \beta_j}=\sum_{i=1}^n\frac{\partial L_i}{\partial \beta_j}=0, \:\: (j=1,\cdots,p)\tag{2}$$
우도방정식을 구하기 위해서 각 $i$에 대하여 $\frac{\partial L_i}{\partial \beta_j}$을 계산해야 한다. 하지만 $\beta_j$가 $L_i$에 직접적으로 의존하는 것이 아니기 때문에 바로 미분하기가 어렵다. 따라서 우리가 알고 있는 함수 관계와 Chain Rule을 이용하여 식 (2)를 정리할 것이다. 먼저 Chain Rule을 이용하면 다음과 같음을 알 수 있다.
$$\frac{\partial L_i}{\partial \beta_j}=\frac{\partial L_i}{\partial \theta_i}\frac{\partial \theta_i}{\partial \mu_i}\frac{\partial \mu_i}{\partial \eta_i}\frac{\partial \eta_i}{\partial \beta_j}$$
이제 위 식에서 각 미분요소를 구해보자. 이를 위해서는 다음의 사실을 알아야한다.
$$E\left(\frac{\partial L_i}{\partial \theta_i}\right)=0, \:\: E\left(-\frac{\partial^2 L_i}{\partial \theta_i^2}\right)=E\left(\frac{\partial L_i}{\partial \theta_i}\right)^2\tag{3}$$
$$\frac{\partial L_i}{\partial \theta_i} = \frac{y_i-b'(\theta_i)}{a(\phi)}, \:\: \frac{\partial^2 L_i}{\partial \theta_i^2}=-\frac{-b''(\theta_i)}{a(\phi)}\tag{4}$$
(3)의 첫번째 식과 (4)의 첫번째 식을 이용하면 $\mu_i=E(y_i)=b'(\theta_i)$임을 알 수 있다. 따라서 $\partial L_i/\partial \theta_i = (y_i-\mu_i)/a(\phi)$, $\partial \mu_i/\partial \theta_i = b''(\theta_i)$이 된다. 다음으로 (3)의 두번째 식과 (4)의 두번째 식을 이용하면 $b''(\theta_i)=var(y_i)/a(\phi)$임을 알 수 있다. $\eta_i = \sum_{j=1}^p\beta_jx_{ij}$ 이므로 $\partial \eta_i/\partial \beta_j=x_{ij}$이다. 이제 앞서 구한 미분요소를 이용하면 식 (3)은 다음과 같다는 것을 알 수 있다.
$$\begin{split}
\frac{\partial L_i}{\partial \beta_j}&=\frac{\partial L_i}{\partial \theta_i}\frac{\partial \theta_i}{\partial \mu_i}\frac{\partial \mu_i}{\partial \eta_i}\frac{\partial \eta_i}{\partial \beta_j} \\
& = \frac{(y_i-\mu_i)}{a(\phi)}\frac{a(\phi)}{var(y_i)}\frac{\partial \mu_i}{\partial \eta_i}x_{ij} \\
& = \frac{(y_i-\mu_i)x_{ij}}{var(y_i)}\frac{\partial \mu_i}{\partial \eta_i}
\end{split}$$
따라서 우도방정식은 다음과 같다.
$$\frac{\partial L}{\partial \beta_j}=\sum_{i=1}^n\frac{(y_i-\mu_i)x_{ij}}{var(y_i)}\frac{\partial \mu_i}{\partial \eta_i}=0, j=1, 2, \cdots, p\tag{5}$$
식 (5)를 행렬을 이용해서 표현하면 다음과 같다.
$$\begin{align}X^tDV^{-1}(\textbf{y}-\pmb{\mu}) &= X^tWD^{-1}(\textbf{y}-\pmb{\mu})\\ &=0\tag{6}\end{align}$$
여기서 $X=(x_{ij})$, $D$는 $i$번째 대각원소가 $\partial \mu_i/\partial \eta_i$인 대각행렬, $V$는 $i$번째 대각원소가 $var(y_i)$인 대각행렬, $W$는 $i$번째 대각원소가 $w_i=\frac{(\partial \mu_i/\partial\eta_i)^2}{var(y_i)}$인 대각행렬, $\textbf{y}=(y_1, y_2, \cdots y_n)^t$, $\pmb{\mu}=(\mu_1, \mu_2, \cdots , \mu_n)^t$이다.
이제 우리는 우도방정식을 만족하는 $\beta_j$를 찾으면 된다. 우도방정식에서 $\beta_j$는 $\pmb{\mu}$속에 포함되어 있음을 주목하자(자연스럽게 $D$에도 포함된다).
우도방정식이 왜 중요한지 생각해보았다. 아마도 우리가 풀려는 $\beta_j$에 대한 방정식 형태를 한눈에 알아보기 위해 필요한 것이 아닐까 생각한다. 만약 저 방정식이 $\beta=(\beta_1, \cdots, \beta_p)$에 대하여 Closed Form 형태로 정리할 수 있다면(예: 정규방정식) 한방에 최대우도추정량을 구할수 있을테고 정리가 잘안된다면 Iterative Method와 같은 수치적(Numerical)인 방법으로 최대우도추정량의 근사값을 구해야 한다. 즉, 우도방정식은 최대우도추정량을 어떻게 구해야할지 알려주는 고마운 친구인 것이다.
참고로 우리가 잘알고 있는 등분산 선형 회귀모형, 즉 $var(y_i)=\sigma^2$이고 $\mu = X\beta$(이 때 연결함수는 Identity 함수가 되므로 $\partial \mu_i/\partial \eta_i=1$이다)인 경우가 된다. 그러면 식 (6)은 우리가 잘알고 있는 정규방정식
$$X^ty=X^tX\beta$$가 된다.
2. 일반화 선형 모형 적합하기
이제 주어진 데이터를 이용하여 모형을 어떻게 적합하는지 알아보자. 여기서 모형을 적합한다는 것은 $\beta$의 최대우도추정량을 계산한다는 것이다. 즉, 앞서 살펴본 우도방정식을 만족하는 $\beta$를 구한다는 뜻이다. 일반적으로 우도방정식은 $\beta$에 대하여 직접적으로 풀 수 없는 경우가 많아 수치적(Numerical)으로 구해야하는데 Newton-Raphson 방법과 Fisher-Scoring 방법을 많이 쓴다.
1. Newton-Raphson 방법
먼저 $$\begin{equation}
\textbf{u} = \left(\frac{\partial L}{\partial \beta_1}, \cdots, \frac{\partial L}{\beta_p}\right)^t=X^tWD^{-1}(\textbf{y}-\pmb{\mu}) \\
H = (h_{ab}) = (\frac{\partial^2L}{\partial \beta_a\partial \beta_b})
\end{equation}$$라고 하자. 여기서 $H$는 Hessian 행렬이라고 한다. 이때 $\textbf{u}$의 두 번째 등식은 $y_i$의 분포가 exponential dispersion family인 경우의 성립한다.
$k$번째 스텝에서 계산된 최대우도추정량의 근사값을 $\beta^{(k)}$라 하고 $\beta^{(k)}$에서 계산된 $\textbf{u}$와 $H$를 $\textbf{u}^{(k)}$, $H^{(k)}$라 하자($\textbf{u}$와 $H$은 $\beta$에 의존한다는 것을 기억하자).
이제 Taylor 정리를 이용하면 로그우도함수의 $\beta^{(k)}$에서 이차근사식이 다음과 같다는 것을 알 수 있다.
$$L(\beta) \approx L(\beta^{(k)}) + \textbf{u}^{(k)t}(\beta-\beta^{(k)})+\frac{1}{2}(\beta-\beta^{(k)})^tH^{(k)}(\beta-\beta^{(k)})\tag{7}$$
여기서 $L(\beta)$는 $\beta$가 주어졌을 때로그 우도함수값이다.
식 (7)에서 양쪽을 다음과 같이 $\beta$에 대하여 미분한 뒤 0으로 놓고 풀어보자.
$$0 = \frac{\partial L}{\partial \beta} \approx \textbf{u}^{(k)}+H^{(k)}(\beta-\beta^{(k)})\tag{8}$$
식 (8)을 이용하면 $k+1$번째 최대우도추정량의 근사값 $\beta^{(k+1)}$을 구할 수 있다.
$$\beta^{(k+1)} = \beta^{(k)} - (H^{(k)})^{-1}\textbf{u}^{(k)}$$
대부분의 Iterative Method가 그러하듯이 Newton-Raphon 방법 또한 local maxima 문제가 존재하므로 초기값을 잘 설정해야 한다. 하지만 $X$가 full rank인 경우 GLM 모형에서 다루는 우도함수들은 Hessian 행렬이 음정치(Negative Definite)이다. 따라서 로그 우도함수는 Strictly concave가 된다. 이 경우 특정 조건하에서(이 조건은 로지스틱 회귀, Poisson Log linear 모형 등 웬만한 경우에는 다 성립한다) 최대우도추정량은 유일하게 존재한다는 사실이 알려져있으며 $k\rightarrow\infty$일 때 $\beta^{(k)}$는 최대우도추정량으로 수렴하게 된다(다시말하면 local maxima 문제는 발생하지 않는다).
2. Fisher-Scoring 방법
Fisher-Scoring 방법은 Newton-Raphson 방법에서 $H$대신 $H$의 기대값의 음수 $-E(H)$를 사용한다. $J=-E(H)$라고 했을 때 $J$는 다음과 같이 구할 수 있다.
$$\begin{split}
E\left(-\frac{\partial^2 L_i}{\partial \beta_h \partial\beta_j}\right) &= E\left[\left(\frac{\partial L_i}{\partial \beta_h}\right)\left(\frac{\partial L_i}{\partial \beta_j}\right)\right]\\
&= E\left[\frac{(y_i-\mu_i)x_{ih}}{var(y_i)}\frac{\partial \mu_i}{\partial \eta_i}\frac{(y_i-\mu_i)x_{ij}}{var(y_i)}\frac{\partial \mu_i}{\partial \eta_i}\right] \\ &= \frac{x_{ih}x_{ij}}{var(y_i)}\left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 \:\: i=1, 2, \cdots, n
\end{split}$$
$L=\sum_{i=1}^nL_i$ 이므로
$$E\left(-\frac{\partial^2 L}{\partial \beta_h\partial \beta_j}\right)=\sum_{i=1}^n\frac{x_{ih}x_{ij}}{var(y_i)}\left(\frac{\partial\mu_i}{\partial\eta_i}\right)^2$$
$w_i=\frac{(\partial \mu_i/\partial\eta_i)^2}{var(y_i)}$라고 할 때 $i$번째 대각원소가 $w_i$인 대각행렬을 $W$라고 하자. 그러면 $J$는 다음과 같다.
$$J=E(-H)=X^tWX$$
따라서 $\beta^{(k)}$가 주어졌을때 $W$를 $W^{(k)}$, $J^{(k)}=X^tW^{(k)}X$라고 한다면 $\beta^{(k+1)}$은 다음과 같이 계산한다.
$$\beta^{(k+1)} = \beta^{(k)} + (J^{(k)})^{-1}\textbf{u}^{(k)}$$ 또는
$$J^{(k)}\beta^{(k+1)}=J^{(k)}\beta^{(k)}+\textbf{u}^{(k)}\tag{9}$$
이것이 Fisher-Scoring을 이용하여 최대우도추정량을 구하는 알고리즘이다. 한편 (9)의 오른쪽 부분은 다음과 같이 쓸 수 있다.
$$\begin{align}J^{(k)}\beta^{(k)}+\textbf{u}^{(k)}&=(X^tW^{(k)}X)\beta^{(k)}+X^tW^{(k)}(D^{(k)})^{-1}(\textbf{y}-\pmb{\mu}) \\ &= X^tW^{(k)}\left[X\beta^{(k)}+(D^{(k)})^{-1}(\textbf{y}-\pmb{\mu})\right]=X^tW^{(k)}z^{(k)}\end{align}\tag{10}$$
(9)의 왼쪽 부분과 $J^{(k)}=X^tW^{(k)}X$, 그리고 (10)을 이용하면 $\beta^{(k+1)}$은 다음과 같이 구할 수 있다.
$$\beta^{(k+1)} = (X^tW^{(k)}X)^{-1}X^tW^{(k)}z^{(k)}\tag{11}$$
식 (11)이 의미하는 것은 $z^{(k)}$를 반응 변수, 설명 변수 행렬이 $X$라고 했을 때 가중 최소 제곱 추정량이 $\beta^{(k+1)}$라는 것이다.
3. 초기값 및 알고리즘 비교
- 초기값 설정
알고리즘을 알았다면 초기값은 어떻게 주어야할까? $\mu$에 대한 초기값으로 반응변수 벡터 $\textbf{y}$를 고려한다. 그리고 연결함수 $g$와 설명변수로 이루어진 행렬 $X$이 full-rank임을 이용하면 $\beta$에 대한 초기값이 결정된다.
$$\beta^{(0)}=(X^tX)^{-1}X^tg(\textbf{y})$$
여기서 $g(\textbf{y})=(g(y_1), \ldots g(y_n))^t$
(초기값 설정이 어렵다면 $\beta^{(0)}$을 0으로 세팅해놓고 돌려도 잘 작동하는 경우가 많다.)
- 알고리즘 비교
Newton-Raphson 방법의 경우 Fisher-Scoring 방법보다 빠른 수렴속도를 갖고 있지만 초기값을 잘못 설정할 경우 수렴이 잘안된다고 알려져있다. Fisher-Scoring 방법은 초기값에 대하여 Robust한 방법이다.그렇다면 어떤 알고리즘을 써야할지에 대한 질문이 나올 수있다. 대부분의 GLM 모형에서 연결함수를 Canonical 연결함수를 사용하는데 이경우 Newton-Raphson 방법과 Fisher-Scoring 방법은 정확히 일치하므로 어떤 것을 사용해야할지에 대한 문제는 사라진다. 그 외의 경우에서는 계산이 쉽고 초기값에 민감하지 않은 Fisher-Scoring 방법을 사용할 것을 추천한다.
참고자료
Keith Knight - Mathematical Statistics
Alan Agresti - Foundations of Linear and Generalized Linear Models
댓글