반응 변수 중에는 특정 기간 동안에 발생한 특정 사건 횟수 정보가 담긴 경우가 종종 있다. 예를 들어 1주일에 음주 횟수, 담배 흡연 횟수 등이 있다. 이처럼 횟수 정보를 가지는 반응 변수는 포아송(Poisson) 분포를 따른다고 볼 수 있다. 일반화 선형 모형은 반응 변수가 포아송 분포를 따르는 경우에 적합한 모형을 제공한다. 이번 포스팅에서는 포아송 분포를 따르는 반응 변수에 대하여 모형을 적합하는 방법을 소개한다. 여기서 다루는 내용은 다음과 같다.
1. 모형 적합 알고리즘 유도
먼저 모형 적합 알고리즘에 일반적인 내용을 다룬 포스팅이 있으니 반드시 읽어보기 바란다.
우리에게 데이터 $(\tilde{x}_i, y_i), \;i=1,\ldots,n$이 주어졌다. 여기서 $\tilde{x}_i=(1, x_{i1}, x_{i2}, \ldots, x_{ip})^t$이다. 이제 Fisher Scoring 알고리즘을 유도해보자.
1. 모형 설정
$y_i$가 특정 기간동의 사건 발생 횟수라 하면 $y_i$는 평균이 $\mu_i$인 포아송 분포 $Poisson(\mu_i)$를 따른다고 할 수 있다. 여기서 추정해야 하는 모수는 $\mu_i$이고 연결 함수 $g$를 이용하여 다음과 같이 모형을 설정한다.
$$g(\mu_i) = \beta^t\tilde{x}_i\tag{1}$$
2. 우도 함수(Likelihood Function)
$$\begin{align}f(y_i | \mu_i) &= \frac{e^{-\mu_i}\mu_i^{y_i}}{y_i!} \\ \\ &= \exp [y_i\log\mu_i -\mu_i-\log (y_i!)] \\ \\ &= \exp \{[y_i\theta_i-b(\theta_i)]/a(\phi)+c(y_i, \phi)\}\end{align}$$
여기서 $\theta_i = \log \mu_i$, $b(\theta_i) = \exp (\theta_i)$, $a(\phi)=1$ 그리고 $c(y_i, \phi) = -\log (y_i!)$이다. 따라서 $y_i$의 확률분포는 exponential dispersion family가 된다.
3. 연결 함수(Link Function)
$y_i$의 분포가 exponetial dispersion family이므로 canonical 연결 함수를 쓰자.
$$\eta_i=g(\mu_i) = \theta_i=\log \mu_i$$
연결 함수는 log 함수이다. 이러한 점 때문에 포아송 분포를 따르는 반응 변수의 GLM을 Poisson Log Linear Model(포아송 로그 선형 모형)이라고 한다.
4. 알고리즘
모형 행렬(Model Matrix)을 $X=(x_{ij})$라 하자. $i$번째 대각 원소가 $d_i=\partial \mu_i/\partial \eta_i$인 대각 행렬 $D$, $i$번째 대각 원소가 $w_i=\frac{(\partial \mu_i/\partial\eta_i)^2}{\text{var}(y_i)}$인 대각 행렬 $W$를 계산해야 한다.
$\eta_i = \log\mu_i$임을 이용하면 $d_i$를 다음과 같이 계산할 수 있다.
$$d_i = \frac{\partial \mu_i}{\partial \eta_i} = \left(\frac{\partial \eta_i}{\partial \mu_i}\right)^{-1} = \mu_i$$
$\text{var}(y_i)=\mu_i$임을 이용하면 $w_i$를 다음과 같이 계산할 수 있다.
$$w_i=\frac{(\partial \mu_i/\partial\eta_i)^2}{\text{var}(y_i)}=\frac{\mu_i^2}{\mu_i}=\mu_i$$
이제 Fisher Scoring 알고리즘을 유도할 수 있다(연결 함수가 canonical인 경우 이는 Newton-Raphson 알고리즘과 동일하다).
$\beta^{(k)}$가 주어졌을 때 $\mu_i$의 추정값 $\mu_i^{(k)}$는 (1)과 연결 함수 $g$를 이용하여 다음과 같이 계산할 수 있다.
$$\mu_i^{(k)} = \exp(\beta^{(k)t}\tilde{x}_i)\tag{2}$$
이제 식 (2)를 이용하여 $W$ 와 $D$에서 $\mu_i$를 $\mu_i^{(k)}$로 대체한 행렬을 각각 $W^{(k)}$ 와 $D^{(k)}$라고 하자. 이제 다음과 같이 새로운 반응 변수 $z^{(k)}$를 정의한다.
$$z^{(k)}=X\beta^{(k)}+(D^{(k)})^{-1}(\textbf{y}-\pmb{\mu})$$
이로부터 $\beta^{(k+1)}$을 다음과 같이 구할 수 있다.
$$\beta^{(k+1)} = (X^tW^{(k)}X)^{-1}X^tW^{(k)}z^{(k)}$$
2. 실제 데이터 적용
실제 데이터를 이용하여 Poisson Log Linear Model을 적합해보자. 이번 포스팅에서 사용할 데이터를 다운받는다.
필요한 모듈을 임포트하고 데이터를 불러오자.
import pandas as pd
import numpy as np
import statsmodels.api as sm
from functools import reduce
from scipy.stats import norm
df = pd.read_csv('./num_awards.csv')
데이터 분석 목적은 고등학생을 대상으로 수학 성적과 참여한 프로그램 유형에 따라서 1년에 상을 받은 횟수가 어떻게 변하는지 확인하고자 했다. 모형을 적합하기 전에 간단한 데이터 전처리가 필요하다.
먼저 필요한 열만 추출하고 'prog' 변수는 범주형 변수이므로 이를 더미 변수로 바꾸어주었다. 또한 1로 이루어진 절편항도 추가하였다.
## 전처리
df = df[['num_awards','prog','math']] ## 열 추출
df = pd.get_dummies(df,columns=['prog'],drop_first=True) ## 더미변수 변환
df = sm.add_constant(df) ## 절편항 추가
다음은 모형 적합 과정이다. 반복 횟수가 100을 넘거나 회귀 계수의 현재값과 다음 값의 차이를 계산하고 이 차이값이 0.0001보다 작으면 반복문을 멈추게 하였다.
y = df['num_awards']
X = np.array(df[['const','math','prog_2','prog_3']])
## initial beta
y_mean = y.mean()*np.ones(len(y))
X_tX = np.matmul(X.T,X)
X_tX_inv = np.linalg.inv(X_tX)
beta = reduce(np.dot,[X_tX_inv, X.T, np.log(y_mean)])
# beta = [0]*X.shape[1] ## 베타의 초기값을 0으로 해도 좋다.
epsilon = 0.0001 ## tolerance
max_iter = 100
iter_count = 0
while iter_count <= max_iter:
iter_count += 1
eta = np.matmul(X,beta)
mu = np.exp(eta)
D = np.diag(mu)
D_inv = np.linalg.inv(D)
z = eta + np.matmul(D_inv,y-mu)
W = np.diag(mu)
X_tWX = np.matmul(np.matmul(X.transpose(),W),X)
X_tWX_inv = np.linalg.inv(X_tWX)
X_tWz = np.matmul(np.matmul(X.transpose(),W),z)
next_beta = np.matmul(X_tWX_inv,X_tWz)
diff = np.linalg.norm(next_beta-beta)
beta = next_beta
if diff < epsilon:
break
최종 회귀 계수 추정값을 출력해보자.
각 회귀 계수에 대한 95% 신뢰구간을 구해보자.
## 신뢰구간
alpha = 0.05
columns = ['const','math','prog_2','prog_3']
for i in range(len(beta)):
upper_limit = beta[i] + norm.ppf(1-alpha/2)*np.sqrt(X_tWX_inv[i][i])
lower_limit = beta[i] - norm.ppf(1-alpha/2)*np.sqrt(X_tWX_inv[i][i])
print(columns[i], lower_limit, upper_limit)
prog_3에 대한 신뢰구간이 0을 포함하므로 이에 대응하는 회귀 계수 추정값은 유의하지 않다고 판단할 수 있다. 즉, Vocational 프로그램은 수상 횟수에 영향을 크게 주지 않는 것 같다.
이번엔 statsmodels를 이용하여 Poisson Log Linear 모형을 적합해보자. 아래의 모듈을 추가적으로 임포트해준다.
from statsmodels.genmod.generalized_linear_model import GLM
from statsmodels.genmod.families import Poisson
모형을 적합하기 위해서 GLM의 첫 번째 인자와 두 번째 인자에 반응 변수와 설명변수를 넣고 family 인자에 Poisson()을 넣어준다. 그러고 나서 fit을 통하여 모형을 적합하면 된다.
y = df['num_awards']
X = df[['const','math','prog_2','prog_3']]
model = GLM(y,X,family=Poisson())
results = model.fit()
모형 적합 결과는 summary를 이용하여 확인할 수 있다.
results.summary()
앞에서 본 것과 동일한 회귀 계수 추정값과 신뢰구간을 얻을 수 있다.
참고자료
Alen Agresti - Foundations of Linear and Generalized Linear Models
댓글