응용통계학 (13주차) 5월 26일
GLM, 일반화선형모형 3번째
선형모형, 로지스틱모형, count regression 등에 이르기까지 다양한 종류의 response를 (1) response 자체의 분포의 특성을 고려하여 (2) conditional mean function을 추정하기 위한 적절한 연결함수를 고려하여 모형화 및 추정을 시도하였다.
이러한 것들을 보다 일반화하여 앞서 다루었던 모형 및 추정법등을 포함하는 하나의 모형을 생각할 수 있는데 이를 일반화선형모형 (Generalized linear model : GLM)이라 한다.
GLM은 크게 두 가지 요소를 특정하는 것으로 정의된. 첫째 response의 분포로, 기본적으로 지수족 (exponential family)을 따르는 것으로 가정한다. 둘째는 link 함수로, response의 conditional mean function을 linear predictor와 어떻게 연결할 것인지이다.
- binary - logistic link function
- poisson - log link function
다음과 같은 분포를 따를 때 $Y$가 지수족을 따른다고 한다.
$$f(y|\theta,\phi) = exp \big[ \frac{y\theta - b(\theta)}{a(\phi)} + c(y,\theta) \big]$$ $\theta$ : canonical parameter, $\phi$: disperison parameter
- 정규분포($N(\mu, \sigma^2)$): $$\theta = \mu, \phi = \sigma^2, a(\phi) = \phi, b(\theta) = \theta^2, c(y,\phi) = -(y^2 /\phi + log(2\pi\phi))/2$$
- $b'(\theta) = \theta = \mu$
- $f(\theta) = b'(\theta) = \theta$
- $f^{-1}(\theta) = \theta = g(\theta)$ : 정준 연결 함수 $\to$ 안 쓰는게 자연스런 $\theta$ function 쓰는 것!
- Poisson: $$\theta = log(\mu), \phi = 1, a(\phi) = 1, b(\theta) = exp(\theta), c(y,\phi) = -logy!$$
- $b'(\theta) = exp(\theta) = exp(log(\mu)= \mu$
- $ f(\theta)= b'(\theta) = e^{\theta}$
- $f^{-1}(\theta) = log(\theta) = g(\theta)$ : 정준 연결 함수 $\to$ $log$쓰는 게 자연스런 $\theta$ function 쓰는 법
- Binomial: $$\theta = log(\mu/(1-\mu)), b(\theta) = -nlog(1-\mu) = nlog(1+exp(\theta)), c(y,\phi) = log\begin{pmatrix} n\\ y\end{pmatrix}$$
- $\mu = P(Y=1)$
- $b'(\theta) = \frac{exp(\theta)}{1+exp(\theta)}$
- $f(\theta) = \frac{exp(\theta)}{1+exp(\theta)}$
- $f^{-1} = log\frac{exp(\theta)}{1+exp(\theta)} = g(\theta)$ : 정준 연결 함수 $\to$ $logit$ 쓰는 게 자연스런 $\theta$ function 쓰는 법
감마 및 역감마분포도 지수족에 포함된다.
일반적으로 지수족을 따르는 경우 response의 평균과 분산이 어떻게 표현되는지는 다음으로부터 이끌어 낼 수 있다.
먼저 log-likelihood for a sigle $y$는 다음과 같이 주어지므로:
$$l(\theta) = (y\theta)-b(\theta))/a(\phi) + c(y,\phi)$$
이 함수의 일차/이차 미분은 다음과 같이 계산된다.
$$l'(\theta) = (y-b'(\theta))/a(\phi), l''(\theta) = -b''(\theta)/a(\phi)$$
- 분자가 0이 되면 0이 되므로 $y=b'(\theta)$만 만족하면 된다.
그러면 score function 이 mean이 0이라는 성질과 이차미분에 대한 성질로부터 ;
$El'(\theta) = 0$과 $El''(\theta) = -E[l'(\theta)^2]$을 얻고 이로부터 다음 식을 얻을 수 있다.
$$E(Y) = \mu b'(\theta), Var(Y) = b''(\theta)a(\phi)$$
- $E(Y)$에서는 $\theta$만 역할을 하지만, $Var(Y)$에서는 $\theta$랑 $\phi$ 모두 역할이 있다.
위로부터 평균은 $θ$이 함수임을 알 수 있다. 또한 분산은 $θ$와 $ϕ$의 함수의 곱 형태로 주어진다. 여기서 $V(μ)=b′′(θ)/w$를 variance function이라 지칭하는데 이는 분산이 평균과 어떻게 연결되는지를 기술해주는 것으로 이해할 수 있다.
$$\eta = g(\mu) = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p = x^\top \beta$$
- $g(\mu) = \theta$
- $\to$ $g(b'(\theta) = \theta$
- $\to$ $(g \bullet b')(\theta) = \theta$
- $\to$ $g(\theta) = b'^{-1}(\theta)$
-
- 분포가 정해진 후
-
- link function 이 정해진다.
-
$E(Y) = \mu$ 이다. 원론적으로는 어떤 monotone continuous and differentiable function도 연결함수의 역할을 할 수 있으나 GLM에서 보통 선택하는 연결함수에는 몇 가지 유용한 성질을 가지게 된다. 이를 정준연결함수(canonical link function)이라 한다.
정준연결함수 $g$는 $η=g(μ)=θ$를 만족하는 함수를 말한다. 이는 곧 $g(b′(θ))=θ$를 함의하게 된다. Poisson 분포에서는 log함수가 Binomial 분포에서는 logit 함수가 정규분포에서는 항등함수가 정준연결함수가 됨이 알려져 있다.
이러한 정준연결함수는 수학적으로 또한 추정량 계산상에 있어서 여러가지 장점이 있어 널리 쓰인다.
- 적어도 이론적으로 추정량 찾는데 많은 장점이 있다.
최대가능도추정법이 GLM의 추정에서 기본적으로 쓰인다. 분포가 정규분포일 때를 제외하면 가능도방정식의 해가 명시적으로 표현되지 않기 때문에 수치적인 접근법을 이용해야 한다. 보통 Newton-Raphson에 근거한 방법들을 이용하게 된다.
- Null model : 가장 간단한 모형으로 설명변수가 하나도 포함되지 않은 모형
- Full (Saturated) model :가장 복잡하고 큰 모형으로 $n$개의 관측치에 대하여 $n$개의 parameter를 사용. 이 경우 $\hat{μ}=y$가 만족되게 된다.
앞서 언급했듯이 현재 모형이 데이터를 얼마나 잘 적합하는지를 측정하고 할 때, 현재 모형을 full model과 비교하게 된다. 즉, 다음과 같은 두 모형의 차이를 생각한다.
$$1[log(y,\phi|y) - log(\hat{\mu},\phi|y)]$$
- 적합도가 좋으면 이 값이 작아진다.
이를 deviance라 하고 다음과 같이 쓴다.
$$D(y,\hat{\mu})$$
- $\sqrt{d_i}$ : Deviance type residual
| GLM | Deviance |
|---|---|
| Gaussian_Normal | $\sum_i (y_i - \hat{\mu}_i)^2$ $\approx$ SSE |
| Poisson | $2\sum_i[y_i log(y_i/\hat{\mu}-i) - (y_i - \hat{\mu}_i)]$ $ =2 생략하고 \sum d_i = \sum (\sqrt{d_i})^2$,제곱하는 개념으로 보고 root씌운 것을 잔차라고 놓는다. $\sqrt{d_i}$는 양수/음수 가능 |
| Binomial | $2\sum_i [y_i(log(y_i/\hat{\mu}_i) + (m-y_i)log((m-y_i)/(m-\hat{\mu}_i))]$ |
| Gamma | $2\sum_i[-log(y_i/\hat{\mu}_i) + (y_i - \hat{\mu}_i)/\hat{\mu}_i]]$ |
| Inverse Gaussian | $\sum_i (y_i - \hat{\mu}_i)^2/(\hat{\mu}^{2}_{i} y_i)$ |
이는 선형모형에거 SSE를 GLM 하에서 일반화시킨 것으로 간주할 수 있다.
적합도를 측정하는 또 다른 방법 중 하나는 Pearson’s $X^2$ 통계량을 고려하는 것이다.
$$X^2 = \sum_i \frac{(y_i - \hat{\mu_i})^2}{V(\hat{\mu}_i)}$$
- 일종의 잔차 표준화이다.
- 표준화 $\big(\frac{X-\mu}{\sigma}\big)^2 = \frac{(X-\mu)^2}{\sigma} \sim \chi^2_{(1)}, Z \sim \chi^2_{(1)}$
- 값이 클 수록 적합도가 안 좋다는 뜻, 값이 작을수록 적합도가 좋다는 뜻
더 큰 모형을 $Ω$ 작은 모형을 $ω$라 할 때, scaled deviance의 차이 $D_ω−D_Ω$는 점근적으로 $χ^2$ 분포를library(faraway) 따르게 된다. 이 때 자유도는 두 모형의 모수 개수의 차이이다.
- 적합도 평가 뿐만이 아니라 차이의 카이제곱분포를 따름을 보이고 유의성 평가까지 할 것.
library(faraway)
- conc는 살충제 농도이고,
- binomial로 놓은 이유는 생존을 이진수로 나타내니까
data(bliss, package="faraway")
modl <- glm(cbind(dead,alive) ~ conc, family=binomial, bliss)
summary(modl)
- conc의 estimate가 양수이다. 1.1619
1-pchisq(deviance(modl),df.residual(modl)) # test for goodness of fit
- 함수 2개 쓰는 2차 함수 사용
- test='chi' 없어도 된대
modl2 <- glm(cbind(dead,alive) ~ conc+I(conc^2), family=binomial, bliss)
anova(modl,modl2,test="Chi") # test for additional parameter
anova(modl2,test="Chi") # test for additional parameter
- 0.6685914로 유의하지 않음.
- 유의미하게 적합도를 높이지 않는다.
- 굳이 모델 2로 확장하지 않아도 된다.
- 모델 2가 안 좋은 모델은 아님
선형모형에서 잔차를 중심으로 모형진단을 했듯이 GLM에서도 비슷한 절차가 필요할 수 있다.
- 이 모형 진단은 중요도가 떨어지긴 한다.
잔차는 response 상에서 관측치와 모형에 의한 적합치 사이의 차이를 의미한다. 선형모형에서 정의했던 것을 확장하여 잔차를 GLM에서 정의하기도 한다.
- 사실 로지스틱하고 진단을 잘 안 하긴 하다.
선형모형에서의 표준화잔차와 유사한 개념으로 이해할 수 있다.
$$r_P = \frac{y-\hat{\mu}}{\sqrt{V(\mu})}$$
- 피어슨 타입 잔차
여기서 $V(\mu)=b''(\theta)$이다. 또한 $\sum r^{2}_{P}=X^2$가 된다.
- $y = X\beta' + \epsilon$에서 $\epsilon \sim N(,)$임을 확닉하나? $|to$ $y$랑 분포 맞추자.
- 연속 $y \to \sqrt{y}. log y$
- 이진 $y \to $ ? $\epsilon$(Bernoulli) or $\epsilon$(Binary)
- count $y \to $ ? $\epsilon$(poisson)
Deviance로부터 정의되는 타입의 잔차이다.
- Deviance는 SSE의 일반화 버전
- $\sum^{n}_{i=1} e^{2}_{i}$ 잔차제곱합 일반화 = 총합의 개념 $y-\hat{\mu}>0$
$$sum r^{2}_{D} = Deviance = \sum d_i$$
- $d_i$를 Deviance의 제곱이라 보자
따라서
$$r_D = sign(y - \hat{\mu})\sqrt{d_i}$$
- $y>\hat{\mu} \to$ residual >0
- $y<\hat{\mu} \to$ residual <0
residuals(modl) # deviance residual
residuals(modl,"pearson") # Pearson residual
residuals(modl,"response") # response residual
bliss$dead/30 - fitted(modl)
QQ plot은 선형모형에서 정규성을 확인하기 위해 사용하는 시각화 방법 중 하나였다. GLM에서는 분포 자체를 확인하기 위한 목적보다는 이상점의 존재를 확인하기 위한 목적으로 보통 활용된다.