Exploring Regression Models for Regression Analysis (3): GEE, GLMM, M-statistics, Robust (sandwich) estimation

들어가며

3장에서는 2장에서 다룬 Generalized Linear Model (GLM)에서 더 나아가, 데이터 내에 군집(Clustered) 구조가 존재하거나, 반복측정(Repeated measures) 데이터로 인해 독립성 가정이 깨지는 경우를 다루는 방법론인 GEE (Generalized Estimating Equation)와 GLMM (Generalized Linear Mixed Model)을 다루며, 그 전에 중요한 추정 방법론 중 하나인 M-estimation과 Robust(Sandwich) estimation에 대해서 다루겠습니다.

1. M-estimation

1.1. M-estimation 정의

---

통계 분석에서 통계 모델이 비모수(non-parametric)가 아니라 모수(parametric)인 경우, 우리는 model의 모수, 즉 parameter ($\mathit{\theta }$, $\mathit{\beta }$ 등)를 추정해야 합니다. 이는 생각보다 어려운 일이 될 수 있으며, 이전에 언급한 MLE, OLS, Method of Moments (MOM) 등 다양한 model의 estimation 방법이 제안되어 왔습니다. 그런데, 이러한 추정 방법들은 사실상 하나의 추정방정식(estimating equation, 통계 모델의 parameter 추정 방향을 제시하는 모든 식)을 세운 뒤, 그 방정식을 만족하는 $\hat{\mathit{\theta }}$를 찾는 과정으로 해석할 수 있습니다. 예를 들어, MLE에서는 log likelihood를 parameter로 미분한 함수(score function)가 0이 되는 parameter point를 추정하는 과정이었고, OLS에서는 cost function (SSR, 오차 제곱합)을 parameter로 미분한 함수가 0이 되는 parameter point을 찾는 과정이었습니다. M-estimation은 이러한 공통된 개념을 일반화하여 공통된 parameter의 성질을 제시해줍니다.

즉, M-estimation은 다음 과 같은 형태를 지닌 추정 방정식을 세우고, 이를 만족하는 모수의 값을 해로 삼습니다:

$$\sum _{i=1}^n{\psi }_i(\mathit{\theta })=\mathbf{0},$$

여기서 ${\psi }_i(\mathit{\theta })$는 (i)번째 관측치에 대해 정의된 estimating function (e.g. score function), $\mathit{\theta }$는 추정하고자 하는 parameter(위 식에서는 우항이 scalar 0이 아닌 벡터 0이므로 $\mathit{\theta }$ 또한 벡터.)입니다. 위에서 언급한 것처럼 MLS와 OLS, OLS의 일반화 버전인 Non-linear Least Squares 모두 M-estimation입니다. 1, 2장에 걸쳐 이미 익숙하시겠지만 다시 확인해보면, MLE를 M-estimation 형태로 작성해보면 다음과 같고,

$${\psi }_i(\mathit{\theta })=\frac{\partial }{\partial \mathit{\theta }}\log f(Y_i;\mathit{\theta })⟹\sum _{i=1}^n{\psi }_i(\mathit{\theta })=\sum _{i=1}^n\frac{\partial }{\partial \mathit{\theta }}\log f(Y_i;\mathit{\theta })=\mathbf{0}.$$

OLS는 다음과 같습니다. $${\psi }_i(\mathit{\beta })=(Y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }){\mathbf{x}}_i,⟹\sum _{i=1}^n{\mathbf{x}}_i(Y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta })=\mathbf{0}.$$

여기서 OLS의 estimation equation이 다음과 같은 이유는,
$$\begin{gathered} {\mathrm{\nabla }}_{\beta }J(\beta )={\mathrm{\nabla }}_{\beta }\frac{1}{2}(X\beta -y{)}^{\mathrm{\top }}(X\beta -y) \\ =\frac{1}{2}{\mathrm{\nabla }}_{\beta }((X\beta {)}^{\mathrm{\top }}X\beta -(X\beta {)}^{\mathrm{\top }}y-y^{\mathrm{\top }}(X\beta )+y^{\mathrm{\top }}y) \end{gathered}$$

$$=\frac{1}{2}(2X^{\mathrm{\top }}X\beta -2X^{\mathrm{\top }}y)$$

$$=X^{\mathrm{\top }}X\beta -X^{\mathrm{\top }}y$$

$$=\sum _{i=1}^n{\mathbf{x}}_i({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }-Y_i)=\mathbf{0}$$ 에서 유도된 식입니다.

1.2. M-estimation 특징

---

M-estimation의 가장 중요한 특징은 일반성과 확장성입니다. 즉, parameter estimation 문제를

$$\sum _{i=1}^n{\psi }_i(\mathit{\theta })=\mathbf{0}$$

의 해로서 바라보면, 여러 기존 추정법들을 하나의 큰 이론적 틀에서 이해할 수 있고, 이로부터 발생하는 성질들은 해당되는 방법론 모두에 적용됩니다. 이렇게 M-estimation을 강조하여 설명하는 이유는, M-estimation은 아래 두 가지 수렴 이론(Asymptotic theory)을 제공하기 때문입니다.

(1) 적절한 정규성 조건(regularity conditions) 하에서(종속변수의 정규 분포 가정이 아니며, 언급드린 적이 없지만 아주 general한 조건이라고 생각해주시면 됩니다.), 위 M-estimation의 estimating equation의 추정해 $\hat{\mathit{\theta }}$가 참 모수 ${\mathit{\theta }}_0$에 대해 일치성(Consistency)과 점근정규성(Asymptotic Normality)을 가집니다.

(2) 또한, 정규성을 갖는 모수의 점근분포가 중심극한정리(CLT)의 연장선상에 있다고 볼 수 있으며, 그 결과 위 $\hat{\mathit{\theta }}$의 asymptotic Normality에서 parameter의 분산은 Asymptotically&robust하게 추정 가능하고, 이 Robust Estimator의 형태가 샌드위치(sandwich) 형태로 생겼기 때문에 Sandwich Estimation(or)이라고도 부릅니다.

즉, M-estimation으로부터 얻는 의의를 살펴보자면, 우리가 Regression Model의 parameters를 추정하는 과정에서 estimating equation이 위 M-estimation의 형태를 만족한다면, 어떠한 methods를 사용하든 이를 통해 추정한 parameter $\hat{\mathit{\theta }}$는 참 모수 ${\mathit{\theta }}_0$에 대해 consistent함과, robust한 parameter의 분산을 얻을 수 있다는 것입니다. 이제 (1)과 (2)에 대한 수학적 증명을 걸친 뒤, 이들의 의미를 살펴보겠습니다.

1.3. M-estimation의 Asymptotic Normality 증명

---

M-estimation 추정량 $\hat{\theta }$의 점근적 성질을 유도하기 위해 1차 Taylor 전개를 사용합니다. 아래와 같은 M-estimation 추정 방정식 (증명의 편리를 위해 양변에 $\frac{1}{n}$을 나누었으며, 나누지 않아도 똑같이 증명 가능하고, 등식에서 상수 term을 곱하고 나누는 것은 당연히 문제되지 않습니다. ) $$\frac{1}{n}\sum _{i=1}^n{\psi }_i(\hat{\theta })=0$$ 을 참 모수 ${\theta }_0$ Taylor 식으로 전개하면,

$$\frac{1}{n}\sum _{i=1}^n{\psi }_i(\hat{\theta })\approx \frac{1}{n}\sum _{i=1}^n{\psi }_i({\theta }_0)+\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}(\hat{\theta }-{\theta }_0)=0.$$

가 됩니다. 이때 좌항 $\frac{1}{n}\sum _{i=1}^n{\psi }_i(\hat{\theta })$은 우리가 위 estimating equation에서 보았듯이, 이 항이 0이 되도록 하는 parameter를 추정한 결과가 $\hat{\theta }$였기 때문에 당연히 0일 것이고, 따라서 중앙항 (${\theta }_0$에 대한 Taylor 1차 식 전개) 또한 0이 되는 것입니다. 이제 $\theta$에 대한 식을 도출하기 위해 $\frac{1}{n}\sum _{i=1}^n{\psi }_i({\theta }_0)$ term을 넘기고 양변에 $\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}$을 inverse하여(Matrix 이므로) 곱해주고 $\sqrt{n}$을 곱하면,

$$\sqrt{n}(\hat{\theta }-{\theta }_0)\approx -{\left(\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\right)}^{-1}\cdot \frac{1}{\sqrt{n}}\sum _{i=1}^n{\psi }_i({\theta }_0).$$

가 됩니다. 여기서 다음 두 Matrix들을 정의하겠습니다:

• $$\mathbf{A}=\mathbb{E}[-\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}]$$ (2차 도함수 또는 score function의 미분의 기댓값)
• $$\mathbf{B}=\mathbb{E}[{\psi }_i({\theta }_0){\psi }_i({\theta }_0{)}^T]$$ (score function의 분산의 기댓값)

이제 대수의 법칙(LLN)과 중심극한정리(CLT)를 각각 적용하면 다음 두 식을 얻을 수 있습니다.

$$-\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\stackrel{p}{\to }\mathbf{A}$$ $$\frac{1}{\sqrt{n}}\sum _{i=1}^n{\psi }_i({\theta }_0)\stackrel{d}{\to }\mathcal{N}(0,\mathbf{B})$$

각 정리를 간단하게 설명드리자면, 대수의 법칙(LLN, Law of Large Numbers)은 표본 크기 $n$이 충분히 크면, 표본 평균이 모평균에 점근적으로 수렴한다는 법칙으로, 확률 변수 $X_i$가 동일 분포이고 기대값 $\mathbb{E}[X_i]=\mu$를 가지면, 수학적으로 표현하면 아래 식과 같습니다.

$$\frac{1}{n}\sum _{i=1}^n{X}_i\stackrel{p}{\to }\mu .$$

즉, 위 식에서는 $$-\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\approx -\mathbb{E}\left[\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\right]=\mathbf{A}$$ 가 되는 것입니다. 중심극한정리(CLT, Central Limit Theorem)는 독립이고 동일한 분포를 따르는 확률변수들의 표본 평균이 정규 분포를 따른다는 정리로, 분산이 ${\sigma }^2$인 확률변수 $X_i$들에 대해,

$$\frac{1}{\sqrt{n}}\sum _{i=1}^n(X_i-\mu )\stackrel{d}{\to }\mathcal{N}(0,{\sigma }^2).$$ 입니다. (양변에 $\sqrt{n}$이 나눠진 식이 더 친숙하실 겁니다.) 즉, 위 식에서는 ${\psi }_i({\theta }_0)=0$이고, 따라서 $$\mathbf{B}=\mathbb{E}[{\psi }_i({\theta }_0){\psi }_i({\theta }_0{)}^T]=Var({\psi }_i({\theta }_0))$$ 이므로,

$$\frac{1}{\sqrt{n}}\sum _{i=1}^n{\psi }_i({\theta }_0)\stackrel{d}{\to }\mathcal{N}(0,\mathbf{B})$$ 임을 확인할 수 있습니다. 정리하자면, 대수의 법칙이 평균값으로의 수렴을 보장한다면, 중심극한정리는 표본 평균이 정규성을 띤다는 것을 보장하고, 이를 통해 우리가 고려하던 아래 식 $$\sqrt{n}(\hat{\theta }-{\theta }_0)\approx -{\left(\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\right)}^{-1}\cdot \frac{1}{\sqrt{n}}\sum _{i=1}^n{\psi }_i({\theta }_0).$$ 의 우항에 대한 두 정보를 얻을 수 있었습니다. 최종적으로 점근정규성을 보이기 위해선 이 두 수렴하는 분포의 곱을 나타낼 수 있는 Slutsky 정리를 보고, 최종적으로 식을 도출하겠습니다. Slutsky 정리는 두 개의 점근적 확률 분포를 결합하는 방법으로, 만약 $X_n\stackrel{d}{\to }X$ (약한 수렴)과, $Y_n\stackrel{p}{\to }c$ (확률적 수렴)이면,

$$X_i{Y}_i\stackrel{d}{\to }Xc.$$

입니다. 즉, 확률적으로 수렴하는 변수와 분포적으로 수렴하는 변수를 곱하면, 여전히 위 식과 같이 분포적으로 수렴한다는 것이 증명된 정리이고, 위 식에서는

$$\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\stackrel{p}{\to }A$$

$$\frac{1}{\sqrt{n}}\sum _{i=1}^n{\psi }_i({\theta }_0)\stackrel{d}{\to }\mathcal{N}(0,B)$$ 이므로, Slutsky 정리를 사용하면 최종적으로

$$\sqrt{n}(\hat{\theta }-{\theta }_0)\approx -{\left(\frac{1}{n}\sum _{i=1}^n\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\right)}^{-1}\cdot \frac{1}{\sqrt{n}}\sum _{i=1}^n{\psi }_i({\theta }_0).$$

$$\stackrel{d}{\to }{A}^{-1}\mathcal{N}(0,B)=\mathcal{N}(0,A^{-1}B{A}^{-1})$$

입니다.(deteminant한 값은 분산 term에서 제곱된다는 것은 몇 번 보았었습니다.) 결국 M-estimation의 추정을 통해 얻은 $\hat{\mathit{\theta }}$가 참 모수 ${\mathit{\theta }}_0$에 대해 일치성(Consistency)을 갖고, $\hat{\mathit{\theta }}$는 점근정규성(Asymptotic Normality)을 갖으며 그 식은 $\mathcal{N}(0,A^{-1}B{A}^{-1})$입니다. 또한, 이 샌드위치(sandwich) 형태($A^{-1}$ 빵 사이에 껴있는 고기 $B$)처럼 생긴 점근적 분산 식이 바로 Sandwich estimator의 general version입니다. 이제 이 sandwich estimator에 대해 좀더 설명드리겠습니다.

1.4. Sandwich(Robust) Estimator

---

Sandwich(Robust) Estimator의 식은 써보면 다음과 같습니다: $$\widehat{\mathrm{Var}}(\hat{\mathit{\theta }})=A^{-1}B{A}^{-1}$$ $$where,\mathbf{A}=-\mathbb{E}\left[\frac{\partial {\psi }_i({\theta }_0)}{\partial {\theta }^T}\right],\mathbf{B}=\mathbb{E}[{\psi }_i({\theta }_0){\psi }_i({\theta }_0{)}^T]$$

2장에서 GLM case에 대해 log likelihood의 1차 도함수를 score function, 이의 negative 2차 도함수를 Fisher Information matrix라고 언급한 적이 있습니다. 이의 general한 버전이 위와 같으며, 여기에서는 이 score function ${\psi }_i({\theta }_0)$의 분산을 $\mathbf{B}$, 2차 도함수를 $\mathbf{A}$로 표기하고 있습니다.

$\mathbf{A}$는 모형의 곡률(curvature)을, $\mathbf{B}$은 모형의 분산을 반영합니다. 또한, Estimating equation이 log likelihood로부터 MLE 철학으로 나온 parameter라면, $\mathbf{A}=\mathbf{B}$입니다. 이 이유는, 2장에서 2차 도함수가 정의되는 임의의 distribution을 따르는 종속변수 $Y$와 그의 parameter $\theta$에 대해서 $${\ell }^{″}=\frac{d^2\ell }{d{\theta }^2}$$ 가 참임을 보였고, ${\ell }^{″}$은 $\mathbf{A}$, $\frac{d^2\ell }{d{\theta }^2}$는 $\mathbf{B}$와 같기 때문입니다. (${\ell }^{\prime }=\psi$ 이므로.)

즉 철학적으로 해석해보면, Regression Model의 selection이 정확한 경우 Fisher Information 행렬 동일성에 의해 $\mathbf{A}=\mathbf{B}$가 성립하게 되고, 이에 따라 parameter의 분산은 $A^{-1}$ 만으로 추정될 수 있습니다. 그러나 Regression Model이 정확하지 않은 경우, consistent한 parameter estimation을 한다고 하더라도 이 모델의 추정 분산 $A^{-1}$은 더이상 신뢰할 수 없으며, 이때 Sandwich estimtor는 경험적 분산(empirical variance) $\mathbf{B}$를 통해 robust하게 이를 추정할 수 있습니다. 즉, Regression Model의 몇 가지 가정이 의심될 때, 심지어는 의심되지 않더라도 Sandwich estimtor는 robust하게 parameter의 분산을 추정할 수 있는 것입니다.

또한 이전에 스포한대로, 이전 장들에서 다루어 왔던 robust한 parameter variance estimator인 Heteroskedasticity-Consistent SE, Cluster-robust SE는 모두 이 Sandwich estimator의 special한 case입니다.(생김새부터 짐작할 수 있으셨을 겁니다.) LM version에서만 이를 증명한 뒤(GLM 버전도 같습니다.), GLM을 복습하고 GEE, GLMM에 대해서 설명드리겠습니다.

1. Prove HC0 is Sandwich estimator. (LM version)

OLS의 score function은은 위에서 보았듯 다음과 같습니다:

$${\psi }_i(\beta )=x_i(Y_i-x_i^T\beta ).$$ 그렇다면, $A$는 베타로 미분 후 -를 씌워주면 다음과 같이 계산되며,

$$A=\mathbb{E}[-\frac{\partial {\psi }_i(\beta )}{\partial {\beta }^T}]=\mathbb{E}\left[x_i{x}_i^T\right].$$

$A$의 추정치는 결국

$$\hat{A}=\frac{1}{n}\sum _{i=1}^n{x}_i{x}_i^T=\frac{1}{n}{X}^{T}X.$$ 가 됩니다. 마찬가지로 $B$를 계산하면,

$$B=\mathbb{E}[{\psi }_i(\beta ){\psi }_i(\beta {)}^T]=\mathbb{E}[x_i{x}_i^T(Y_i-x_i^T\beta {)}^2].$$

이며, 추정치는

$$\hat{B}=\frac{1}{n}\sum _{i=1}^n{x}_i{x}_i^T{e}_i^2=\frac{1}{n}{X}^T\text{diag}(e_i^2)X.$$ 입니다. 결국

$$\widehat{\text{Var}}(\hat{\beta })={\hat{A}}^{-1}\hat{B}{\hat{A}}^{-1}={\left(\frac{1}{n}{X}^{T}X\right)}^{-1}(\frac{1}{n}{X}^T\text{diag}(e_i^2)X){\left(\frac{1}{n}{X}^{T}X\right)}^{-1}.$$

$$=(X^{T}X{)}^{-1}{X}^T\text{diag}(e_i^2)X(X^{T}X{)}^{-1}.$$

가 되고, 이 식은 1장에서 보았던 HC0의 식과 동일함을 확인할 수 있습니다.

2. Prove Clustered-Robust SE is Sandwich estimator. (LM version)

이 또한 OLS와 같은 환경이므로(LM, cluster가 $g$개로 구성되어 있다고 할 때, score function은 다음과 같습니다:

$${\psi }_g(\beta )=\sum _{i\in g}{x}_i(Y_i-x_i^T\beta ).$$

$A$의 식과 추정치 또한 비슷하게 구해지고,

$$A=\mathbb{E}[-\frac{\partial {\psi }_g(\beta )}{\partial {\beta }^T}]=\mathbb{E}\left[\sum _{i\in g}{x}_i{x}_i^T\right].$$ $$\hat{A}=\frac{1}{n}\sum _{g=1}^G\sum _{i\in g}{x}_i{x}_i^T=\frac{1}{n}{X}^{T}X.$$

$B$도 비슷하게 계산되며, cluster간의 independent는 여전히 가정됩니다.

$$B=\mathbb{E}[{\psi }_g(\beta ){\psi }_g(\beta {)}^T].$$

$$\hat{B}=\frac{1}{n}\sum _{g=1}^G\left(\sum _{i\in g}{x}_i{e}_i\right){\left(\sum _{i\in g}{x}_i{e}_i\right)}^T=\frac{1}{n}\sum _{g=1}^G{X}_g^T{\hat{u}}_g{\hat{u}}_g^T{X}_g.$$

이에 따라 분산의 Sandwich estimator를 구하면

$$\widehat{\text{Var}}(\hat{\beta })={\hat{A}}^{-1}\hat{B}{\hat{A}}^{-1}={\left(\frac{1}{n}{X}^{T}X\right)}^{-1}\left(\frac{1}{n}\sum _{g=1}^G{X}_g^T{\hat{u}}_g{\hat{u}}_g^T{X}_g\right){\left(\frac{1}{n}{X}^{T}X\right)}^{-1}.$$

$$=(X^{T}X{)}^{-1}\left(\sum _{g=1}^G{X}_g^T{\hat{u}}_g{\hat{u}}_g^T{X}_g\right)(X^{T}X{)}^{-1}.$$

이고, 이는 Cluster-robust SE의 식과 동일합니다.

1.5. GLM 복습

---

Generalized Linear Model (GLM)의 모델 식은 다음과 같이 표현됩니다: $$\begin{gathered} g(\mathbb{E}[Y_i|X_i])=g({\mu }_i)={\eta }_i=X_i^T\beta \\ where,Y_i\sim \text{Exponential Family}({\mu }_i,\varphi ). \end{gathered}$$ 이때 $g(\cdot )$은 링크 함수(link function)로 logit, log의 예시를 보았고, ${\mu }_i=\mathbb{E}[Y_i|X_i]$는 반응 변수의 기대값으로 모델의 mapping의 목적이 되는 값(예측하고자 하는 값), $\varphi$ 분산과 관련된 parameter(dispersion parameter) 로 정규 분포의 경우 ${\sigma }^2$였습니다.

간략하게 복습하면 링크 함수(link function) $g(\cdot )$를 통해 $E(Y_i)={\mu }_i$와 ${\eta }_i={\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }$를 연결하고, 분산 함수(variance function) $V({\mu }_i)$를 이용해 $\mathrm{Var}(Y_i)$를 표현하며, 추정방정식(estimating equation)을 세워

$$\sum _{i=1}^n\frac{\partial {\mu }_i}{\partial {\mathit{\beta }}^{\mathrm{\top }}}\frac{(Y_i-{\mu }_i)}{\mathrm{Var}(Y_i)}=\mathbf{0}$$

을 푸는 방식으로 $\hat{\mathit{\beta }}$를 구합니다.

이 해석 또한 M-estimation의 한 사례로 볼 수 있습니다. GLM에서 score 함수(추정방정식)는 ${\psi }_i(\mathit{\beta })=\frac{\partial {\mu }_i}{\partial {\mathit{\beta }}^{\mathrm{\top }}}\frac{(Y_i-{\mu }_i)}{\mathrm{Var}(Y_i)}$ 꼴로 정의되며, 이를 0으로 만드는 $\hat{\mathit{\beta }}$가 우리가 구하고자 하는 파라미터 추정치가 됩니다. 2장에서는 GLM의 parameter $\hat{\mathit{\beta }}$를 추정하는 방법으로 IRLS(Iteratively Reweighted Least Squares)이나 Newton-Raphson/Fisher Scoring을 소개하였으며, 이는 결국 M-estimation에서 구체적으로 어떻게 “estimating equation을 수치적으로 풀어낼지” 알고리즘으로 구현한 예시 중에 하나였다고 이해할 수 있습니다. 또한, 2장에서 예고한대로, 왜 parameter $\hat{\mathit{\beta }}$의 분산이

$$Var(\hat{\mathit{\beta }})=-fisher$$ 라고 했었는지 이제 살펴보면, GLM의 추정 또한 M-estimation에 해당하므로 GLM의 estimating equation을 만족하는 estimator에 대해서 위에서 확인한 점근정규성이 만족함을 알 수 있고, 때문에 consistent한 parameter estimator에 대해서 다음과 같은 robust한 Sandwich 분산을 갖습니다.
$$\mathrm{Var}(\hat{\mathit{\beta }}-{\mathit{\beta }}_0)\approx {(\mathcal{I}({\mathit{\beta }}_0))}^{-1}\mathrm{Var}(\mathbf{S}({\mathit{\beta }}_0)){(\mathcal{I}({\mathit{\beta }}_0))}^{-1}$$

그리고, 계속 보아왔던 것처럼 여기서 $\mathrm{Var}(\mathbf{S}({\mathit{\beta }}_0))=\mathcal{I}({\mathit{\beta }}_0)$가 만족(스코어 함수의 분산의 기댓값과 Fisher information matrix가 같습니다.) 하기 때문에 이를 대입하면 위 식이 아래 처럼 소거되고,

$$\mathrm{Var}(\hat{\mathit{\beta }}-{\mathit{\beta }}_0)\approx \mathcal{I}({\mathit{\beta }}_0{)}^{-1}$$

가 됩니다. 결국 GLM의 모형 기반 분산은 다음과 같습니다:

$$\mathbb{V}a{r}_{\text{모형}}(\hat{\beta })={\mathbf{A}}^{-1}={\left(\sum _{i=1}^n\frac{{\partial }^2\log f(Y_i;\beta )}{\partial \beta \partial {\beta }^T}\right)}^{-1}.$$

또한, 이때 경험적 분포를 고려하여 Sandwich로 추정한 분산은,

$$\mathbb{V}a{r}_{\text{robust}}(\hat{\beta })={\mathbf{A}}^{-1}\mathbf{B}{\mathbf{A}}^{-1},$$

$$where,\mathbf{B}=\sum _{i=1}^n{\psi }_i(\hat{\beta }){\psi }_i(\hat{\beta }{)}^T,{\psi }_i(\beta )=(Y_i-{\mu }_i)x_i/V({\mu }_i).$$

로, 이는 이전에 확인한 HC0의 형태와도 이어집니다.

다시 돌아와서.. (for clustered data)

일반적으로 선형 모델(Linear Model)과 일반화 선형 모델(Generalized Linear Model, GLM)은 독립 동일 분포(i.i.d.)를 가정합니다. 즉, 기존의 GLM은 관측치(observations, data points)들이 서로 독립이며(Independent)일 때 동일한 분산 구조에서 잘 작동합니다. 그러나 학교나 병원 등 군집(클러스터) 단위로 샘플이 묶여 있는, 비슷한 특성을 지닌 대상들을 클러스터(cluster)로 묶은 패널 데이터(panel data)나 동일한 실험 대상(피험자)에게서 반복 측정된 데이터(longitudinal data)의 경우, 같은 cluster(또는 group: 같은 피험자, 같은 단위 등)에 속한 data간에는 correlation이 존재합니다. 때문에 더이상 data들이 독립이 아니게 되며, GLM만으로는 이 상관구조를 모델 자체에서 고려할 수 없기에, GEE와 GLMM 와 같은, 더욱 general한 Regression Model이 개발되었습니다. 이제 아래에서 위 두 model에 대해서 살펴보겠습니다.

2. Generalized Estimating Equation (GEE)

2.1. GEE 정의

---

GEE (Generalized Estimating Equation)는 GLM이 독립성 가정을 전제로 하는 한계마저 뛰어넘어, 군집(Clustered) 자료나 반복측정 자료 등 상관구조가 존재하는 데이터에 적용될 수 있도록 확장한 방법론입니다. 가장 critical하게 다른 점을 보면, GLM은 $\mathrm{Var}(Y_i)=\varphi V({\mu }_i)$로 종속변수의 분산을 표현할 때 diagonal matrix로 두어 data points 간에는 correlation이 없음을 표현하였다면, GEE는 $\mathrm{Var}({\mathbf{Y}}_i)=\varphi {\mathbf{A}}_i^{1/2}{\mathbf{R}}_i(\alpha ){\mathbf{A}}_i^{1/2}$와 같은 working correlation 행렬 ${\mathbf{R}}_i(\alpha )$를 사용하여, 관측치들 간의 상관관계를(반복 측정, 클러스터 내 상관) 모델에 반영합니다. 또한, 이러한 가정을 위해 종속변수의 확률 모델(공동 확률 분포)을 완전히 명시하지 않아도, Quasi-Likelihood(준 우도) 접근법을 통해 점근적(score) 방정식을 확장하여 모델을 적합합니다. 간단히 말하면, GEE는 “평균 모형은 GLM처럼 유지하되, 상관 구조를 적절히 지정하여 군집성이나 반복측정을 고려하자”라는 접근입니다.

LM이나 GLM은 서로 독립적인(i.i.d.) 표본을 가정하여 이를 기반으로 추정하는 반면, GEE에서는 상관 구조(correlation structure) $R(\alpha )$를 추가하여 이러한 독립 가정을 완화하고, 평균 모형과 분산-공분산 구조에 대한 가정을 분리해서 Quasi 형태로 추정합니다. 이 때, Quasi-likelihood에 대한 적용을 짧게 설명하자면, GEE는 확률 모델을 직접 설정(완전한 공동 확률 밀도 함수 명시)하지 않고, GLM의 log likelihood function에 상관구조를 추가하는 형태로 접근합니다. 즉 GLM에서 종속 변수의 Exponential family 가정 -> 독립 가정 후 모든 data point의 확률(likelihood)를 곱해서 얻은 likelihood finction -> 로그 씌워서 log likelihood -> model parameter로 미분한 결과인 score function 순으로 추정 과정을 설명했었다면, GEE는 처음부터 직접적인 종속변수의 가정으로 시작하는 대신 log likelihood에서 시작하고, 이 때 독립이 아님을 고려하기 위해 variance function에 상관 행렬 ${\mathbf{R}}_i(\alpha )$을 추가하여 $\mathrm{Var}({\mathbf{Y}}_i)=\varphi {\mathbf{A}}_i^{1/2}{\mathbf{R}}_i(\alpha ){\mathbf{A}}_i^{1/2}$로 둔 후 추정하는 것입니다. 이러한 접근은, 실제로 종속변수의 완전한 확률적 기반(joint PDF)이 존재하지 않아도, 점근적 성질을 활용하여 일관된 추정량을 얻을 수 있고, 오류 항이 독립적이지 않은 경우에도 GLM과 유사한 방식으로 추정할 수 있는 장점이 있습니다. 마지막으로, GEE는 상관행렬과 Quasi의 개념을 통해 GLM과 같이 data points들을 marginal하게 고려하여 fit하기 때문에 Population-Average GEE(or 모델) 이고, 이는 무작위 효과(Random Effect)를 통해 각 실험 단위(피험자)에 특화된 효과를 추정하는 GLMM(Generalized Linear Mixed Model)과 철학이 다르며,이 GLMM은 Subject-Specific GEE(or 모델)이라고도 부릅니다.

2.2. GEE 수학적 표현 및 추정

---

위에서 언급하였듯, GLM과 동일하게 GEE는 아래와 같은 marginal 모델입니다:
$$g({\mathit{\mu }}_i)={\mathbf{X}}_i\mathit{\beta },$$ 여기서 ${\mathbf{Y}}_i$는 (i)번째 클러스터(또는 피험자)에서 나온 $n_i$개의 관측치 벡터, ${\mathit{\mu }}_i=E({\mathbf{Y}}_i)orE({\mathbf{Y}}_i|X_i)$입니다. 또한, 언급한 대로 Working correlation을 아래와 같이 설정하며,

$$\mathrm{Var}({\mathbf{Y}}_i)=\varphi {\mathbf{A}}_i^{1/2}{\mathbf{R}}_i(\alpha ){\mathbf{A}}_i^{1/2}.$$ 이때 ${\mathbf{A}}_i$가 기존 $V({\mu }_{ij})$의 역할 이었다면 이에 루트를 씌워 A라고 두고 (행렬에서의 square root, 또는 1/2 승은 기존 $V$가 diagonal 이었으므로 이때는 단순히 각 대각 성분을 루트 씌운 값입니다.) 그 사이에 클러스터 당 상관관계를 ${\mathbf{R}}_i(\alpha )$로 표현합니다. 이때 상관행렬 $\mathbf{R}(\alpha )$의 종류로는 크게 아래와 같은 예시들이 있습니다.

(1) Independent (기존 GLM)

$$R(\alpha )=I,V_k=V_k^{\prime }$$

(2) Exchangeable Correlation (동일 상관 구조)

$$R(\alpha )=\left(\begin{array}{ccc}1& \alpha & \alpha \\ \alpha & \ddots & \alpha \\ \alpha & \alpha & 1\end{array}\right)$$

(3) Autoregressive (AR-1)

$$\text{Corr}(y_{ki},y_{kj})={\alpha }^{|i-j|}$$

$$R(\alpha )=\left(\begin{array}{ccccc}1& \alpha & {\alpha }^2& \dots & {\alpha }^{n_k}\\ \alpha & 1& \alpha & \dots & {\alpha }^{n_k-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\alpha }^{n_k}& {\alpha }^{n_k-1}& {\alpha }^{n_k-2}& \dots & 1\end{array}\right)$$

(4) Unstructured Form

$$R(\alpha )=\left(\begin{array}{cccc}1& {\alpha }_1& {\alpha }_2& {\alpha }_3\\ {\alpha }_1& 1& {\alpha }_4& {\alpha }_5\\ {\alpha }_2& {\alpha }_4& 1& {\alpha }_6\\ {\alpha }_3& {\alpha }_5& {\alpha }_6& 1\end{array}\right)$$

이러한 상관행렬 $\mathbf{R}(\alpha )$는 사전에 정의되어야 하므로 분석 대상인 data의 성질에 따라 선정해야 하며, 이러한 관계의 구조를 어떻게 선택하는 지에 따라 분산이 다르게 나오므로, 위에서 다룬 Sandwich를 통한 robust한 추정이 GEE에서 대게 사용됩니다.

GEE’s Estimating Equation

이전에 GLM에서는 다음과 같이 score functions로 부터 estimating equation을 세웠습니다:

$$\mathrm{\Psi }=\sum {\psi }_i=\sum _{i=1}^N\frac{y_i-{\mu }_i}{a(\varphi )V({\mu }_i)}\left(\frac{\partial {\mu }_i}{\partial {\eta }_i}\right)x_i=0$$

이제 이를 GLM 때와 다르게 각 cluster $k$에 대해 $D_k=diag(\left(\frac{\partial {\mu }_i}{\partial {\eta }_i}\right)x_i)$, $V_k={\mathbf{A}}_i^{1/2}{\mathbf{R}}_i(\alpha ){\mathbf{A}}_i^{1/2}$라고 하면,

$$\sum _{k=1}^K\frac{1}{a(\varphi )}{D}_k{V}_k^{-1}(y_k-{\mu }_k)=0$$ 으로 식을 GLM의 score function 으로부터 변형해서 얻을 수 있고, 최종적으로 벡터와 행렬 연산으로 모든 클러스터 $k$에 대해 block diagonal로 한 번에 표현하면($V$), $$\mathrm{\Psi }=\frac{1}{a(\varphi )}D{V}^{-1}(y-\mu )=0$$ 가 되며, $$\mathbf{U}(\mathit{\beta })=\sum _{i=1}^n{\mathbf{D}}_i^{\mathrm{\top }}{\mathbf{V}}_i^{-1}({\mathbf{Y}}_i-{\mathit{\mu }}_i)=\mathbf{0},$$ 로 표현할 수 있습니다. 이는 GLM의 score 함수와 같은 시작이지만, GEE는 ${\mathbf{V}}_i$가 군집/반복측정 상관을 반영하도록 조작합니다. 결국 이 추정방정식을 품으로써 GEE의 추정이 가능할 것입니다.

가장 중요한 GEE에서 분산 term의 변형을 다시 한 번 강조하자면, $V_k$ 는 cluster 별 (Co)variance 행렬로, data가 independent하다면 $V_k$와 이를 모두 합친 $V$가 diagonal matrix가 되지만, 그룹 내 상관을 고려할 경우 $V_k$가 더 이상 diagonal하지 않고, 이에 따라 $V$는 block diagonal matrix 형태를 갖습니다. (block diagonal한 이유는 cluster끼리 독립이고 cluster안은 상관관계가 있는 1차 clustered data에서 다룬 2장의 cluster-robust를 떠올리면 좋을 것 같습니다.)

$$V=\left(\begin{array}{ccc}{V}_1& 0& 0\\ 0& \ddots & 0\\ 0& 0& V_K\end{array}\right)$$

GEE parameter 추정(IRLS)

GEE의 parameter 추정 또한 GLM에서 비롯된 만큼, 이전에 다루었던 방식과 유사한 반복 알고리즘으로 $\hat{\mathit{\beta }}$를 추정할 수 있습니다. 하나의 스텝을 예시로 들어보면,

1. 현재 추정치 ${\hat{\mathit{\beta }}}^{(t)}$에서, 각 클러스터 $i$에 대해 ${\mathbf{D}}_i$ (편미분 행렬), ${\mathit{\mu }}_i=\mu ({\hat{\mathit{\beta }}}^{(t)})$, ${\mathbf{V}}_i$ (working correlation ${\mathbf{R}}_i(\alpha )$, dispersion parameter $\varphi$ 포함)을 계산합니다. 이때, working correlation ${\mathbf{R}}_i(\alpha )$와 $\varphi$도 반복적으로 업데이트됩니다. 예컨대, gee나 geepack 패키지에서는 각 반복 단계에서 잔차(residual)를 기반으로 $\alpha$와 $\varphi$를 재추정하여 새로운 ${\mathbf{V}}_i$를 구합니다.

2. 아래 식을 만족하도록 ${\hat{\mathit{\beta }}}^{(t+1)}$를 업데이트합니다:

   $${\hat{\mathit{\beta }}}^{(t+1)}={\hat{\mathit{\beta }}}^{(t)}-{\left(\sum _{i=1}^n{\mathbf{D}}_i^{\mathrm{\top }}{\mathbf{V}}_i^{-1}{\mathbf{D}}_i\right)}^{-1}\sum _{i=1}^n{\mathbf{D}}_i^{\mathrm{\top }}{\mathbf{V}}_i^{-1}({\mathbf{Y}}_i-{\mathit{\mu }}_i).$$

3. 이전처럼 parameter의 변화량(distance between ${\hat{\mathit{\beta }}}^{(t+1)}$and ${\hat{\mathit{\beta }}}^{(t)}$)가 특정 threshold 아래로 수렴할 때까지 이 과정을 반복합니다.

R에서는 geepack이나 gee 라이브러리에서 내부적으로 이러한 절차를 수행합니다. 1에서 어떻게 잔차로부터 $\alpha$를 추정할 수 있는지 간단하게 예시를 보면 아래와 같습니다. 이는 이전과 원리는 같으며, 상관행렬에서 추정해야 하는 parameter 개수에 따른 자유도를 고려하기 때문에 식이 좀더 복잡해진 것입니다.

잔차는 아래와 같이 계산됩니다(Pearson): $${\hat{r}}_{ki}=\frac{y_{ki}-{\hat{\mu }}_{ki}}{\sqrt{V({\hat{\mu }}_{ki})}},$$ where $y_{ki}$는 observed response for cluster $k$ and observation $i$, ${\hat{\mu }}_{ki}$는 predicted mean for cluster $k$ and observation $i$, $V({\hat{\mu }}_{ki})$는 variance function evaluated at ${\hat{\mu }}_{ki}$. 이제 $n_k=n$라고 가정하면 다음과 같습니다. (이는 클러스터 당 data point 개수가 같다는 가정이며, 이를 만족하지 않아도 식이 복잡해질 뿐 똑같이 계산됩니다.)

(2) Exchangeable Correlation: $$\hat{\alpha }=\frac{1}{a(\varphi )}\sum _{k=1}^K\sum _{i>j}\frac{{\hat{r}}_{ki}{\hat{r}}_{kj}}{K\cdot \frac{1}{2}n(n-1)-p},$$

(3) Autoregressive (AR-1): $$\hat{\alpha }=\frac{1}{a(\varphi )}\sum _{k=1}^K\sum _{i=1}^{n_k-1}\frac{{\hat{r}}_{ki}{\hat{r}}_{k(i+1)}}{K(n-1)-p},$$

(4) Unstructured Form: $${\hat{a}}_{ij}=\frac{1}{a(\varphi )}\sum _{k=1}^K\frac{{\hat{r}}_{ki}{\hat{r}}_{kj}}{K-p},$$

위 식들은 그저 잔차로부터 (co)variance를 추정하는 것일 뿐이고, 분산 term은 degree of freedom을 고려하기 때문에 그저 각각의 상관 행렬 속 미지수(parameter)의 개수에 따른 반영입니다.

2.3. GEE parameter’s Variance

---

GEE의 모수 추정치 ($\hat{\mathit{\beta }}$)도 M-estimation의 범주에 속하므로, 점근 분산은 Sandwich 형태를 갖습니다.

$$\widehat{\mathrm{Var}}(\hat{\mathit{\beta }}{)}_{\text{robust}}={\left(\sum _{i=1}^n{\mathbf{D}}_i^{\mathrm{\top }}{\mathbf{V}}_i^{-1}{\mathbf{D}}_i\right)}^{-1}(\sum _{i=1}^n{\mathbf{D}}_i^{\mathrm{\top }}{\mathbf{V}}_i^{-1}({\mathbf{Y}}_i-{\mathit{\mu }}_i)({\mathbf{Y}}_i-{\mathit{\mu }}_i{)}^{\mathrm{\top }}{\mathbf{V}}_i^{-1}{\mathbf{D}}_i){\left(\sum _{i=1}^n{\mathbf{D}}_i^{\mathrm{\top }}{\mathbf{V}}_i^{-1}{\mathbf{D}}_i\right)}^{-1}.$$

이를 robust 또는 empirical 표준오차라고 하며, 실질적으로 상관구조 (${\mathbf{R}}_i(\alpha )$)가 부정확하게 지정되었을지라도 일관성을 보장해 줍니다. 즉, 위 M-estimation으로 GLM을 해석할 때와 일치하게, Model-based SE 는 ${\mathbf{A}}^{-1}$만을 사용해서 계산하는 것이고, 이는 설정한 working correlation(상관행렬) 가정이 정확하다고 믿을 때이기 때문에, 이를 신뢰할 수 없는 경우 거의 무조건 Robust SE${\mathbf{A}}^{-1}\mathbf{B}{\mathbf{A}}^{-1}$ 를 사용합니다. 이는 R에서 GEE를 계산할 때 기본 SE(model-based)와 robust SE(empirical) 두 가지가 함께 리포팅되는 이유이기도 합니다.

3. Generalized Linear Mixed Model (GLMM)

3.1. GLMM 정의

---

GLMM(Generalized Linear Mixed Model)은, 우리가 이미 익숙하게 다뤄온 GLM(Generalized Linear Model)을 GEE와는 다른 방식으로 (Mixed model) “군집(cluster) 또는 계층적 구조를 가지는 자료”에까지 확장하기 위한 방법론입니다. 즉, GLMM은 이러한 내재된 상관(또는 군집성)을 모델화하기 위해서 고정 효과 + 무작위 효과의 결합으로 모형을 설정합니다. 즉, GLMM은

• 고정 효과(fixed effects): 전체 모집단에 공통적으로 적용되는 회귀계수(예: 전체 평균 경향)에 해당,
• 무작위 효과(random effects): 피험자(또는 군집, 클러스터)별로 달라지는 편차(“개인별 random intercept” 혹은 “개인별 random slope” 등)를 도입

을 둘다 고려하는 모델이며, 즉 “Generalized Linear Model + Linear Mixed Model(Random Effects)”의 결합이라고 요약할 수 있습니다. GEE와 비교하여 이 GLMM은 각 cluster(또는 group)마다 직접적인 고려를 모델에 넣기 때문에(random effect) Subject-Specific 모델(또는 GEE)라고도 불리며, 이는 Population-Average GEE와 대비되는 특징입니다. 무작위 효과는 정규분포로 가정하는 것이 일반적이며, 경우에 따라서는 다른 분포(예: Gamma)로 가정하기도 하고, GLMM에서은 이러한 LMM을 GLM으로 확정한 것이기 때문에 종속변수의 분포를, Exponential Family로 확장합니다.

3.2. LMM 수학적 표현 및 추정

---

GLMM을 이해하기 위해서는 먼저 선형혼합모형(LMM; Linear Mixed Model)을 확실하게 이해할 필요가 있습니다. (이 LMM과 GLMM을 완벽하게 이전처럼 분석하려면 내용이 산만해지기 때문에 여기선 중요한 점을 위주로 짚고, 추가적인 공부가 필요하신 분들은 위키피디아에서 비롯되는 교재 및 논문 내용들을 집중적으로 살펴보시면 좋을 것 같습니다.) LMM은 종속변수 $Y$의 분포가 정규분포라는 전제하에서, 고정 효과(fixed effects)와 무작위 효과(random effects)가 동시에 존재한다고 보는 모형입니다.

LMM 수학적 표현

---

가장 단순한 형태의 LMM(임의절편 모형, random intercept model)을 생각해 보겠습니다. (이때 LMM에서 모형을 나누는 기준은 random effect, 즉 group을 어느 정도로 복잡하게 고려하는 지에 따른 설계의 차이입니다. random effect의 분포, 차원 등을 다양하게 고려할 수 있겠지요.) 예를 들어, $i$번째 클러스터(또는 피험자) 내에서 $j$번째 관측값을 나타내는 $Y_{ij}$를 다음과 같이 모델링합니다:

$$Y_{ij}={\beta }_0+{\beta }_1{X}_{ij}+b_i+{\epsilon }_{ij},i=1,\dots ,K,j=1,\dots ,n_i.$$

이때 ${\beta }_0,{\beta }_1$은 고정 효과(fixed effects) parameter, 즉 모든 클러스터에 공통 적용되는 평균적인 효과)이고, $b_i$는 무작위 효과(random effect) parameter로, 클러스터 $i$마다 서로 다른 절편(intercept) 편차를 갖는 것을 모델링하고 있습니다. ${\epsilon }_{ij}$는 흔히 오차항(residual)으로 간주하고, 대게 ${\epsilon }_{ij}\sim N(0,{\sigma }^2)$로 가정합니다.

추가적으로, 무작위 효과 $b_i$는 다음과 같은 분포로 가정합니다:

$$b_i\sim N(0,{\tau }^2).$$

이는 “(피험자마다) 임의로 달라지는 절편(intercept)”이 정규분포를 따른다는 것을 의미합니다. 모든 $b_i$를 독립으로 가정하면,

$$\mathrm{Var}(b_i)={\tau }^2,\mathrm{Var}({\epsilon }_{ij})={\sigma }^2.$$

결국, 어떤 $Y_{ij}$에 대해서는

$$Y_{ij}={\beta }_0+b_i+{\beta }_1{X}_{ij}+{\epsilon }_{ij},$$

이고,

$$\mathrm{Var}(Y_{ij})={\tau }^2+{\sigma }^2.$$

이먀, 더 일반화 된 모델로 무작위 절편 + 무작위 기울기(random intercept + random slope)를 도입하여 독립변수 $X$에 대해서도 개인별로 기울기가 달라지도록 만들 수 있습니다. 이 경우,

$$Y_{ij}=({\beta }_0+b_{0i})+({\beta }_1+b_{1i})X_{ij}+{\epsilon }_{ij},b_{0i}\sim N(0,{\tau }_{00}),b_{1i}\sim N(0,{\tau }_{11}),\mathrm{Cov}(b_{0i},b_{1i})={\tau }_{01}.$$

가 될 것입니다. 이처럼 무작위 효과를 하나 혹은 여러 개 갖는다는 것은, “클러스터마다 고유하게 발생하는 변동”을 모델에 포함하는 방식으로, LMM은 이러한 방식로 상관구조를 모델링 해낸다고 생각할 수 있습니다. 이를 벡터와 행렬 형태로 표현해보면, 각 클러스터(또는 피험자) $i$에 대해

$${\mathbf{Y}}_i={\mathbf{X}}_i\mathit{\beta }+{\mathbf{Z}}_i{\mathbf{b}}_i+{\mathit{\epsilon }}_i,$$

• ${\mathbf{Y}}_i$: $i$번째 클러스터에서의 $n_i$차원 응답벡터
• ${\mathbf{X}}_i$: $n_i\times p$ 차원의 고정 효과 설계 행렬(fixed effect parameter $\mathit{\beta }$와 매칭)
• ${\mathbf{Z}}_i$: $n_i\times q$ 차원의 무작위 효과 설계 행렬(random effect parameter ${\mathbf{b}}_i$와 매칭)
• ${\mathbf{b}}_i\sim N(\mathbf{0},\mathit{G})$ 이며 $\mathit{G}$는 $q\times q$ 공분산 행렬
• ${\mathit{\epsilon }}_i\sim N(\mathbf{0},{\sigma }^2{\mathbf{I}}_{n_i})$로 일반적으로 가정(독립 동일 분포)

여기서 설계 행렬이란, 1장의 LM에서부터 사용하였지만, data point(observation)당 미리 input으로 지정되는 행렬로, 정확한 의미는 “일련의 개체에 대한 설명 변수 값을 나열한 행렬로 각 행은 개별 개체를 나타내며, 열은 해당 개체에 대한 변수 및 특정 값에 해당한다”입니다. X는 계속 봐왔지만 Z는 이번에 처음 나온 설계 행렬인데, 이는 각 data point마다 해당되는 cluster에는 1, 해당되지 않는 나머지 cluster는 0의 값을 갖는, cluster를 선택하는 스위치 느낌으로, input으로 정해지는 행렬이라고 생각하시면 됩니다.

이 LMM의 (Co) variance matrix는 단순하게 분산 term을 씌우면 random한 (determinant하지 않은) 항만 남아 다음과 같이 계산 될 것입니다:

$$\mathrm{Var}({\mathbf{Y}}_i)={\mathbf{Z}}_i\mathit{G}{\mathbf{Z}}_i^{\mathsf{T}}+{\sigma }^2{\mathbf{I}}_{n_i}.$$

LMM’s parameter 추정(Maximum Likelihood, REML)

---

이 LMM에서 $\mathit{\beta }$, $\mathit{G}$ (또는 ${\tau }^2$ 등), ${\sigma }^2$ 모두 미지수입니다. 이를 추정하기 위해 주로 최대우도추정법(ML; Maximum Likelihood) 또는 제한최대우도추정법(REML; Restricted Maximum Likelihood)을 사용합니다. 각각이 어떻게 계산될지 설명드리면 다음과 같습니다.

(1) With ML(MLE).
$b_i$가 정규분포라는 가정 하에, ${\mathbf{Y}}_i$의 joint distribution도 다변량 정규분포로 표현할 수 있습니다. 모든 $i$에 대해 $b_i$ 또는 ${\mathbf{Y}}_i$가 독립이라 가정하면(cluster간은 독립), 전체 자료의 joint density를 곱해서 likelihood 함수를 정의할 수 있고, 이를 최대화하는 $\hat{\mathit{\beta }}$와 $\hat{\mathit{G}}$, ${\hat{\sigma }}^2$를 찾으면 됩니다. 단점으로는, ML은 $\mathit{\beta }$ 추정에서 편향(bias)이 발생할 수 있어, 표본 크기가 작거나 모형 구조가 복잡해질 때 문제될 수 있습니다. 수식을 중요 부분만 전개해보면, LMM에서 모든 $b_i$를 각각 적분하여(즉 클러스터 마다 적분) 얻은 ${\mathbf{Y}}_i$의 분포는

$${\mathbf{Y}}_i\sim N({\mathbf{X}}_i\beta ,V_i)\text{with}{V}_i={\mathbf{Z}}_i\mathbf{G}{\mathbf{Z}}_i^T+{\sigma }^2{I}_{n_i}.$$

입니다. 이제 $i$번째 클러스터 자료 ${\mathbf{Y}}_i$에 대한 log-likelihood를 계산해보면, 다차원 정규분포이므로 다음과 같이 나옵니다:

$${\ell }_i(\beta ,\mathbf{G},{\sigma }^2)=-\frac{1}{2}[n_i\log (2\pi )+\log |V_i|+({\mathbf{Y}}_i-{\mathbf{X}}_i\beta {)}^T{V}_i^{-1}({\mathbf{Y}}_i-{\mathbf{X}}_i\beta )].$$

이를 전체 $K$개의 cluster에 대해 모두 합하면, 전체 자료에 대한 log-likelihood 함수 $\ell (\beta ,\mathbf{G},{\sigma }^2)$가 도출될 것입니다.

$$\ell (\beta ,\mathbf{G},{\sigma }^2)=\sum _{i=1}^K{\ell }_i(\beta ,\mathbf{G},{\sigma }^2).$$

MLE(${\hat{\beta }}_{ML},{\hat{\mathbf{G}}}_{ML},{\hat{\sigma }}_{ML}^2$)를 구하기 위해서는 위의 log-likelihood를 $\beta ,\mathbf{G},{\sigma }^2$에 대해 미분하여 0이 되게 하는 해를 찾으면 됩니다. 그러나 일반적으로 $\mathbf{G}$와 ${\sigma }^2$에 대한 미분은 해석적으로 단순화하기 어렵고, 또한 $\mathbf{G}$가 양의 정부호(positive definite)가 되어야 하는 제약이 있으므로, 수치적 최적화(EM 알고리즘, Newton-Raphson, Fisher scoring 등)를 사용해야 합니다.

(2) With REML.
이는 ML를 직접 바로 계산하는 대신, 고정 효과 $\mathit{\beta }$와 관련이 없는 term을 이용해 $\mathit{G}$와 ${\sigma }^2$를 먼저 추정한 후 모델을 추정하는 방식입니다. 일반적으로 LMM 을 추정할 때는 REML이 고정 효과 추정에 대한 편의를 줄여주고, 분산 요소에 대한 추정이 좀 더 안정적이기 때문에 ML보다 선호됩니다. 이는 모델에서 무작위 효과를 적분(marginal likelihood)하는 접근을 통해 $\mathit{\beta },\mathit{G},{\sigma }^2$에 대한 우도 함수를 세우고(restricted likeli hood), 이 함수를 최대화하는 방식으로 진행됩니다. 실제 계산은 Iterative 알고리즘(EM 알고리즘, 또는 Fisher scoring, Newton-Raphson 등)을 사용합니다. 식을 보면,

$${\ell }_{REML}(\mathbf{G},{\sigma }^2)=-\frac{1}{2}[\sum _{i=1}^K\log |V_i|+\log |{\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{X}|+(\mathbf{Y}-\mathbf{X}\hat{\beta }{)}^T{\mathbf{V}}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\beta })]+\text{const}.$$

이고, 이때 $\mathbf{V}=\text{blockdiag}(V_1,\dots ,V_K)$, $\hat{\beta }$는 $\mathbf{G},{\sigma }^2$가 주어졌을 때의 일반화최소제곱(GLS) 해입니다. 여기서 $\log |{\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{X}|$가 REML에서 추가로 나타나는 항으로, 이것이 $\beta$를 제거(또는 $\beta$에 무관한 부분만 모아놓음)하여 우선적으로 식을 전개한 효과라고 이해할 수 있습니다. REML에서는 이 ${\ell }_{REML}(\mathbf{G},{\sigma }^2)$를 $\mathbf{G},{\sigma }^2$에 대해 최대화한 뒤, 그 해 ${\hat{\mathbf{G}}}_{REML},{\hat{\sigma }}_{REML}^2$를 이용해 최종적으로 ${\hat{\beta }}_{REML}$ 을 구합니다. REML은 ML보다 fixed effect estimator의 편향 문제가 덜하며, 분산 성분 $\mathbf{G},{\sigma }^2$에 대해 좀 더 안정적인 추정을 제공합니다. 특히 샘플이 작거나 무작위효과 구조가 복잡할 때 일반적으로 더욱 안정적인 REML을 권장하는 편입니다.

이를 통해 얻은 $\hat{\mathit{\beta }}$, $\hat{\mathit{G}}$, ${\hat{\sigma }}^2$는 LMM의 MLE(or REML) 추정치이며, R에서는 lme4 패키지 등에서 이 과정을 내부적으로 수행합니다.

3.3. GLMM의 수학적 표현 및 추정

---

이제 LMM에서 정규 오차항을 일반화하여, 종속변수가 이항, 포아송, 혹은 다른 지수분포족을 따를 수 있도록 확장하면, GLMM으로 이어집니다. GLMM은

$$g({\mu }_{ij})={\mathbf{x}}_{ij}^{\mathsf{T}}\mathit{\beta }+{\mathbf{z}}_{ij}^{\mathsf{T}}{\mathbf{b}}_i$$

$$where,{\mathbf{b}}_i\sim N(\mathbf{0},\mathit{G}),$$

$$Y_{ij}\mid {\mathbf{b}}_i\sim \text{Exponential Family}({\mu }_{ij},\varphi ),$$

의 구조입니다. 직관적으로도 GLMM은 LMM+GLM임을 볼 수 있고, 당연히 $g(\cdot )$는 link function으로, GLM과 마찬가지로 ${\mu }_{ij}=E(Y_{ij}\mid {\mathbf{b}}_i)$를 적절한 링크 함수 $g$로 mapping하는 함수이며, 예시로 binomial case에서 로짓 링크(logit)를 사용하면 $Y_{ij}$는 0 또는 1 값을 가지는 이항분포가 될 수 있고, 이는 $\log ({\mu }_{ij}/(1-{\mu }_{ij}))$를 회귀식을 표현하는 것이었습니다.

이때, 위 식의 likelihood는

$${\mathbf{Y}}_i\mid {\mathbf{b}}_i\sim \prod _{j=1}^{n_i}f(Y_{ij}\mid {\mu }_{ij}({\mathbf{b}}_i)),$$

로 쓸 수 있으며, ${\mathbf{b}}_i$를 적분(marginalizing over ${\mathbf{b}}_i$)하면,

$$p({\mathbf{Y}}_i)=\int \prod _{j=1}^{n_i}f(Y_{ij}\mid {\mu }_{ij}({\mathbf{b}}_i))\phi ({\mathbf{b}}_i)d{\mathbf{b}}_i,$$

가 최종적으로 cluster $i$에 대한 (marginal) 분포를 만들어냅니다.

$$p({\mathbf{Y}}_i)=\int [\prod _{j=1}^{n_i}f(Y_{ij}|{\mu }_{ij}(b_i),\varphi )]\phi (b_i)d{b}_i,\phi (b_i)=\frac{1}{\sqrt{(2\pi {)}^q|\mathbf{G}|}}\exp (-\frac{1}{2}{b}_i^T{\mathbf{G}}^{-1}{b}_i).$$

모든 ${\mathbf{Y}}_i$가 (조건부) 독립이라면, 전체 자료에 대한 marginal likelihood는

$$L(\beta ,\mathbf{G},\varphi )=\prod _{i=1}^K\int \prod _{j=1}^{n_i}f(Y_{ij}|{\mu }_{ij}(b_i),\varphi )\phi (b_i)d{b}_i.$$

입니다. 문제는 ${\mu }_{ij}(b_i)$가 비선형이기 때문에 적분이 closed-form으로 풀리지 않는 경우가 대부분이라는 것이고, 따라서 실제로는 이 적분을 수치적(또는 근사적)으로 계산한 뒤, 그 결과(근사치)를 최대화하여 $\hat{\beta },\hat{\mathbf{G}},\hat{\varphi }$를 구하게 됩니다.

GLMM’s parameter 추정(Marginal Likelihood & Approximation)

---

다시 한 번 언급하자면 문제는 ${\mathbf{b}}_i$가 랜덤효과이므로 이를 적분해야 한다는 것인데, ${\mu }_{ij}({\mathbf{b}}_i)$가 비선형이기 때문에 이 적분이 closed-form으로 표현되지 않는 것이고, 다음과 같은 근사화 기법을 사용하여 계산합니다.

• Laplace Approximation
  ${\mathbf{b}}_i$ 주변에서 2차 근사를 수행하여 적분을 근사화하는 방법입니다. 한 번(1차) 또는 고차(AGQ, Adaptive Gauss-Hermite Quadrature) 버전으로 더 정확하게 시도할 수 있습니다.

• Gauss-Hermite Quadrature
  적분을 수치적(Numerical)으로 가까운 근사값으로 계산합니다. 무작위 효과 차원이 높아질수록 계산량이 기하급수적으로 늘어날 수 있으므로, 실무에서는 차원이 작은 랜덤 효과 구조(예: 랜덤 인터셉트만)에서 자주 사용합니다.

• Penalized Quasi-Likelihood (PQL)
  고전적으로 제안된 근사 기법으로, GLM의 IRLS 절차를 변형하여 무작위효과를 순차적으로 추정합니다. 데이터가 크거나, 근사 정밀도가 크게 중요하지 않은 상황에서 가볍게 쓰일 수 있습니다.

최종적으로, (1) 적분으로 정의된 marginal likelihood를 (2) 수치적 근사화를 통해 (3) 최적화(예: Newton-Raphson, EM 등)하여, $\hat{\mathit{\beta }}$, $\hat{\mathit{G}}$, $\hat{\varphi }$ 등을 찾습니다.

$$\hat{\mathit{\beta }},\hat{\mathit{G}},\hat{\varphi }=\underset{\mathit{\beta },\mathit{G},\varphi }{\mathrm{argmax}}\{\prod _{i=1}^K\int \prod _{j=1}^{n_i}f(Y_{ij}\mid {\mu }_{ij}({\mathbf{b}}_i),\varphi )\phi ({\mathbf{b}}_i)d{\mathbf{b}}_i\}.$$

GLMM vs. GEE

---

이 data간 상관관계를 고려하기 위해 개발된 두 모델을 짧게 정리해보면, GEE는 “Population-Average” 접근으로 군집 내 상관을 working correlation 방식으로 모델링하며, 완전한 joint PDF를 명시하지 않고 Quasi-likelihood처럼 추정하는 기법이었고, GLMM은 “Subject-Specific” 접근으로 군집/클러스터 효과를 무작위 효과로 모델링하여 종속변수를 (조건부) Exponential Family distribution으로 가정하고, 이 likelihood를 marginal하게 적분함으로써 추정합니다.

3.4. GLMM parameter’s Variance

---

마지막으로, GLMM에서의 추정된 파라미터(고정 효과 $\mathit{\beta }$, 무작위 효과 분산-공분산 행렬 $\mathit{G}$ )의 분산 추정(표준 오차, 신뢰구간 등) 방법을 보겠습니다. GLMM의 경우, 근사화하여 최대화한 marginal log-likelihood에서의 헤시안 행렬(Hessian)을 기반으로 고정효과 $\hat{\mathit{\beta }}$ 의 분산을 추정할 수 있습니다. 구체적으로, 아래와 같은 일반적 형식을 취합니다:

$$\widehat{\mathrm{Var}}(\hat{\mathit{\beta }})=[-{\mathrm{\nabla }}_{\mathit{\beta },\mathit{\beta }}^2\ell (\hat{\mathit{\beta }},\hat{\mathit{G}},\hat{\varphi }){]}^{-1},$$

여기서 $\ell$은 GLMM의 (근사) marginal log-likelihood, ${\mathrm{\nabla }}_{\mathit{\beta },\mathit{\beta }}^2$는 고정효과 파라미터 $\mathit{\beta }$에 대한 2차 미분(Hessian)으로, 이 헤시안을 (적절한 수치 방법으로) 근사화하여 얻고, 그 역행렬이 분산 추정의 결과에 해당합니다. 이는 여전히 log likelihood을 통한 추정이기 때문에 Fisher information matrix로 분산을 추정한다고 생각하면 될 것 같습니다. GEE(2장에서)와 마찬가지로, GLMM에서도 모델의 설계에서의 작은 misspecification이 있을 가능성을 고려하여 안정적으로 Sandwich estimator를 통해 추정할 수 있는지 고민할 수 있습니다. GLMM의 경우, 군집 간 독립 이나 군집 내 random effect의 정규성 가정과 같은 가정이 크게 벗어나지 않는다고 믿으면 위 모델 기반(model-based) 추정 분산을 사용하면 되고, 그렇지 않은 “무작위 효과 분포가 정규가 아닐 가능성” 혹은 “link/variance function 형태가 부정확할 가능성” 등을 고려하기 위해 적절한 샌드위치 추정(sandwich-type variance) 기법을 시도할 수도 있습니다. 다만, GEE와 달리 GLMM에서의 robust variance estimation은 쉽게 구현되지 않으며, 근사기법, 부트스트랩(bootstrap) 등을 통해 대안적으로 접근하는 사례도 많습니다.

random effect의 분산-공분산 행렬 $\mathit{G}$ 또한 우도(또는 제한 우도)에서 편미분이 0 조건을 이용하여 추정하지만, 그 표준 오차(불확실성)를 추정하는 과정 역시 (1) 2차 미분, (2) 프로파일(profile) likelihood, (3) 수치적 근사화 등을 거쳐야 합니다.

4. R 코드 예제: GEE, GLMM

아래 R 코드를 복사하여 로컬 환경에서 돌려보세요.

#library(nlme)
#data(Orthodont)  # 치아 성장 데이터 (클러스터: Subject)
#Orthodont$binary <- ifelse(Orthodont$distance > 25, 1, 0)  # 이항 변환
#
## GEE 모델 적합 (Exchangeable 상관 구조)
#library(geepack)
#gee_fit <- geeglm(binary ~ age + Sex,
#                 id = Subject,          # 클러스터 변수
#                 data = Orthodont,
#                 family = binomial,
#                 corstr = "exchangeable")  # 동일 상관 가정
#summary(gee_fit)  # 결과 출력
#
## GLMM 모델 적합 (랜덤 절편 모델)
#library(lme4)
#glmm_fit <- glmer(binary ~ age + Sex + (1|Subject),  # 랜덤 절편
#                 data = Orthodont,
#                 family = binomial)
#summary(glmm_fit)  # 결과 출력

마무리하며

이번 장에서는 M-estimation 개념부터 시작하여, GLM이 어떻게 “estimating equation”의 한 사례로 해석되는지, GEE가 GLM을 확장하여 상관구조를 모델링하고, robust 분산을 제공함으로써 군집/반복측정 데이터를 다루는 과정을, GLMM이 임의효과를 통해 계층적 구조를 명시적으로 모델링하는 방식을 자세히 살펴보았습니다. 그리고 샌드위치 추정량(robust variance) 형태가 M-estimation의 일반 이론에서 비롯된다는 점도 수식과 함께 설명했습니다.

정리하자면, M-estimation은 MLE, OLS, GEE, GLMM 모두를 포괄하는 추정 이론적 틀로서, 샌드위치 분산은 그 점근 정규성(Asymptotic Normality)의 결과물이며, GEE는 marginal mean에 주목하고 robust한 표준오차를 산출해주는 반면, GLMM은 임의효과를 통해 개체별(군집별) 차이를 직접 모델링합니다. 실제 데이터 분석에서는 연구 목적(개체별 효과 추정 vs 전체 평균 효과 추정), 데이터 특성(정확한 상관 구조 가정 vs 모형 가정의 유연성) 등을 종합하여 GEE와 GLMM 중 적절한 접근을 택하거나 비교하는 것이 중요합니다. 사실 Regression Model에는 이번 블로그 “Exploring Regression Models for Regression Analysis”에서 다룬 모델들을 제외하고도 아주 다양한 철학과 수식을 가진 모델들이 있습니다. 다만 여기서는 의학 분석에서 자주 사용되는 모델을 다루었으며, 이를 어느 정도 이해하셨다면 이외의 모델을 이해하는 데에 부족함이 없을 것이라고 생각합니다.