Exploring Regression Models for Regression Analysis (2): GLM, Exponential Family, Link Function, IRLS(Fisher scoring), Cluster-robust standard error

들어가며

2장에서는 1장에서 다룬 기본 linear regression에서 link function을 도입하여 regression의 개념을 outcome of single yes/no, outcome of single K-way, count 등 non-normal한 종속변수로 확장한 Generalized linear model의 개념을 Exponential Family, Link Function와 같은 핵심개념과 함께 깊게 살펴보며, model의 parameter를 estimate하는 알고리즘인 IRLS(Fisher scoring) 집중적으로 estimation methods를 소개합니다. 이후, HC standard errors의 clustered data 버전인 Cluster-robust standard error를 보고 마치겠습니다.

1. Generalized Linear Models (GLMs)

1.1. Linear Model 한계

---

1장에서 본 Linear Regression Model은 (1) 선형성(Linearity) 가정, (2) 오차항의 정규성(Normality) 가정, (3) 오차항의 독립성(Independence) 가정, (4) 오차항의 등분산성(Homoscedasticity) 가정에서 비롯된 모델이었고, Heteroskedasticity-Consistent Standard Errors (HC Standard Errors)를 통해 오차항의 등분산성(Homoscedasticity) 가정이 깨진 data에 대해서도 Linear model로부터 얻은 모델 parameter의 분산을 robust하게 추정할 수 있었습니다. 그러나, 위에서 언급하였듯 outcome of single yes/no, outcome of single K-way, count 등 많은 data는 반응 변수 Y가 정규분포를 따르지 않거나 등분산성 가정, 선형성에 위배됩니다. 각각에 대해서 좀 더 설명하자면, 어떤 사건이나 행동이 일어나거나 그렇지 않은 경우를 고려하는 이진 데이터(binary data, outcome of single yes/no)의 경우, $Y\in \{0,1\}$로 제한되며 이를 $Y\in \mathbb{R}$인 정규분포로 가정하는 것은 옳지 않습니다. 특정 기간 동안 발생하는 사건의 횟수 등, 이진 분류처럼 discrete한 종속변수 값을 가지는 카운트 데이터(count data) 또한 discrete(정수) 값만 갖으며, 이 두 경우는 종종 분산이 평균(모델의 예측)에 비례하는 형태를 갖을 수 있고, 이는 당연하게도 등분산성 가정을 위배합니다.

이러한 데이터의 경우 단순히 독립변수의 선형결합 형태, 또는 기하학적으로는 Hyper plane 형태로 모델을 만들면, 비선형적인 (이진 데이터 등) 위 같은 경우에 대해서는 올바르게 고려하지 못할 것입니다. 이러한 기존의 Linear Regression 모델의 한계를 극복하고, (종속변수의) 다양한 형태의 데이터를 모델링하기 위해 여러 함수를 설계함으로써 유연성을 확장한 Generalized Linear Models이 개발되었습니다. Generalized Linear Models(GLMs)의 중요 구성 요소들과 원리를 간략히 설명해보자면, 선형 결합으로 바로 종속변수를 예측하는 대신, non-linear한 Link Function에 넣어 최종적으로 예측함으로써 non-linear한 종속변수에도 fit 할 수 있고, 이에 따라 종속변수의 분포가 정규분포가 아닌 다른 분포(Exponential Family)도 포함할 수 있도록 하였으며, 이 Exponential Family와 Variance function구성은 종속변수의 분산이 모델의 예측값(종속변수의 mean)마다 다를 수 있도록 합니다. 이를 통해 Generalized Linear Models는 위 네 개의 Linear Regression 가정 중 (1) 선형성(Linearity) 가정, (2) 오차항의 정규성(Normality) 가정, (4) 오차항의 등분산성(Homoscedasticity) 가정을 깼으며, 위에서 Linear Model의 한계로 언급한 데이터들을 고려할 수 있는 모델입니다.

1.2. GLM 정의 및 수학적 표현

---

GLM은 세 가지 구성 요소 (Random component, Systematic component, Link function)으로 정의 되며, 이때 Random component는 Y를 Exponential Family로, Systematic component는 Linear predictor와 Link function으로 구성됩니다. 어떻게 Generalized Linear Models가 설계되었는지 component들을 하나하나 자세히 다뤄보겠습니다.

Linear predictor

---

Linear predictor $\eta$는 말그대로 Linear Model처럼 모델 parameter $\mathit{\beta }$와 독립변수 $\mathbf{X}$의 선형 결합으로, 기존에는 $\eta$로 바로 $\mathbf{y}$를 추정하여 non-linear한 종속변수를 고려하지 못하였었다면, GLM은 $\eta$를 계산한 후, 이 값을 non-linear한 Link function에 input하여 최종적으로 종속변수를 예측합니다. 중요한 점은, 이는 단순히 종속변수를 transform한 뒤(로그 등) 이전처럼 선형적으로 추정하는 Transformation (with LM)과 다르다는 점입니다. 가장 큰 차이점은 Transformation을 함으로써 종속변수의 sample space에서 boundaries에 있는 값들은 정의가 되지 않고(로그는 0에서 정의되지 않음.), 이후 바로 Linear Model을 사용하기 위해선 종속변수가 변형 이후 반드시 linearity와 variance의 homogeneity가 거의 보장되어야 합니다.(기존 LM을 사용하기 때문에 이때 사용한 가정 또한 필요하게 됩니다.) GLM의 Linear predictor (선형 예측자) $\eta$의 식은 다음과 같습니다: $${\eta }_i={\beta }_0+{\beta }_1{x}_{1i}+{\beta }_2{x}_{2i}+\cdots +{\beta }_p{x}_{pi}$$

Link function (링크 함수)

---

Link function $g({\mu }_i)$는 non-linear하고 미분 및 inverse(역)가 가능한 함수로 정의되며, 종속변수의 평균 $E(Y_i)={\mu }_i$를 선형 예측자 ${\eta }_i$와 연결하여 간접적으로 독립변수 및 모델 parameter의 선형결합과 종속변수를 mapping하는 역할을 합니다.
$$g({\mu }_i)={\eta }_i$$

Variance function (분산 함수)

---

분산 함수는 평균 ${\mu }_i$에 따라 종속변수의 분산이 어떻게 변하는지를 나타냅니다. 이를 통해 간접적으로 독립변수에 따라 분산이 다르게 나오는 것을 반영할 수 있으며 식은 아래와 같고,
$$\mathrm{Var}(Y_i)=\varphi V({\mu }_i)$$

여기서 $\varphi$는 dispersion parameter로, 일반적으로 특정 분포에 따라 다르게 정의됩니다. (예: Poisson 분포에서는 $\varphi =1$).

Exponential Family

---

GLM은 종속변수의 분포로 Gaussian(또는 정규분포)를 포함한, 더욱 general한 Exponential Family을 고려하며, 이 분포는 다음과 같은 일반 형태를 가집니다.

$$f(y;\theta ,\varphi )=\exp \{\frac{y\theta -b(\theta )}{\varphi }+c(y,\varphi )\}$$

즉 GLM은 LM과 다르게 linear predictor, link function, variance function을 설계함으로써 종속변수가 더욱 general한 분포인 exponential family distribution인 경우에도 잘 mapping할 수 있도록 하는 모델이라고 볼 수 있습니다. 위 식에서 의미론적으로 각 parameters를 해석하면 $\theta$는 canonical parameter로 분포의 위치를 나타내는 파라미터, $\varphi$는 dispersion parameter로 분산과 관련된 파라미터, $b(\theta )$는 평균과 분산 관계를 정의하는 함수입니다. 이 분포에 대해 $E(Y)=b^{\prime }(\theta )=\mu$, $\mathrm{var}(Y)=\varphi b^{″}(\theta )=\varphi V(\mu )$이라는 특성이 증명 가능하고, 이는 “2. GLMs 추정”에서 모델 $\beta$를 추정하는 과정에 필요하기 때문에 아래에서 증명할 것입니다. 이보다 더욱 general한 분포로 (dispersion parameter 관련) exponential dispersion family가 있습니다.

다음으로 넘어가기 전에 간단하게 잘 알려져있는 Exponential Family의 예시인 정규분포, 이항분포, 포아송분포, 감마분포가 이에 포함됨을 확인해보겠습니다.

(1) 정규분포 (Normal Distribution)

정규분포의 확률밀도함수는 다음과 같습니다: $$f(y;\mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \{-\frac{(y-\mu {)}^2}{2{\sigma }^2}\}.$$ 이를 Exponential Family 형태로 변환하면: $$f(y;\mu ,{\sigma }^2)=\exp \{\frac{y\mu -\frac{{\mu }^2}{2}}{{\sigma }^2}-\frac{y^2}{2{\sigma }^2}-\frac{1}{2}\log (2\pi {\sigma }^2)\}.$$

• Canonical parameter: $\theta =\mu$,
• Dispersion parameter: $\varphi ={\sigma }^2$,
• $b(\theta )=\frac{{\theta }^2}{2}$,
• $b^{\prime }(\theta )=\theta =\mu$,
• $b^{″}(\theta )=1$,
• $c(y,\varphi )=-\frac{y^2}{2\varphi }-\frac{1}{2}\log (2\pi \varphi )$.

(2) 이항분포 (Binomial Distribution)

이항분포의 확률질량함수는 다음과 같습니다: $$f(y;n,p)=(\genfrac{}{}{0ex}{}{n}{y})p^y(1-p{)}^{n-y}.$$ 이를 Exponential Family 형태로 변환하면: $$f(y;n,p)=\exp \{y\log \left(\frac{p}{1-p}\right)+n\log (1-p)+\log (\genfrac{}{}{0ex}{}{n}{y})\}.$$

• Canonical parameter: $\theta =\log \left(\frac{p}{1-p}\right)$,
• Dispersion parameter: $\varphi =1$,
• $b(\theta )=n\log (1+e^{\theta })$,
• $b^{\prime }(\theta )=e^{\theta }=\lambda$,
• $b^{″}(\theta )=e^{\theta }=\lambda$,
• $c(y,\varphi )=\log (\genfrac{}{}{0ex}{}{n}{y})$.

(3) 포아송분포 (Poisson Distribution)

포아송분포의 확률질량함수는 다음과 같습니다: $$f(y;\lambda )=\frac{{\lambda }^y{e}^{-\lambda }}{y!}.$$ 이를 Exponential Family 형태로 변환하면: $$f(y;\lambda )=\exp \{y\log \lambda -\lambda -\log (y!)\}.$$

• Canonical parameter: $\theta =\log \lambda$,
• Dispersion parameter: $\varphi =1$,
• $b(\theta )=e^{\theta }$,
• $b^{\prime }(\theta )=\frac{n{e}^{\theta }}{1+e^{\theta }}=np$,
• $b^{″}(\theta )=\frac{n{e}^{\theta }}{(1+e^{\theta }{)}^2}=np(1-p)$,
• $c(y,\varphi )=-\log (y!)$.

(4) 감마분포 (Gamma Distribution)

감마분포의 확률밀도함수는 다음과 같습니다: $$f(y;k,\theta )=\frac{1}{\mathrm{\Gamma }(k){\theta }^k}{y}^{k-1}{e}^{-\frac{y}{\theta }}.$$ 이를 Exponential Family 형태로 변환하면: $$f(y;k,\theta )=\exp \{-\frac{y}{\theta }+(k-1)\log y-k\log \theta -\log \mathrm{\Gamma }(k)\}.$$

• Canonical parameter: $\theta =-\frac{1}{\theta }$,
• Dispersion parameter: $\varphi =\frac{1}{k}$,
• $b(\theta )=-\log (-\theta )$,
• $b^{\prime }(\theta )=\frac{1}{\theta }=\mu$,
• $b^{″}(\theta )=\frac{1}{{\theta }^2}={\mu }^2$,
• $c(y,\varphi )=(k-1)\log y-\log \mathrm{\Gamma }(k)$.

위에서 $b(\theta )$는 한 번 미분하면 mean, 두 번 미분하면 variance의 term과 관련됨을 언급하였고 위 4개의 분포에서 원래 알고 계신 mean, variance와 $b^{\prime }(\theta )$, $b^{″}(\theta )$가 dispersion parameter를 고려하면 일치한 것을 확인하실 수 있습니다. 이는 cumulant generating function의 일부이기 때문이며, 따라서 평균과 분산 관계를 정의하는 항이라고 언급하였던 것입니다.

Canonical Link

---

Canonical link는 GLM에서 통계적 성질을 최적화하기 위해 사용되는 링크 함수(link function)로, 다음과 같이 정의됩니다:

$$g({\mu }_i)=g(b^{\prime }({\theta }_i))={\theta }_i={\eta }_i$$ 이 식의 의미는 결국 아래 식과 같습니다. $$g=(b^{\prime }{)}^{-1}$$

아래에서 확인하겠지만, Binomial 분포의 경우 canonical link는 logit 함수이고, Poisson 분포의 경우 canonical link는 log 함수이며, Canonical link를 사용하면 MLE(Maximum Likelihood Estimation) 과정이 단순화되고, efficient한 추정치를 얻을 수 있기 때문에 link function은 거의 항상 Canonical link로 정의합니다. 또한, canonical하지 않은 link function의 경우에도 위 Exponential Family distribution에서 식 조작을 통해 canonical link 형태를 만들 수 있습니다.

1.3. GLM 예시

---

위 철학에 따라, data가 따르는 Exponential Family 중 특정 분포가 정해지면, 이에 해당하는 보통 사용하는 Link function(Canonical link), Variance function가 정해져 있고 결국 모델이 특정되며, GLM은 이렇게 특정될 수 있는 모든 모델에서 공통적으로 parameter와 그 variance를 추정해내는 general한 모델이라고 생각할 수 있습니다. 여기에서 다루지는 않겠지만, 사실 특정한 형태의 data에서 가능한 link function은 여러 개이며, 이에 따라 variance function도 여러 가지가 가능할 수 있습니다. 그러나 효율성과 computation cost를 고려하여 보통 사용되는 function forms는 정해져 있다고 알아두시면 좋을 것 같습니다. 아래 예시 중 대표적으로 Binomial 예시에서는 link function이 0 이상 1 이하의 정의역에서 실수 전체(for linear predictor)를 map할 수 있는 미분 및 역이 가능한 함수이면 되지만, 보통 logit function이 사용됩니다. 헷갈릴 수 있지만 아래 Exponential Family 중 친숙한 분포의 예시를 직관적인 관점에서 고려하여 위에서 얻은 Exponential Family의 form과 같은 결과가 나옴을 보시면 좋을 것 같습니다.

아래의 예시에서 linear predictor는 공통이므로 미리 정의하고 각각의 link function, variance function은 어떻게 특정되는지를 보겠습니다.

$${\eta }_i={\beta }_0+{\beta }_1{x}_{1i}+\cdots +{\beta }_p{x}_{pi}$$

where ${\eta }_i$는 linear predictor, ${\beta }_0,{\beta }_1,...,{\beta }_p$는 regression coefficients(parameters), $x_{1i},...,x_{pi}$는 predictor variables(독립변수) 입니다.

Binomial Case

---

Binomial Data, 즉 종속변수가 $Y_i\sim \text{Binomial}(n_i,p_i)$인 data의 경우, 종속변수의 기댓값의 sample space 또한 0~1이며, 우리가 모델링하고 싶은 값이 $Y_i/n_i$인 경우를 상정해보겠습니다. 직관적으로 의미로부터 functions가 어떻게 되어야 할 지 생각해보면, $E(Y_i/n_i)=p_i$이고, 분산은 $\frac{1}{n_i}{p}_i(1-p_i)$ 입니다. $Y_i/n_i$의 variance 식에 $Y_i/n_i$의 mean이 들어감을 알 수 있고, 기존 LM에서는 이렇게 관측치에 따라 다르게 variance를 고려할 수 없었지만, GLM에서는 이 관계를 variance function을 통해 고려할 수 있으며, 식은 다음과 같습니다: $$V({\mu }_i)={\mu }_i(1-{\mu }_i)$$ 또한, $Y_i/n_i$와 linear predictor를 matching 해줄 수 있는 미분가능한 function을 link로 고려해야 하고, Binomial에서는 non-linear link function으로 대부분 logit 함수를 사용합니다. (이때, logit function의 inverse는 sigmoid function입니다.)

$$g({\mu }_i)=\log \left(\frac{{\mu }_i}{1-{\mu }_i}\right)$$

이 식은 위에서 확인한 이항분포의 canonical parameter와 같은 형태임을 알 수 있습니다.

Poisson Case

---

Poisson Data, 즉 종속변수가 $Y_i\sim \text{Poisson}({\lambda }_i)$인 data이고 우리가 모델링하고 싶은 값이 종속변수 $Y_i$인 경우를 상정해보겠습니다. 직관적으로 의미로부터 functions가 어떻게 되어야 할 지 생각해보면, $E(Y_i)={\lambda }_i$이고 분산은 ${\lambda }_i$ 이므로, 마찬가지로 $Y_i$의 variance 식에 $Y_i$의 mean이 들어감을 알 수 있고, variance function은 다음과 같습니다:

$$V({\mu }_i)={\mu }_i$$ 또한, $Y_i$의 sample space는 0 이상의 실수로, 이와 linear predictor를 matching 해줄 수 있는 미분가능한 non-linear function을 link로 고려해야 하고, Binomial에서는 link function으로 대부분 log 함수를 사용합니다. (inverse는 지수 함수.) $$g({\mu }_i)=\log ({\mu }_i)$$이 식은 위에서 확인한 포아송분포의 canonical parameter와 같은 형태임을 알 수 있습니다. 위 두 예시에서는 어떻게 GLMs의 구성 요소들이 선택되는지 직관적으로 보았고, 이는 이해를 돕기 위한 해석이었으며 이미 위에서 Exponential Family에 포함됨을 보일 때 같은 결과가 나왔다는 것을 보시면 됩니다. 위 Binomial Data의 GLM은 Logistic Regression, Poisson Data의 GLM은 Poisson Regression으로도 불립니다.

2. GLMs 추정

위 내용들을 통해서 GLMs가 어떻게 비정규분포를 갖는 종속변수를 고려해서 잘 작동하며, 어떠한 함수(Link function, Variance function, Exponential Family, Canonical Link)가 어떠한 수식과 철학으로 GLM을 구성하고 있는지 확인할 수 있었습니다. GLMs은 거의 대부분의 고려 가능한 data 분포가 Exponential Family를 따르며, 이에 대해 일관적인 form과 parameter estimation이 가능하기 때문에 아주 powerful한 Regression Model입니다. 그러나 어떻게 Exponential Family를 따르는 data를 다룰 수 있는지는 확인할 수 있었지만, 어떻게 Regression Model의 parameter와 그 분산을 추정할 수 있는지는 다루지 않았습니다. Linear Model에서는 closed-form solution을 쉽게 찾을 수 있었지만, GLM은 대부분의 경우(있는 경우도 있습니다.) 이러한 closed-form이 없어 컴퓨터 프로그램으로 여러 번에 걸쳐 추정할 수 있도록 알고리즘을 구현하여 이를 추정합니다. 실제로 이 estimation의 수식과 실제 구현 과정을 다루기 위해서는 긴 증명 과정을 거치는데, 최대한 중요한 부분은 빠지지 않으면서 증명해보겠습니다. 우선, MLE로 모델을 추정하는 과정을 증명하기 위해 필요한 두 가지 유용한 성질을 살펴보겠습니다. (Derivatives of Log Likelihood’s, Exponential Family 성질)

2.1. Derivatives of Log Likelihood’s 성질

---

확률변수 $Y$의 밀도 함수 $f(y;\theta )$가 주어지며, 이때 $\theta$는 스칼라 매개변수라고 가정하겠습니다. 또, $\ell$이 $\theta$에 대해 최소 두 번 미분 가능하다고 가정하면, 단일 관측치 $Y$에 대한 로그 가능도(log likelihood) 함수 $\ell (\theta ;Y)$에 대해서 함수의 첫 번째 및 두 번째 도함수는 다음과 같습니다.

• 첫 번째 도함수: $${\ell }^{\prime }=\frac{d\ell }{d\theta }$$
• 두 번째 도함수: $${\ell }^{″}=\frac{d^2\ell }{d{\theta }^2}$$

이때, 이 두 함수들은 다음 관계식이 성립합니다.

$$E\{{\ell }^{\prime }(\theta ;Y)\}=0$$

$$E[\{{\ell }^{\prime }(\theta ;Y){\}}^2]=-E\{{\ell }^{″}(\theta ;Y)\}$$

이를 증명해보겠습니다. 의지를 잃지 않기 위해 위 식들의 의미를 스포하자면, MLE를 통해 모델을 추정할 때 보통 log likelihood를 모델의 parameter로 미분한 식이 0(또는 영벡터)이 되도록 하는 parameter를 찾음으로써 이를 수행하는데, 첫 번째 식은 이 미분한 식(score function)의 기댓값(평균)이 0이라는 의미이고, 두 번째 식은 첫 번째 식에서 mean이 0이었으므로 왼쪽항의 제곱 안에 -0을 넣어주면 $$E[\{{\ell }^{\prime }(\theta ;Y)-E[{\ell }^{\prime }(\theta ;Y)]{\}}^2]$$

가 되어 분산 term이 되고, 따라서 분산은 이차 도함수의 기댓값의 음수와 같다는 의미입니다. 이러한 성질들을 이용해서 앞으로 Likelihood 기반의 다양한 모델 추정을 수행할 수 있게 됩니다.

(1) Prove $E\{{\ell }^{\prime }(\theta ;Y)\}=0$.

우선, 확률 분포는 모든 범위에서의 적분 또는 누적합이 1이므로,

$$1=\int f(y;\theta )dy$$

입니다. 이제 양변을 $\theta$에 대해 미분한 후 미분과 적분의 순서를 바꾸면,

$$\begin{gathered} 0=\frac{d}{d\theta }\int f(y;\theta )dy \\ =\int \frac{d}{d\theta }f(y;\theta )dy \end{gathered}$$

입니다. 여기서 수학적 증명 과정에서 굉장히 자주 사용되는 skill $\frac{d}{d\theta }f(y;\theta )=\frac{d}{d\theta }\{\log f(y;\theta )\}f(y;\theta )$를 사용하면

$$\begin{gathered} 0=\int \frac{d}{d\theta }f(y;\theta )dy \\ =\int \frac{d}{d\theta }\{\log f(y;\theta )\}f(y;\theta )dy \\ =\int {\ell }^{\prime }(\theta ;Y)f(y;\theta )dy \\ =E\{{\ell }^{\prime }(\theta ;Y)\} \end{gathered}$$ 입니다. 의미를 다시 해석해보면, 어떠한 분포를 따르는 $Y$와 이의 매개변수 $\theta$에 대해서, 우리는 MLE를 통해 log likelihood 함수 $\ell (\theta ;Y)$를 $\theta$로 미분하였을 때 0이 나오도록 하는 $\hat{\theta }$를 찾음으로써 parameter를 estimate합니다. 위 (1)은 이러한 ${\ell }^{\prime }(\theta ;Y)$의 기댓값은 $Y$의 분포가 이계도함수가 존재한다면 어떤 분포이건 관계 없이 0임을 보인 것입니다.

(2) Prove ${\ell }^{″}=\frac{d^2\ell }{d{\theta }^2}$.

동일한 논리를 따라 위 식을 한 번 더 미분하면,

$$0=\frac{d}{d\theta }[\int \frac{d}{d\theta }\{\log f(Y;\theta )\}f(y;\theta )dy]$$

입니다. 두 함수의 곱 형태의 미분이며 둘 다 $\theta$를 포함하므로 이를 전개하고 마찬가지로 기댓값으로 표기하면,

$$\begin{gathered} 0=\int \frac{d^2}{d{\theta }^2}\{\log f(y;\theta )\}f(y;\theta )dy+\int \frac{d}{d\theta }\{\log f(y;\theta )\}\frac{d}{d\theta }f(y;\theta )dy \\ =E\{{\ell }^{″}(\theta ;Y)\}+E[\{{\ell }^{\prime }(\theta ;Y){\}}^2] \end{gathered}$$

이고, 따라서 아래 식이 증명되었습니다.

$$E[\{{\ell }^{\prime }(\theta ;Y){\}}^2]=-E\{{\ell }^{″}(\theta ;Y)\}$$

위 증명 과정에서 미분과 적분 연산자의 교환을 정당화하는 과정의 설명이 생략되었지만 exponential family에선 문제가 없고, Y가 discrete한 경우는 적분을 누적합으로 바꿔주면 된다고 얘기해두며 마무리 하겠습니다. 또한, 위 증명에서는 $\theta$가 1차원 스칼라 변수라고 가정했지만, 다차원 매개변수에 대해서도 동일한 결과가 성립됩니다. 식의 의미를 마지막으로 되짚어보면, log likelihood의 일차 도함수는 기대값이 0이고, 이 일차 도함수의 공분산 행렬은 이차 도함수 행렬의 기대값의 음수에 해당합니다. 이 값은 피셔 정보 행렬(Fisher Information Matrix)이라고도 불립니다. (함수의 기댓값이라는 말이 어색하게 들릴 수도 있는데, 애초에 모든 랜덤(확률)변수는 어떠한 관측치에 대해서 실수를 output하는 함수임을 되새기면 좋을 것 같습니다.)

2.2 Exponential Family 성질

---

이번에는 위에서 증명한 수식을 통해서 Exponential Family를 소개할 때 언급한 $E(Y)=b^{\prime }(\theta )=\mu$, $\mathrm{var}(Y)=\varphi b^{″}(\theta )=\varphi V(\mu )$을 증명할 것입니다. Exponential Family distribution은 다음과 같은 일반적인 형식으로 정의됩니다.

$$f(y;\theta ,\varphi )=\exp \{\frac{y\theta -b(\theta )}{a(\varphi )}+c(y,\varphi )\}$$ 때문에 log likelihood는 단순하게 아래와 같이 도출됩니다.

$$\ell (y;\theta ,\varphi )=\frac{y\theta -b(\theta )}{a(\varphi )}+c(y,\varphi )$$ 이제 이 log likelihood의 derivatives를 계산하면 다음과 같습니다:

• 첫 번째 도함수 (Score Function): $${\ell }^{\prime }(y;\theta ,\varphi )=\frac{y-b^{\prime }(\theta )}{a(\varphi )}$$

• 두 번째 도함수 (Observed Information): $${\ell }^{″}(y;\theta ,\varphi )=\frac{-b^{″}(\theta )}{a(\varphi )}$$

이 두 함수를 통해서 위에서 유도한 두 공식을 활용하면 다음 두 수식을 얻을 수 있습니다.

$$E\left\{\frac{Y-b^{\prime }(\theta )}{a(\varphi )}\right\}=0$$

$$E\left[{\left(\frac{Y-b^{\prime }(\theta )}{a(\varphi )}\right)}^2\right]=\frac{b^{″}(\theta )}{a(\varphi )}$$ 이때, 식을 잘 보면 첫 번째 식은 결국 $$E[Y-b^{\prime }(\theta )]=E[Y]-b^{\prime }(\theta )=0$$ 이 되어 $E\{Y\}=b^{\prime }(\theta )$을 얻을 수 있고, 두 번째 식에서 $$E\left[{\left(\frac{Y-b^{\prime }(\theta )}{a(\varphi )}\right)}^2\right]=\frac{E[(Y-b^{\prime }(\theta ){)}^2]}{E[a(\varphi {)}^2]}=\frac{\mathrm{Var}(Y)}{a(\varphi {)}^2}=\frac{b^{″}(\theta )}{a(\varphi )}$$ 이므로, $\text{Var}(Y)=b^{″}(\theta )a(\varphi )$임을 보일 수 있습니다.

증명한 수식을 다시 한 번 확인하면,

$$E(Y)=b^{\prime }(\theta )=\mu$$

$$Var(Y)=a(\varphi )V({\mu }_i)=a(\varphi )b^{″}(\theta )$$

2.3. GLMs’ parameter 추정식 유도

---

이제 필요한 식이 준비되었으니, 위에서 계속 다루고 있는 log likelihood을 이용해서 MLE estimation으로 GLMs’ parameter을 추정하는 과정을 살펴볼 것입니다. 이때, 추정 과정은 계속 언급한대로 Exponential Family distribution을 따르는 종속변수에 대해서 log likelihood의 model parameter에 대해 미분한 식이 (parameter가 벡터이므로, 좀더 엄밀하게 정의해야 하지만, 의미는 같으니 이렇게 얘기하겠습니다.) 0이 되게 하는 parameter를 찾음으로써 수행되며, 이 때의 함수 (log likelihood의 1차 도함수)를 앞으로는 score function이라고 부르겠습니다.

우리는 MLE estimation을 통해 여러 Exponential Family distributions에 대해 통일된 estimation algorithm으로 parameter를 추정할 수 있습니다. (이러한 분포 가정 마저 없다면, 3장 GEE에서 보겠지만 분포에 대한 직접적 가정없이 cumulative generating function 등 몇 함수 만으로 Likelihood를 고려하는 Quasi-likelihood Estimation의 개념으로 이어집니다.)

주어진 data가 $(y_1,...,y_n)$일 때, 위에서부터 계속 사용해왔던 log-likelihood function은 다음과 같습니다.

$$l=\sum _{i=1}^n(\frac{y_i{\theta }_i-b({\theta }_i)}{{\varphi }_i}+c(y_i,{\varphi }_i))$$

지금까지 우리는 $\theta$로 log likelihood를 다뤘지만, 추정해야 하는 parameter는 $\mathit{\beta }$입니다. $\mathit{\beta }$는 벡터이기 때문에 이 중 하나의 파라미터 ${\beta }_j$에 대해 먼저 log likelihood를 미분해보면 아래와 같은 식이 나옵니다. ($\text{Var}(y)={\varphi }_{i}V({\mu }_i)$, $\frac{\partial \mu }{\partial \eta }=\frac{1}{g^{\prime }({\mu }_i)}$)임은 위에서 보았습니다. g는 link function이었습니다.)

$$\frac{\partial l}{\partial {\beta }_j}=s({\beta }_j)$$

$$=\left(\frac{\partial l}{\partial \theta }\right)\left(\frac{\partial \theta }{\partial \mu }\right)\left(\frac{\partial \mu }{\partial \eta }\right)\left(\frac{\partial \eta }{\partial {\beta }_j}\right)$$

$$=\frac{y_i-{\mu }_i}{\text{Var}(y_i)}\left(\frac{\partial \mu }{\partial \eta }\right)x_{ij},\text{or}$$

$$=\sum _{i=1}^n\frac{y_i-{\mu }_i}{{\varphi }_{i}V({\mu }_i)}\times \frac{x_{ij}}{g^{\prime }({\mu }_i)}=0$$

식이 혼란스러울 수 있는데, 이는 단지 chain rule을 이용해서 $\beta$에 대한 $\ell$의 기울기를 구하는 과정이며, 위에서 GLM을 구성하는 과정을 차근차근 복기하면 각각의 변화율은 다음과 같이 구할 수 있기 때문에 최종 식을 얻을 수 있었음을 알 수 있습니다.

$$\frac{\partial l}{\partial \theta }=\frac{y-b^{\prime }(\theta )}{a(\varphi )}=\frac{y-\mu }{a(\varphi )}$$

$$\frac{\partial \theta }{\partial \mu }=\frac{1}{b^{″}(\theta )}=\frac{1}{V(\mu )}=\frac{a(\varphi )}{\text{Var}(y)}$$

$$\frac{\partial \eta }{\partial {\beta }_j}=x_{ij}$$

이제 이 score function의 음의 미분(또는 분산)의 기댓값을 전개하면 다음과 같습니다.

$$-E\left(\frac{{\partial }^{2}l}{\partial {\beta }_j\partial {\beta }_k}\right)=E\left[\left(\frac{\partial l}{\partial {\beta }_j}\right)\left(\frac{\partial l}{\partial {\beta }_k}\right)\right]$$

$$=E\left[{\left(\frac{y-\mu }{\text{Var}(y)}\right)}^2{\left(\frac{\partial \mu }{\partial \eta }\right)}^2{x}_{ij}{x}_{ik}\right]$$

$$=E\left[\frac{\text{Var}(y)}{\text{Var}(y{)}^2}{\left(\frac{\partial \mu }{\partial \eta }\right)}^2{x}_{ij}{x}_{ik}\right]$$

$$=\frac{1}{\text{Var}(y)}{\left(\frac{\partial \mu }{\partial \eta }\right)}^2{x}_{ij}{x}_{ik}$$

위 term은 Fisher Information matrix라고도 부르며, 이 term이 분산의 기댓값인 이유를 생각해보면, 이전에 구한 derivatives of log likelihood의 성질들에 의해 $\ell$의 negative 2차 도함수의 기댓값은 1차 도함수의 기댓값의 square와 같고, 이 1차 도함수의 기댓값이 0이므로 이는 분산과 같습니다. 정리하자면, score function의 미분식이 score function의 분산과 기댓값이 같으므로, 미분을 직접하는 대신 분산으로 근사 후 식을 전개한 것이며, 이 때 근사한 이 Matrix를 Fisher Information matrix라고 부릅니다.

이들을 한 번에 벡터와 행렬 연산으로 표현하면 다음과 같습니다:

$$\frac{\partial l}{\partial \beta }=X^{\mathrm{\top }}A(y-\mu )$$

$$E\left(\frac{{\partial }^{2}l}{\partial \beta \partial {\beta }^T}\right)=-X^{T}WX$$

$$where,W\text{is diagonal matrix comes from }\frac{1}{\text{Var}(y)}{\left(\frac{\partial \mu }{\partial \eta }\right)}^2,$$

$$\text{and }A\text{is diagonal matrix comes from }\frac{1}{\text{Var}(y)}\left(\frac{\partial \mu }{\partial \eta }\right).$$

이제 우리는 score function와 그 미분, 또는 log likelihood의 1차 도함수와 (negative) 2차 도함수의 추정식을 얻었습니다. 사실, $X^{\mathrm{\top }}A(y-\mu )$가 0이 되도록 하는 parameter $\beta$만 찾으면 parameter 추정이 끝나지만, 이는 closed-form solution이 존재하지 않기 때문에 그 미분(2차 도함수)식을 이용해 근사적으로 구할 수 있는 알고리즘을 최종적으로 다룰 것입니다. (물론 후에 보겠지만 이 Fisher Information matrix는 분산과도 관련이 있습니다.) 즉, 우리는 추정식을 유도하는 것은 완성했지만, 실제로 알고리즘을 설계하여 어떻게 이를 추정할 지에 대해서는 모르는 상태이고, 때문에 최종적으로 GLM을 제안한 학자가 소개하였으며 대부분의 패키지에서 이 GLM을 estimate하기 위해 사용하고 있는 method인 IRLS(Iteratively Reweighted Least Squares) Algorithm을 살펴볼 것입니다. (negative) 2차 도함수의 추정식(or Fisher Matrix) $X^{T}WX$을 유도한 이유는, 이 알고리즘에서 필요로 하기 때문이고, 위 식에서 $W,A$는 식이 복잡해보이지만, 그저 observations(data) 하나 당 GLM 모델의 구성요소를 통해 determinant하게 미리 계산되어 대각성분으로 각각 들어가는 term임을 명심하시면 좋을 것 같습니다. (이전 설명에서, 종속변수의 분포로 Exponential Family 중 특정 분포가 정해지면, 이에 따라 Link function(Canonical link), Variance function가 정해져 모델이 특정된다고 설명드린 적이 있고, 위 $W,A$모두 이 두 함수로 이루어진 식이기 때문에 관측치마다 각각 넣으면 determinant하게 하나의 값이 나오는 식인 것입니다.)

2.4. GLMs’ parameter 추정 (IRLS)

---

위에서 언급하였듯, GLM의 MLE estimation은 비선형 최적화 문제이기 때문에 공통된 framework에서 사용할 수 있는 closed-form solution이 존재하지 않으며, 대신 여러 최적화 방법을 사용할 수 있습니다. 뉴턴-랩슨 방법(Newton-Raphson Method)은 2차 도함수(Hessian Matrix)를 사용하여 score function을 수렴시키지만, Hessian Matrix $H=\frac{{\partial }^{2}l}{\partial \beta \partial {\beta }^T}$를 직접 구해야 하고, 이는 계산이 복잡하여 computation cost가 큽니다. Fisher Scoring은 Newton Method에서 Hessian Matrix 대신 이를 근사하는 Fisher Information Matrix를 사용합니다. 이는 이전에 Derivatives of log likelihood 의 성질이랑 추정식 유도에서 모두 보았던대로 2차 도함수가 1차 도함수의 분산(또는 제곱)와 기댓값이 같다는 수학적 성질을 토대로 $$E\left(\frac{{\partial }^{2}l}{\partial \beta \partial {\beta }^T}\right)=E\left[\left(\frac{\partial l}{\partial {\beta }_j}\right)\left(\frac{\partial l}{\partial {\beta }_k}\right)\right]=X^{T}WX$$ 가 만족함을 확인하였기 때문에, MLE추정에서 Newton-Raphson Method보다 computation cost가 합리적인 method라고 생각해 볼 수 있습니다. 이외에도 경사 하강법(Gradient Descent) 알고리즘은 1차 도함수(Gradient)만 사용하여 특정 값만큼 조금씩 점진적으로 parameter를 움직여 최적점을 찾는 방법입니다. 모델이 매우 복잡해서 2차 도함수를 계산하기 힘든 딥러닝에서는 많은 경우 이를 발전시킨 여러 methods로 iterativaly하게 parameter를 추정합니다. (이렇게 iterative하게 model’s parameter를 움직이면서 추정하는 과정이 AI에서 얘기하는 learning입니다.)

위에서 얘기한 IRLS(Iteratively Reweighted Least Squares) Algorithm는 이 Fisher Scoring의 알고리즘적 변형으로, Fisher Scoring의 식을 가중 최소제곱(Weighted Least Squares) 문제로 치환하여, 이 문제에서 사용하는 IRLS 알고리즘으로 GLM의 parameter 해를 구하는 method입니다. 우선 Newton-Raphson Method, Fisher Scoring에 대한 이야기를 간단하게 하고, 자세하게 어떻게 IRLS가 GLM의 parameter를 추정하는지 보겠습니다.

Newton-Raphson Method & Fisher Scoring

---

GLM (Generalized Linear Model)의 파라미터 추정을 위한 최적화 과정은 우선 log likelihood function의 최대화를 목적으로 합니다. 이때, Newton-Raphson Method와 그 변형인 Fisher Scoring은 모두 이를 위한 알고리즘입니다.

Newton-Raphson 방법은 다음과 같은 일반적인 업데이트 식을 갖습니다:

$${\mathit{\beta }}^{(t+1)}={\mathit{\beta }}^{(t)}-{\left[\frac{{\partial }^2\ell }{\partial \mathit{\beta }\partial {\mathit{\beta }}^{\mathrm{\top }}}\right]}^{-1}\frac{\partial \ell }{\partial \mathit{\beta }},$$ 이후의 method들은 모두 이 method에서 기인하므로, 위 식의 의미를 이해하는 것은 매우 중요합니다. 위 식이 어떻게 안정적으로 $\mathit{\beta }$를 수렴시킬 수 있는지 2차 테일러 전개를 통해 수식적으로 좀 더 명확히 볼 수 있지만, 여기서는 좀더 직관적으로 가볍게 이해해보겠습니다.

Newton-Raphson 방법은 어떠한 함수 $f(x)$에 대해서 함수의 해, 즉 $f(x)=0$을 만족하는 $x$를 찾기 위한 반복적인 근사 방법입니다.(informal한 증명이므로 1-dimension case로 보겠습니다.) 이 방법은 현재 점에서 함수의 접선을 그려 $x$축과 만나는 점을 다음 근사해로 사용합니다. 예를 들어, 초기값 $x_0$에서 접선을 그리면 그 접선의 방정식은 $y=f^{\prime }(x_0)(x-x_0)+f(x_0)$입니다. 이 접선이 $x$ 축과 만나는 점은 $y=0$일 때이므로, 이를 대입하여 접선이 0이 되는 점을 구해보면 $x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$가 됩니다. 그러나 접선은 함수의 선형 approximation이므로 이 접선이 0이 되는 점 $x_1$이 실제 함수에서도 곧바로 0이 되지는 않습니다. 따라서 이 과정을 반복하여 특정 단계에서 $f(x_n)$이 0에 충분히 가까워지면, $x_n$을 근사해로 채택합니다.

이 방법이 작동하는 이유는 접선의 기울기 $f^{\prime }(x_n)$이 함수의 곡률을 반영하기 때문입니다. 곡률이 클수록(기울기가 가파를수록) 업데이트의 크기가 작아지고, 곡률이 작을수록 업데이트의 크기가 커집니다. 또한, Newton-Raphson 방법은 2차 수렴(Quadratic Convergence) 속도를 가집니다. 이는 오차가 반복마다 제곱으로 줄어들기 때문에 매우 빠르게 해에 수렴한다는 의미입니다.

이정도로 간단하게 Newton-Raphson method를 이해할 수 있고, 다시 돌아와서 위 식에서는 multi-dimention 상황에서 해를 찾고 싶은 함수가 score function, 즉 $\frac{\partial \ell }{\partial \beta }$인 경우이기 때문에 위처럼 식이 구성되었다는 것을 알 수 있습니다.(f의 미분이 분모로 들어간 term은 행렬에서 역행렬 -1과 같은 의미라고 보시면 됩니다.) 첨언하자면, 이 경우 Newton-Raphson method는 두 가지 중요 조건이 붙는데, 언급만 하자면 Hessian matrix가 Positive-definite (볼록) 해야 하며, 초기값이 최적점에 충분히 가까워야 합니다.

Fisher Scoring은 Hessian 행렬 대신 Fisher Information 행렬 $\mathcal{I}(\mathit{\beta })$를 사용하여 다음과 같이 업데이트합니다:

$${\mathit{\beta }}^{(t+1)}={\mathit{\beta }}^{(t)}+{(\mathcal{I}({\mathit{\beta }}^{(t)}))}^{-1}\mathbf{S}({\mathit{\beta }}^{(t)}),$$

여기서 $$\mathbf{S}(\mathit{\beta })=\frac{\partial \ell }{\partial \mathit{\beta }}$$ 는 Score function으로 구한 (Fisher) score 벡터이고, $$\mathcal{I}(\mathit{\beta })=E\left[\left(\frac{\partial l}{\partial {\beta }_j}\right)\left(\frac{\partial l}{\partial {\beta }_k}\right)\right]=E[-\frac{{\partial }^2\ell }{\partial \mathit{\beta }\partial {\mathit{\beta }}^{\mathrm{\top }}}]$$ 는 Fisher information matrix입니다.

즉, Fisher Scoring 방법은 뉴턴-랩슨 방법(Newton-Raphson Method)에서 Hessian Matrix $H=\frac{{\partial }^{2}l}{\partial \beta \partial {\beta }^T}$ 를 사용하는 대신, Fisher information matrix를 사용해서 업데이트를 수행함으로써 parameter를 estimation하는 매커니즘입니다. IRLS는 GLM에서 이 Fisher Scoring와 거의 일치하다 봐도 무방하며, 단순히 위에서 추정한 $\mathbf{S}(\mathit{\beta })$와 $\mathcal{I}(\mathit{\beta })$를 Fisher Scoring 공식에 넣으면 weighted least squares problme(가중 최소제곱 문제)와 완전히 유사해지기 때문에, 이 문제를 해결하는 방식으로 parameter를 추정한다는 의미라고 생각하시면 될 것 같습니다. 좀 더 수식과 같이 자세하게 설명드리겠습니다.

IRLS (Iteratively Reweighted Least Squares) Algorithm

---

앞서, 2.3.에서 log likelihood의 gradient(1차 도함수, Score function)와 expected Hessian(Fisher Information matrix)가 각각

$$\frac{\partial l}{\partial \beta }=X^{T}A(y-\mu )$$

$$E\left(\frac{{\partial }^{2}l}{\partial \beta \partial {\beta }^T}\right)=X^{T}WX$$

$$where,A=\frac{1}{\text{Var}(y)}\left(\frac{\partial \mu }{\partial \eta }\right)andW=\frac{1}{\text{Var}(y)}{\left(\frac{\partial \mu }{\partial \eta }\right)}^2$$

임을 보았습니다. 이 결과들을 Fisher Scoring method의 식에 대입하면,

$${\beta }^{(t+1)}={\beta }^{(t)}+{\left(X^{T}WX\right)}^{-1}{X}^{T}A(y-\mu )$$

이고, 이 equation은 가중 최소제곱 문제(Weighted Least Squares, WLS)의 parameter 추정식과 일치하기 때문에 비선형 GLM의 parameter 추정을 WLS problem으로 치환할 수 있음을 알 수 있습니다.

위의 업데이트 식은 아래 working response $z$를 정의하면

$$z=X{\beta }^{(t)}+W^{-1}A(y-\mu ),$$

아래와 같이 표현할 수 있습니다.

$${\beta }^{(t+1)}={\beta }^{(t)}+{\left(X^{T}WX\right)}^{-1}{X}^{T}A(y-\mu )$$ $$={\left(X^{T}WX\right)}^{-1}\left(X^{T}WX\right){\beta }^{(t)}+{\left(X^{T}WX\right)}^{-1}{X}^{T}W{W}^{-1}A(y-\mu )$$ $$={\left(X^{T}WX\right)}^{-1}{X}^{T}W(X{\beta }^{(t)}+W^{-1}A(y-\mu ))$$

$$={\left(X^{T}WX\right)}^{-1}{X}^{T}Wz$$

또 강조하지만, 이는 가중치 행렬 $W$에 따라 각 관측치의 기여도를 달리하는 선형 회귀 문제(가중 최소제곱 문제)의 정규방정식과 동일합니다:

$$(X^{T}WX)\beta =X^{T}Wz.$$

따라서 해당 정규방정식의 $\beta$를 가중 최소제곱 문제 방식으로 풀어냄으로써 추정치 ${\beta }^{(t+1)}$를 구할 수 있으며, 이는 현재 단계의 추정치 ${\beta }^{(t)}$에서의 예측값과 오차를 반영한 새로운 업데이트가 됩니다.

IRLS 구체적 절차

---

정리하자면, IRLS(Iteratively Reweighted Least Squares) 알고리즘은 위의 아이디어를 바탕으로 GLM의 최대우도추정 문제를 반복적으로 가중 최소제곱 문제로 전환하여 해결합니다. 구체적인 단계는 다음과 같습니다:

(1) 초기화: 초기 파라미터 ${\beta }^{(0)}$를 설정합니다.

(2) 현재 단계 계산:

(2.1) 예측값 계산: 현재 추정치 ${\beta }^{(t)}$을 이용하여 선형 예측치 ${\eta }^{(t)}=X{\beta }^{(t)}$를 구하고, link 함수의 역함수를 통해 ${\mu }^{(t)}=g^{-1}({\eta }^{(t)})$를 계산합니다.

(2.2) 가중치 및 보조 행렬 계산: 정의한 식에 따라

$$A^{(t)}=\frac{1}{\text{Var}(y)}{\left(\frac{\partial \mu }{\partial \eta }\right)}^{(t)}$$

$$W^{(t)}=\frac{1}{\text{Var}(y)}{\left(\frac{\partial \mu }{\partial \eta }\right)}^{(t)2}$$

을 계산합니다.

(2.3) Working Response 구성:

$$z^{(t)}=X{\beta }^{(t)}+{\left(W^{(t)}\right)}^{-1}{A}^{(t)}(y-{\mu }^{(t)}).$$

(3) 가중 최소제곱 문제 해결:

위의 working response와 가중치 행렬을 사용하여 정규방정식

$$(X^T{W}^{(t)}X){\beta }^{(t+1)}=X^T{W}^{(t)}{z}^{(t)}$$ 을 풀어 새로운 추정치 ${\beta }^{(t+1)}$를 구합니다.

(4) 수렴 판단 및 반복:

$||{\beta }^{(t+1)}-{\beta }^{(t)}||$(L1 norm, 쉽게는 절댓값)가 미리 설정한 임계값 이하가 될 때까지 2번과 3번의 단계를 반복하고, 수렴이 되었다면 이 ${\beta }^{(t+1)}$ 값이 GLM의 parameter에 대한 IRLS의 최종 Estimation 결과입니다. 결국, Fisher Scoring에서 대입한 수식이 결국 가중 최소제곱 문제로 귀착됨을 통해, IRLS 알고리즘은 각 반복마다 선형 회귀 문제와 유사한 방식으로 parameter를 업데이트합니다.

이 IRLS의 소프트웨어 구현에 대한 첨언을 하자면, 보통 IRLS에서는 각 단계마다 위 3단계와 같이 아래 정규방정식을 풉니다:

$$(X^{T}WX)\beta =X^{T}Wz.$$

이때 직접 $(X^{T}WX{)}^{-1}$를 구해 업데이트하는 방법은 계산적으로 불안정할 수 있습니다. 특히, 데이터의 규모가 크거나 $X$ 행렬이 ill-conditioned(조건수가 열악한)인 경우에는 직접 역행렬을 계산하는 과정에서 수치적인 문제가 발생할 위험이 큽니다. 따라서, 구현의 영역이기 때문에 더이상 나열하지는 않겠지만 실제는 역행렬을 direct하게 구하는 대신, QR 분해나 Cholesky 분해 같은 선형대수 기법을 활용하여 안정적으로 선형 시스템을 풀 수 있습니다.

2.5. GLMs’ parameter Variance

---

앞서 IRLS(Iteratively Reweighted Least Squares) 알고리즘으로 GLM의 파라미터를 추정하는 과정을 살펴보았습니다. 이때 우리는 MLE(최대우도추정)을 IRLS(반복적으로 가중 최소제곱)문제로 전환하는 과정을 거쳤는데, 최종적으로 구해지는 추정치 $\hat{\mathit{\beta }}$의 분산에 대한 추정까지 마쳐야 유의성 검정 등의 분석을 수행할 수 있을 것입니다. GLM에서 최대우도추정(MLE)을 사용해 얻은 $\hat{\mathit{\beta }}$는, 이론적으로 위에서 언급한 Fisher 정보 행렬(Fisher information matrix)의 역행렬로써 구할 수 있습니다:

$$\widehat{\mathrm{Var}}(\hat{\mathit{\beta }})=(\mathbf{I}(\hat{\mathit{\beta }}){)}^{-1},$$

여기서 $\mathbf{I}(\hat{\mathit{\beta }})$는 $\hat{\mathit{\beta }}$에서의 (observed 혹은 expected) Fisher 정보 행렬입니다.(3장에서 이 이유에 대해 살펴볼 것입니다.) 모델의 parameter를 추정하는 과정에서, 각 관측치 $i$에 대해 $\mathrm{Var}(y_i)$와 $\frac{\partial {\mu }_i}{\partial {\eta }_i}$가 포함된 특정 가중치 $W^{(t)}$가 등장하였었고, 반복(step)마다 업데이트되는 정규방정식을 풀어감으로써 추정치 ${\beta }^{(t+1)}$를 얻었습니다. 이후 최종 수렴 시점($t\to \mathrm{\infty }$)에서, 우리는 $\hat{\mathit{\beta }}={\beta }^{(\mathrm{\infty })}$에 도달하게 되고, 이 시점에서 계산된 Hessian(또는 Fisher information) matrix ${\mathbf{X}}^{\mathrm{\top }}\hat{\mathbf{W}}\mathbf{X}$에 대해서, 이 행렬의 역수가 분산이 되는 것입니다.

즉, 실제 계산 시에는 아래와 같은 모양이 됩니다.

$$\widehat{\mathrm{Var}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\hat{\mathbf{W}}\mathbf{X}{)}^{-1},$$

왜 최종적으로 추정된 모델의 분산이 수렴된 time step에서의 Fisher information matrix로 추정할 수 있는지에 대해서는 3장의 M-estimation에서 증명할 것입니다. 여기서 미리 강조할 것은, 위에서 얻는 분산의 추정값은 GLM에서의 기본 가정들이었던 독립성, 분산 함수의 형태 등이 성립한다는 조건 위에서 도출된 것이므로, 이러한 기본 가정이 어느 정도 엄격히 맞아떨어지는 상황(정확한 포아송 분포를 따르는 count data, 명시적으로 독립적인 개별 관측치 등)이라면 괜찮지만, 현실의 data에서는 이질분산, 클러스터 내 상관, 과산포(overdispersion) 등으로 인해 이 기본 가정들이 깨질 수 있습니다. 다행히도, GLM에서도 모델이 consist할 때(GLM의 위로부터 추정된 parameter 자체는 consist합니다.) HC(Heteroskedasticity-Consistent) se와, 아래에서 clustered data에서 고려할 수 있는 버전인 Cluster-robust standard errors를 사용하여 더욱 robust하게 분산을 추정할 수 있습니다. 때문에 아래에서는 Cluster-robust se를 OLS 버전으로 소개드린 후, GLM에서 사용하기 위해 $\hat{\mathbf{W}}$ 행렬 $\hat{\mathbf{A}}$ 행렬 등의 구조가 어떻게 수식적으로 첨가되어 LM(Linear Model)에서의 식과는 살짝 다른 모습을 취하게 되는지 보고 마치겠습니다.

3. Cluster-Robust Standard Errors

3.1. Clustered Data 정의

---

Clustered data란 데이터 내에서 동일 그룹에 속하는 관측치들이 상관관계를 가지는 경우를 의미합니다. 예를 들어, 한 환자의 여러 진료 기록이 서로 상관되어 있을 수 있습니다. 이 때, cluster간에는 상관관계가 없고 cluster 내의 데이터들은 상관관계가 있다는 가정하에 Cluster-robust standard errors나 3장의 GEE, GLMM model이 개발되었습니다. 의료 분석 상황의 예시로는 대표적으로 데이터의 각 환자 당 여러 시간 또는 주기에 걸쳐 측정한 데이터, 여러 학교나 병원과 같은 단체에서 얻은 데이터들을 한 번에 고려하는 경우가 있을 것입니다. 또한, cluster간에 상관관계가 있거나 cluster 안에 cluster가 있는 hierarchical의 경우도 있지만, 이에 대한 공식들은 위에서 고려하는 1차적인 상황을 이해하면 쉽게 이해할 수 있으며, 의료 분석에서 고려하는 피험자 내 관측치 간 상관관계, 병원 내 관측치 간 상관관계 등을 고려해야 하는 상황은 이번 블로그에서 이야기 할 1차적인 clustered 상황임을 알아두시면 좋을 것 같습니다. 이 observations간의 상관관계에 대한 이야기와, 이때 사용해야 하는 Regression Models에 대한 설명은 3장에서 GEE, GLMM과 함께 더욱 자세하게 다뤄볼 예정입니다.

확실한 것은, 비선형 분포를 추정할 수 있는 GLM이나, OLS에서 안정적인 parameter 분산 추정 method였던 HC(Heteroskedasticity-Consistent) 표준오차는 관측치 간의 독립성을 가정하였었고, 이는 위와 같은 data를 다룰 때에는 깨져야 하는 가정이라는 것입니다. 이제 설명드릴 Cluster-robust standard errors는 HC se와 형태가 매우 비슷하며, 같은 철학으로 clustered data에서 robust한 모델 분산 추정 method입니다. 이때 기억하셔야 할 부분은, 1장에서는 HC se의 안정성을 Linear Regression에 대해서 고려하였고 R의 sandwich 패키지를 통해 구현할 수 있음을 보았는데, 이번 장에서 다루고 있는 GLM에서도 이 HC se, Cluster-robust se를 모두 사용할 수 있다는 것입니다. 이때 식이 LM과 GLM에서 살짝 다른데, 우선 Linear Model에서의 Cluster-robust se에 대해서 설명드리고, GLM에서는 무엇이 다른지 보겠습니다.

실제 소프트웨어의 구현에 대해서 첨언하자면, R의 sandwich 패키지나 대부분의 패캐지에서는 이 두 robust 분산 추정의 계산 및 검정을 LM, GLM 모두에 사용 가능하고, 이 패키지들은 들어오는 모델의 객체가 LM, GLM임을 분류한 뒤 각각에 맞는 살짝 변형된 식으로 추정한다고 생각하시면 될 것 같습니다.

3.2. Cluster-robust standard errors 정의 및 수학적 표현

---

Cluster-robust standard errors는 클러스터 내 상관관계를 고려하여 분산을 추정합니다. 이를 통해 클러스터 간 독립성은 유지하되, 클러스터 내 관측치 간 상관관계가 존재할 때도 일관된 추정치를 제공합니다. LM에서 Cluster-robust standard error를 구하는 식은 다음과 같습니다:

$$\widehat{\text{Var}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}\left(\sum _{g=1}^G{\mathbf{X}}_g^{\mathrm{\top }}{\hat{\mathbf{u}}}_g{\hat{\mathbf{u}}}_g^{\mathrm{\top }}{\mathbf{X}}_g\right)({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

위 식에서 $g$는 클러스터 인덱스, ${\hat{\mathbf{u}}}_g$는 클러스터 $g$의 잔차 벡터입니다.이 식은 1장에서의 HC0과 아주 유사하며, 가운데 meat항 (두 $({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$ 사이에 있는 항)만 달라졌음을 알 수 있습니다. 이는 실제로 HC0에서, 위에서 설명드린 가정인 1차적 clustered 구조(클러스터 간 독립성을 가정하지만, 클러스터 내 관측치들 간 상관관계는 허용)를 고려해서 $\mathrm{\Phi }$ 항만 바뀌었음을 짐작해볼 수 있습니다. 이제 이 Cluster-robust standard errors의 철학에 대해서 구체적으로 살펴보겠습니다.

3.3. Cluster-robust standard errors 수학적 표현

---

LM에서는 이전에 봤던대로 parameter의 분산을 유도하면 다음과 같습니다:

$$\widehat{\text{Var}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{\Phi }\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

여기서 $\mathrm{\Phi }$는 오차항의 공분산 행렬을 나타냈었습니다. HC0 (Heteroskedasticity-Consistent 0)에서는 모든 관측치가 서로 독립임을 가정하였기 때문에 (Heteroskedasticity를 고려하였지 dependent case를 고려하지는 않았었습니다.) 이에 따라 $\mathrm{\Phi }$는 대각행렬로 표현되며,

$${\mathrm{\Phi }}_{\text{HC0}}=\mathrm{diag}({\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_n^2)$$

결과적으로 분산 추정량은 개별 관측치에 대해아래와 같이 계산하였었습니다.

$${\widehat{\text{Var}}}_{\text{HC0}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}\left(\sum _{i=1}^n{\mathbf{x}}_i{\hat{u}}_i^2{\mathbf{x}}_i^{\mathrm{\top }}\right)({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$ 위 식은 1장의 HC0 식과 같은 식입니다. 표현이 어색하다고 느끼시는 분을 위해 이전에 사용한 식을 가져오면 $${\widehat{\text{Var}}}_{\text{HC0}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{diag}(e_i^2)\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

이며, 위에서부터 얘기하고 있는 $\hat{u}$는 이 error term $e$와 비슷한(잔차이기 때문에 사실 의미는 다릅니다) 의미입니다. HC0에서는 각 관측치만을 고려하기 때문에 ${\hat{u}}_i$는 $e_i$와 같고, 길이가 1인 벡터, 즉 scalar이기 때문에 제곱을 사용하였지만 위 cluster-robust 식에서 사용한 ${\hat{u}}_g$는 클러스터 g에 해당하는 모든 관측치를 한 줄로 나열한 임의의 길이의 벡터이기 때문에 제곱이 아니라 ${\hat{\mathbf{u}}}_g{\hat{\mathbf{u}}}_g^{\mathrm{\top }}$ term을 사용한 것입니다.

이에 대한 이해를 바탕으로 Cluster-robust se의 meat term을 생각해보면, Cluster-robust에서는 cluster간은 독립적이고, cluster안의 관측치들은 상관관계를 가질 수 있다고 가정하기 때문에 각 cluser에 대해서 ${\mathrm{\Phi }}_i$를 ${\hat{\mathbf{u}}}_g{\hat{\mathbf{u}}}_g^{\mathrm{\top }}$로 각각 구한 후, 전체 $\mathrm{\Phi }$는 아래와 같이 block diagonal 구조로 넣어준다고 이해할 수 있습니다. (block 행렬은 행렬을 특정한 block으로 나누었을 때 대각선 이외의 모든 행렬 블록이 영행렬인 행렬을 의미하며, cluster간의 독립을 block diagonal 구조로 고려하였다고 이해하면 됩니다.)

$${\mathrm{\Phi }}_{\text{cluster}}=\left(\begin{array}{cccc}{\mathrm{\Phi }}_1& 0& \cdots & 0\\ 0& {\mathrm{\Phi }}_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\mathrm{\Phi }}_G\end{array}\right)$$

즉, 여기서 각 ${\mathrm{\Phi }}_g=E[{\mathbf{u}}_g{\mathbf{u}}_g^{\mathrm{\top }}]$는 클러스터 $g$ 내의 오차의 공분산 행렬이고, 각 cluster에 대해 잔차${\hat{\mathbf{u}}}_g$를 사용하여

$${\widehat{\text{Var}}}_{\text{cluster}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}\left(\sum _{g=1}^G{\mathbf{X}}_g^{\mathrm{\top }}{\hat{\mathbf{u}}}_g{\hat{\mathbf{u}}}_g^{\mathrm{\top }}{\mathbf{X}}_g\right)({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$ 와 같이 추정합니다.

위 Cluster-robust의 meat항에 대한 이해는 비슷하게 3장에서도 필요하기 때문에 예시를 통해 좀더 직관적으로 보여드리겠습니다. 우선 이 중앙항은 다음과 같고,

$$\mathbf{B}=\sum _{g=1}^G{\mathbf{X}}_g^{\mathrm{\top }}{\hat{\mathbf{u}}}_g{\hat{\mathbf{u}}}_g^{\mathrm{\top }}{\mathbf{X}}_g$$

이는 행렬로 보면 위에서 보신 $\mathrm{\Phi }$항과 같이 block diagonal 형태를 갖습니다. 3개의 cluster가 있고, 각 cluster 내 관측치 수가 2, 3, 1개라고 가정하면 각각의 $\mathrm{\Phi }$는 다음과 같고,

$${\mathrm{\Phi }}_1={\mathbf{X}}_1^{\mathrm{\top }}{\hat{\mathbf{u}}}_1{\hat{\mathbf{u}}}_1^{\mathrm{\top }}{\mathbf{X}}_1=\left(\begin{array}{cc}{\sigma }_{11}^2& {\sigma }_{12}\\ {\sigma }_{12}& {\sigma }_{22}^2\end{array}\right)$$

$${\mathrm{\Phi }}_2={\mathbf{X}}_2^{\mathrm{\top }}{\hat{\mathbf{u}}}_2{\hat{\mathbf{u}}}_2^{\mathrm{\top }}{\mathbf{X}}_2=\left(\begin{array}{ccc}{\sigma }_{33}^2& {\sigma }_{34}& {\sigma }_{35}\\ {\sigma }_{34}& {\sigma }_{44}^2& {\sigma }_{45}\\ {\sigma }_{35}& {\sigma }_{45}& {\sigma }_{55}^2\end{array}\right)$$

$${\mathrm{\Phi }}_3={\mathbf{x}}_6{\hat{u}}_6^2{\mathbf{x}}_6^{\mathrm{\top }}$$

로 표현될 수 있으며, 결국 Cluster-robust의 중앙 term은

$$\mathbf{\Phi }=\left(\begin{array}{ccc}{\mathrm{\Phi }}_1& 0& 0\\ 0& {\mathrm{\Phi }}_2& 0\\ 0& 0& {\mathrm{\Phi }}_3\end{array}\right)$$ 가 될 것입니다.

정리하자면, HC0는 $\mathrm{\Phi }$가 대각행렬인 경우로, 개별 관측치의 Heteroskedasticity 만을 고려하며

$${\widehat{\text{Var}}}_{\text{HC0}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}\left(\sum _{i=1}^n{\mathbf{x}}_i{\hat{u}}_i^2{\mathbf{x}}_i^{\mathrm{\top }}\right)({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

Cluster-Robust 분산 추정량은 clusr별 $\mathrm{\Phi }$가 block diagonal 구조로, Heteroskedasticity와 cluster 내의 상관관계를 반영합니다.

$$E[\widehat{\text{Var}}(\hat{\mathit{\beta }})]=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{\Phi }\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

3.4. In GLMs..

---

이미 위(3장에서) Cluster-robust standard errors가 어떤 원리로부터 유도되며, OLS 환경(Linear Model)에서의 공식이 어떻게 생겼는지 살펴보았습니다. 또한 HC(Heteroskedasticity-Consistent) se 역시 기본 가정(등분산, 독립성 등)이 약화되었을 때도 일관된 추정을 제공하기 위해 Robust(샌드위치) 분산 추정량을 쓰게 된다는 것을 보았습니다. GLM에서도 LM과 같이 위 두 robust한 분산 추정치 식을 사용하여 Fisher information matrix의 역행렬로 분산을 추정하는 대신, 더욱 안정적으로 분산을 추정할 수 있습니다. GLM의 경우, 단순 OLS와 달리 $\hat{\mathbf{W}},\hat{\mathbf{A}}$ 등 추가적인 항이 존재하고, 이 행렬들이 실제 분산 추정 과정에 반영됩니다. 이 때문에 “bread”(양쪽에 곱해지는 행렬)와 “meat”(중간에 오는 분산·잔차 구조) 부분이 LM에서의 표기와는 형태가 조금 달라집니다. 즉, 원리는 동일하되, link & variance function으로 부터 비롯된 미분 항($\mathbf{A}$)과 가중치 항($\mathbf{W}$)이 반영되어야 한다는 점만 다릅니다. 1장에서 소개했던 HC0를 떠올리면, LM의 경우

$${\widehat{\mathrm{Var}}}_{\text{HC0}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}\left(\sum _{i=1}^n{\mathbf{x}}_i{\hat{u}}_i^2{\mathbf{x}}_i^{\mathrm{\top }}\right)({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1},$$

로 식이 구성되었습니다. GLM의 경우에는 $\hat{\mathbf{W}}$가 $\mathbf{X}$와 상호작용하여 분산 추정에 들어가므로, 실제로는 다음과 같은 형태를 가집니다. (식은 패키지나 저자별 표기 차이에 따라 다소 달라질 수 있습니다). $${\widehat{\mathrm{Var}}}_{\text{HC0}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\hat{\mathbf{W}}\mathbf{X}{)}^{-1}\left(\sum _{i=1}^n{\mathbf{x}}_i({\hat{u}}_i^2){\mathbf{x}}_i^{\mathrm{\top }}\right)({\mathbf{X}}^{\mathrm{\top }}\hat{\mathbf{W}}\mathbf{X}{)}^{-1}.$$

여기서 ${\hat{u}}_i$는 단순 잔차가 아니라, 펄슨(pearson) 잔차 등 비선형적인 GLM 설정에 맞춰 적절히 조정된 형태일 수 있습니다. 구현별로 이탈도(deviance) 잔차를 사용할 수도 있고, 핵심은 “관측치별 잔차의 크기”를 통해 이질분산성을 추정하는 것입니다. 여기서는 앞뒤의 $({\mathbf{X}}^{\mathrm{\top }}\hat{\mathbf{W}}\mathbf{X}{)}^{-1}$이 bread(빵)이고, 가운데 잔차 ${\hat{u}}_i{\hat{u}}_i^{\mathrm{\top }}$가 meat(고기) 역할을 한다고 보면 됩니다. 철학적으로 해석하면, GLM에서도 LM에서와 동일하게 HC se는 “각 관측치별 오차분산”이 서로 다르더라도 일관된 추정을 제공하기 위하는 목적이며, 식은 (1) 잔차(오차항) 부분은 그대로 meat로 넣고, (2) 정보를 제공하는 bread에는 $\mathbf{X}$에 가중치의 의미를 가진 $\hat{\mathbf{W}}$ 항을 추가하여 구성한 위 형태로 구성됩니다.

클러스터링이 있는 데이터에 대하여, LM과 마찬가지로 GLM에서도 Cluster-robust se가 적용될 수 있습니다. 이미 섹션 3.2~3.3에서 보았듯, 클러스터 간에는 독립이지만 클러스터 내 관측치들 간에는 상관관계가 존재할 수 있으므로, $\mathrm{\Phi }$ 행렬(오차의 공분산 구조)을 block diagonal 형태로 가정하고, 이를 샌드위치 가운데(meat)에 반영합니다.

LM에서의 일반적 식은 다음과 같았습니다.

$${\widehat{\mathrm{Var}}}_{\text{cluster}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}\left(\sum _{g=1}^G{\mathbf{X}}_g^{\mathrm{\top }}{\hat{\mathbf{u}}}_g{\hat{\mathbf{u}}}_g^{\mathrm{\top }}{\mathbf{X}}_g\right)({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}.$$

GLM에서는 동일한 철학으로, 단순히 $\mathbf{X}$ 대신 가중치를 고려해 ${\hat{\mathbf{W}}}^{1/2}\mathbf{X}$와 같은 형태(혹은 관련 도함수 항)가 곱해지게 됩니다. 즉,

$${\widehat{\mathrm{Var}}}_{\text{cluster}}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\hat{\mathbf{W}}\mathbf{X}{)}^{-1}(\sum _{g=1}^G({\mathbf{X}}_g^{\mathrm{\top }}{\hat{\mathbf{W}}}_g^{1/2}){\hat{\mathbf{u}}}_g{\hat{\mathbf{u}}}_g^{\mathrm{\top }}({\hat{\mathbf{W}}}_g^{1/2}{\mathbf{X}}_g))({\mathbf{X}}^{\mathrm{\top }}\hat{\mathbf{W}}\mathbf{X}{)}^{-1},$$

와 같은 꼴이 됩니다(마찬가지로 패키지마다 표기 방식이나 구현 세부가 약간씩 다를 수 있습니다).각 항들 또한 한 번 더 설명하자면, ${\hat{\mathbf{u}}}_g$는 클러스터 $g$ 내 잔차 벡터(pearson 또는 deviance 잔차 등).${\mathbf{X}}_g$는 클러스터 $g$에 해당하는 행만 추출한 $\mathbf{X}$의 서브 행렬. $\hat{\mathbf{W}}$는 클러스터 ($g$)에 해당하는 마찬가지로 가중치 서브 행렬이고, 이때 1/2승을 한다는 의미는 이가 diagonal matrix이므로 이 경우에는 단순히 diagonal 성분들 각각을 루트 씌운 값입니다.

4. R 예시: GLM, Cluster-robust SE

아래 R 코드를 복사하여 로컬 환경에서 돌려보세요. GLM모델의 분산과 cluster-robust 분산을 비교하시면서 해석하면 됩니다.

## 필요한 패키지 설치 (필요시)
## install.packages("sandwich")
## install.packages("lmtest")
## install.packages("nlme")
#
## 데이터 불러오기
#library(nlme)
#data(Orthodont)  # 치아 성장 데이터 (클러스터: Subject)
#Orthodont$binary <- ifelse(Orthodont$distance > 25, 1, 0)  # 이항 변환
#
## 기본 GLM (로지스틱 회귀)
#glm_fit <- glm(binary ~ age + Sex, 
#               data = Orthodont, 
#               family = binomial)
#summary(glm_fit) 
## 클러스터-로버스트 표준오차 (Subject 기준)
#library(sandwich)
#library(lmtest)
#cluster_se <- vcovCL(glm_fit, cluster = ~ Subject)
#coeftest(glm_fit, vcov = cluster_se)  # 결과 출력

마무리하며

이번 2장에서는 1장에서 다룬 Linear Model을 outcome of single yes/no, outcome of single K-way, count 등 non-normal한 종속변수에서도 분석할 수 있도록 확장한 Generalized linear model의 기본 개념과, 실제로 패키지에서 이 GLM의 parameter를 estimate할 때 사용하는 대표적인 알고리즘인 IRLS(Fisher scoring)을 수학적으로 상당히 깊게 살펴보았습니다. 이후, HC standard errors의 clustered data 버전인 Cluster-robust standard error를 보고, GLMs에서도 이 둘을 사용할 수 있다는 것을 밝힌 뒤 그 변형된 수식을 보았습니다. 다음 3장에서는 아직 깨지 못한 가정이었던 오차항의 독립, 즉 data(observations)간의 correlation이 존재하는 경우 자체를 모델에 반영하기 위해 개발된 모델들인 GEE, GLMM에 대하여 어느 정도 살펴보고 (GLMM의 내용은 너무 길어지기 때문에 얕게 다룰 것입니다.), 모델의 분산을 robust하게 추정하기 위한 가장 general한 형태의 Sandwich estimator를 M-estimation의 개념과 함께 공부할 것입니다.