1 들어가며 — GEE 의 모형과 추정 깊이

Ch.8 Overview 와 § 8.1-8.2 GEE 의 토대 sub-post 가 GEE 의 자리와 GLM 토대를 다뤘다. 본 sub-post 는 §8.3·§8.4 — GEE 의 모형 정의 와 추정 절차 — 를 깊이 다룬다.

절	내용	본 sub-post 강조
§ 8.3	GEE 모형 + 5 작동 상관	행렬 구체, robust 성질 정량
§ 8.4	GEE 추정 + Sandwich	IRLS 단계, Naive vs Robust SE

한 줄 요약

“§ 8.3 = GLM 에 작동 상관 추가 (모형 정의), § 8.4 = IRLS 로 추정 + Sandwich 로 SE. 두 절이 GEE 의 작동 메커니즘 전체.”

2 § 8.3 — GEE 모형 정의

2.1 4 가지 명세 (GLM 3 + 작동 상관)

GEE Specifications (식 8.12-8.14, +R)

Ch.8 § 8.2 GLM 복습 의 3 가지 명세에 작동 상관 추가.

1. Linear Predictor (식 8.12):

\[ \eta_{ij} = x_{ij}^\top \beta \tag{8.12} \]

2. Link Function (식 8.13):

\[ g(\mu_{ij}) = \eta_{ij} \tag{8.13} \]

3. Variance Function (식 8.14):

\[ V(y_{ij}) = \phi v(\mu_{ij}) \tag{8.14} \]

4. Working Correlation Matrix (추가):

\[ R_i(a) = \text{Corr}(y_i \mid X_i) \]

\(n \times n\) 상관 행렬 (\(n\) = 시점 수).
\(a\): association parameter 벡터.
“작동” (working) — 분석자 가정, 진짜와 같을 필요 없음.
모든 피험자에서 같은 \(R\) 사용 (시점 다르면 부분 행렬 추출).

2.2 Robust 성질 — 정량적 분석

GEE 의 가장 매력적 성질

핵심 정리 (Liang & Zeger 1986):

작동 상관 \(R_i(a)\) 가 잘못 명세돼도 (true correlation 과 다르더라도): 1. 회귀 계수 \(\hat\beta\) 가 일치 (consistent) 추정. 2. \(\hat\beta\) 의 분포가 점근 정규. 3. Sandwich SE 도 일치 추정 (§ 8.4 에서).

이 성질의 수학적 핵심: 추정 방정식 (식 8.17) 이 평균 함수에만 의존하고 작동 상관은 가중치로만 작용. 일치성 (consistency) 이 평균 함수의 정확성에서 나옴.

효율성 손실의 정량 분석

작동 상관 misspecification 의 효율성 (statistical power) 손실:

True correlation	Working correlation	\(\hat\beta\) 분산 비율 (efficiency)
AR(1)	AR(1)	1.00 (최적)
AR(1)	Independence	약 0.85~0.95
AR(1)	Exchangeable	약 0.90~0.95
Exchangeable	AR(1)	약 0.90~0.95
Exchangeable	Independence	약 0.80~0.90
Unstructured	AR(1)	약 0.85~0.95

→ 효율성 손실은 보통 5~20% 수준. 통계적 검정 결과에는 영향 작음 — 진짜 효과가 있으면 잘못된 작동 상관에서도 잡힘.

표본 크기의 영향: 표본 클수록 효율성 손실 작음. 종단 임상 시험 (\(N \geq 100\)) 에서 작동 상관 선택의 부담은 가벼움.

실무 권고 — 작동 상관 선택

EDA 우선: 관측 잔차의 상관 행렬 그림.
패턴 식별:
- lag 따라 매끄럽게 감소 → AR(1)
- lag 무관 일정 → Exchangeable
- 비단조 → Unstructured (시점 적을 때) 또는 m-dependent
Robust SE 항상 보고: Naive 와 비교로 작동 상관 misspecification 진단.
두 모형 적합: 가장 적합 vs 단순 (Independence) — 결과 비슷하면 단순 모형.
QIC 보조: 모형 비교 (LR 부재).

3 § 8.3.1 — 5 작동 상관 행렬 구체

3.1 Independence

\(R_i(a) = I\)

가장 단순한 형태:

\[ R_i = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix} \]

모수 수: 0 (자유 모수 없음).
가정: 시점들이 서로 독립.

Independence 의 함정

종단 데이터에 보통 비현실적 — 같은 환자의 시점들이 강하게 상관.

효율성 손실: 종단 이항 outcome 에서 시변 공변량 분석 시 큰 손실 (Fitzmaurice 1995).

예외 — 장점이 있는 경우:

시변 공변량 (Pepe & Anderson 1994): 일부 회귀 계수의 식별 가능성.
모수 절약: \(a\) 추정 안 해도 됨 → 작은 표본에서 안정.
시점 매우 다름 (서로 다른 환자가 다른 시점 측정): Independence 가 자연.

3.2 Exchangeable

모든 lag 의 상관 동일

\[ R_i(\rho) = \begin{bmatrix} 1 & \rho & \rho & \cdots & \rho \\ \rho & 1 & \rho & \cdots & \rho \\ \rho & \rho & 1 & \cdots & \rho \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \rho & \cdots & 1 \end{bmatrix} \]

모수 수: 1 (\(\rho\)).
가정: 모든 시점 쌍이 같은 상관.

이는 랜덤 절편 MRM 또는 CS-CPM 와 같은 가정 — § 6.2.1 의 CS 와 직접 매핑.

Exchangeable 의 직관

“모든 시점이 동등한 친밀도” — 시간 lag 무시.

적합 시나리오:

클러스터 데이터 (학교 안 학생, 병원 안 환자) — 시간 개념 없음.
짧은 추적 (시점 2~3 개) — lag 변동 적어 절약.
무작위 효과 모형 가정 — 랜덤 절편이 시점 간 동등 상관 만듦.

한계: 종단 데이터에 보통 부족. Bock WPSS 같은 임상 데이터에서 lag 따라 명확히 감소 → Exchangeable 이 부적합 (§ 6.4 분석 참조).

3.3 AR(1)

지수 감쇠 상관

\[ R_i(\rho) = \rho^{|j-j'|}, \quad j, j' = 1, \ldots, n \]

행렬:

\[ R_i = \begin{bmatrix} 1 & \rho & \rho^2 & \cdots & \rho^{n-1} \\ \rho & 1 & \rho & \cdots & \rho^{n-2} \\ \rho^2 & \rho & 1 & \cdots & \rho^{n-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho^{n-1} & \rho^{n-2} & \rho^{n-3} & \cdots & 1 \end{bmatrix} \]

모수 수: 1 (\(\rho\)).
가정: 정상성, lag 의 지수 함수.

§ 6.2.2 의 AR(1) CPM 또는 § 7.2.1 의 AR(1) AC 오차 와 같은 형태.

AR(1) 의 실무 가치

적합 시나리오:

등간격 추적 (매주, 매월).
lag 따라 매끄럽게 감소하는 상관.
단기 momentum (어제 영향이 오늘에 carry over).

예시 — Bock WPSS 데이터: lag-1 ≈ 0.85, lag-5 ≈ 0.5 — AR(1) 가 일부 적합 (단 시점별 분산 변동으로 UN 보다 부족).

모수 절약성: 1 모수로 모든 lag 결정 — 표본 작을 때 안정.

3.4 m-dependent (Banded)

lag \(\leq m\) 만 비영

\[ R_i(a)_{jj'} = \begin{cases} 1 & \text{if } j = j' \\ \rho_{|j-j'|} & \text{if } 0 < |j-j'| \leq m \\ 0 & \text{if } |j-j'| > m \end{cases} \]

예 — \(n = 6\), \(m = 2\):

\[ R_i = \begin{bmatrix} 1 & \rho_1 & \rho_2 & 0 & 0 & 0 \\ \rho_1 & 1 & \rho_1 & \rho_2 & 0 & 0 \\ \rho_2 & \rho_1 & 1 & \rho_1 & \rho_2 & 0 \\ 0 & \rho_2 & \rho_1 & 1 & \rho_1 & \rho_2 \\ 0 & 0 & \rho_2 & \rho_1 & 1 & \rho_1 \\ 0 & 0 & 0 & \rho_2 & \rho_1 & 1 \end{bmatrix} \]

모수 수: \(m\) (\(\rho_1, \rho_2, \ldots, \rho_m\)).
가정: 정상성 (같은 lag 동일), 고차 lag 0.

m-dependent 의 의미

Toeplitz CPM 의 절약 형태. 가장 일반 (m = n-1) 일 때 full Toeplitz.

적합 시나리오:

시점 많은 종단 데이터 (\(n \geq 6\)).
단기 의존만 — 멀리 떨어진 시점은 거의 무관.
ACF 가 lag \(m\) 에서 0 으로 떨어지는 패턴.

모수 trade-off: \(m\) 클수록 유연성 + 모수. m=1 (lag-1만) ≈ MA(1)-like, m=n-1 ≈ full Toeplitz.

Hedeker §8.3.1 의 표기: m=2 면 “Toeplitz(3)” (CPM 관행에서 분산 포함 카운트).

3.5 Unstructured (Unspecified)

모든 상관 자유

\[ R_i = \begin{bmatrix} 1 & \rho_{12} & \rho_{13} & \cdots & \rho_{1n} \\ \rho_{12} & 1 & \rho_{23} & \cdots & \rho_{2n} \\ \rho_{13} & \rho_{23} & 1 & \cdots & \rho_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho_{1n} & \rho_{2n} & \rho_{3n} & \cdots & 1 \end{bmatrix} \]

모수 수: \(n(n-1)/2\) — 대각 위 모든 원소.
가정: 없음. 가장 유연.

Unstructured 의 모수 폭발

\(n\)	\(n(n-1)/2\)
3	3
4	6
6	15
10	45
20	190

→ 시점 많을 때 추정 불안정. 시점 적고 표본 클 때 가장 적합.

적합 시나리오:

시점 \(\leq 6\), 표본 \(\geq 100\).
분산 구조 자체가 연구 질문.
가정 위반 위험 회피 (보수적 분석).

3.6 5 작동 상관 비교 표

형태	\(R_i(a)\)	모수 수	CPM 대응	적합 시나리오
Independence	\(I\)	0	\(\sigma^2 I\)	시변 공변량, 절약
Exchangeable	모두 \(\rho\)	1	CS	클러스터, 짧은 추적
AR(1)	\(\rho^{\|j-j'\|}\)	1	AR(1)	등간격, 지수 감쇠
m-dependent	lag ≤ m 비영	\(m\)	Toeplitz(s)	단기 의존, 비단조
Unstructured	모든 \(\rho_{jj'}\)	\(n(n-1)/2\)	UN	시점 적음, 표본 큼

CPM 와 동일 5 구조 + Independence

GEE 의 5 작동 상관과 CPM 의 5 구조 가 거의 일치 — 차이는 Independence 추가 (CPM 에는 RE 가 그 자리).

→ GEE 사용자는 CPM 의 분산 구조 이해를 그대로 활용 가능.

4 § 8.4 — GEE 추정

4.1 작동 분산-공분산 (식 8.15)

\(V_i(a) = \phi A_i^{1/2} R_i(a) A_i^{1/2}\)

GLM 의 분산 함수 \(v(\mu_{ij})\) 와 작동 상관 \(R_i(a)\) 를 결합:

\[ V_i(a) = \phi A_i^{1/2} R_i(a) A_i^{1/2} \tag{8.15} \]

여기서 \(A_i\) 는 \(n_i \times n_i\) 대각 행렬, \(j\) 번째 대각 원소가 \(v(\mu_{ij})\).

특수 경우 — 정규 + 동질 분산 (\(v(\mu_{ij}) = 1\)):

\[ V_i(a) = \phi R_i(a) \tag{8.16} \]

분산 함수 + 작동 상관의 결합 의미

\(V_i\) 의 두 부분:

분산 부분 (\(A_i\)): 시점별 marginal 분산 — GLM 분포가 결정.
- 정규: \(v(\mu) = 1\) (분산 동일).
- 이항: \(v(\mu) = \mu(1-\mu)\) (평균에 따라 변동).
- Poisson: \(v(\mu) = \mu\) (평균과 분산 같음).
상관 부분 (\(R_i\)): 시점 간 상관 — 작동 상관 가정.

분산 자체는 GLM 의 분포가 결정 (지수족 이론), 상관만 작동 상관 가정. 이게 GEE 가 비정규 반응까지 자연 처리하는 메커니즘.

4.2 추정 방정식 (식 8.17)

GEE 의 핵심 식

\[ \sum_{i=1}^N D_i^\top [V_i(\hat a)]^{-1} (y_i - \mu_i) = 0 \tag{8.17} \]

\(D_i = \partial \mu_i / \partial \beta\): \(n_i \times p\) Jacobian.
\(\hat a\): association 모수의 일치 추정값.
\(V_i\): 작동 분산-공분산 (식 8.15).

GLM 의 추정 방정식 (식 8.10) 의 직접 일반화 — \(V_i\) 가 작동 상관 포함으로 확장.

식 (8.17) 의 해석

세 항으로 분해:

\(D_i^\top\): 평균이 \(\beta\) 에 어떻게 의존 (Jacobian).
\(V_i^{-1}\): 작동 분산의 역수 (가중치).
\(y_i - \mu_i\): 잔차 벡터.

→ “가중 잔차의 합이 0 인 \(\beta\) 를 찾는다.”

GLM 의 같은 형태이지만 가중치 행렬 \(V_i\) 가 작동 상관을 포함 — 이것이 GEE 의 추가.

4.3 정규 사례 — WLS 형태 (식 8.19)

정규 + identity link 의 닫힌 형태

\(\mu_i = X_i\beta\), \(D_i = X_i\), \(V_i = R_i(\hat a)\) (정규 + \(\phi = 1\)):

\[ \hat\beta = \left[\sum X_i^\top R_i(\hat a)^{-1} X_i\right]^{-1} \left[\sum X_i^\top R_i(\hat a)^{-1} y_i\right] \tag{8.19} \]

→ Weighted Least Squares (WLS) 의 형태. 가중치가 \(R_i^{-1}\).

차이점: 표준 WLS 는 가중치를 알지만, GEE 에서는 가중치가 모수 (\(a\)) 에 의존 → 반복 추정 필요.

4.4 IRLS — GEE 특화 단계

두 단계 반복 (Hedeker §8.4)

Step 1: \(\hat a, \hat\phi\) 고정 → IRLS 로 \(\hat\beta\) 갱신.

Step 2: \(\hat\beta\) 고정 → Pearson 잔차로 \(\hat a, \hat\phi\) 갱신.

수렴까지 1 ↔︎ 2 반복.

Step 1 — IRLS 로 \(\hat\beta\) 갱신

§ 8.2 의 GLM IRLS 와 같은 알고리즘, 단 가중치 행렬에 작동 상관 포함:

각 반복 \(t \to t+1\):

현재 \(\hat\beta^{(t)}\) 에서 \(\hat\mu_i^{(t)}\), \(V_i^{(t)}\) 계산.
작업 응답: \(z_i^{(t)} = \eta_i^{(t)} + (y_i - \mu_i^{(t)}) g'(\mu_i^{(t)})\).
WLS 갱신: \[ \hat\beta^{(t+1)} = \left[\sum X_i^\top V_i^{(t),-1} X_i\right]^{-1} \left[\sum X_i^\top V_i^{(t),-1} z^{(t)}\right] \]

→ GLM IRLS 의 자연 확장. 작동 상관이 가중 행렬에 들어간다.

4.5 Step 2 — Pearson 잔차 (식 8.20)

\(\hat a, \hat\phi\) 갱신

각 시점의 표준화 잔차:

\[ r_{ij} = \frac{y_{ij} - \hat\mu_{ij}}{\sqrt{[V_i(\hat a)]_{jj}}} \tag{8.20} \]

분자: 잔차.
분모: 시점 \(j\) 의 marginal 표준편차.
결과: 평균 0, 분산 1 의 표준화된 양.

이 Pearson 잔차로 association 모수 \(a\) 와 산포 \(\phi\) 추정:

\(\phi\) 추정:

\[ \hat\phi = \frac{1}{N - p} \sum_{i,j} r_{ij}^2 \]

\(a\) 추정 (작동 상관 형태별):

Exchangeable: \(\hat\rho = \frac{1}{N \binom{n}{2}} \sum_i \sum_{j < j'} r_{ij} r_{ij'}\).
AR(1): \(\hat\rho = \frac{\sum r_{ij} r_{i,j+1}}{\sum r_{ij}^2}\) (lag-1 평균).
Unstructured: \(\hat\rho_{jj'} = \frac{1}{N} \sum_i r_{ij} r_{ij'}\).

Pearson 잔차의 직관

Pearson 잔차는 분포의 marginal 분산을 제거 한 표준화된 잔차.

이항 GLM (\(V(y) = \mu(1-\mu)\)): 평균 0.5 부근의 잔차가 가장 큼 → 표준화로 비교 가능한 단위로 변환.

Poisson GLM (\(V(y) = \mu\)): 평균 큰 시점의 잔차가 자연스럽게 큼 → 표준화로 균질화.

→ 다른 시점의 잔차들을 곱해서 lag 별 상관 추정 가능 (분포 효과 제거 후).

5 § 8.4 — Naive vs Robust SE

5.1 두 가지 분산-공분산 추정

수렴 후 \(\hat\beta\) 의 분산-공분산은 두 형태로 계산.

Naive (Model-Based, 식 8.21)

작동 상관이 진짜 상관과 같다고 가정:

\[ V_{\text{naive}}(\hat\beta) = \left[\sum_i D_i^\top \hat V_i^{-1} D_i\right]^{-1} \tag{8.21} \]

형태: 일반 ML 의 정보 행렬 역수.
가정: 모형이 정확.
\(\hat V_i\) = \(V_i(\hat a)\) — 작동 상관 추정값 사용.

Robust (Sandwich, 식 8.22)

작동 상관이 잘못돼도 일치 추정:

\[ V_{\text{robust}}(\hat\beta) = M_0^{-1} M_1 M_0^{-1} \tag{8.22} \]

여기서:

\[ M_0 = \sum_i D_i^\top \hat V_i^{-1} D_i \]

\[ M_1 = \sum_i D_i^\top \hat V_i^{-1} (y_i - \hat\mu_i)(y_i - \hat\mu_i)^\top \hat V_i^{-1} D_i \]

→ Royall (1986) 의 robust 분산 추정량.

5.2 “Sandwich” 의 직관

빵-속재료-빵 구조

식 (8.22) 의 형태:

\[ \underbrace{M_0^{-1}}_{\text{빵 (위)}} \, \underbrace{M_1}_{\text{속재료}} \, \underbrace{M_0^{-1}}_{\text{빵 (아래)}} \]

빵 (\(M_0^{-1}\)): 작동 상관 기반 분산 (model-based). 같은 형태가 위·아래.
속재료 (\(M_1\)): 잔차 기반 empirical 분산. 진짜 상관 정보.

왜 sandwich 가 robust 한가:

작동 상관이 진짜와 같으면:

\[ \hat V_i = (y_i - \hat\mu_i)(y_i - \hat\mu_i)^\top \implies M_1 = M_0 \]

따라서:

\[ V_{\text{robust}} = M_0^{-1} M_0 M_0^{-1} = M_0^{-1} = V_{\text{naive}} \]

작동 상관이 다르면 \(M_1\) 이 잔차의 진짜 분산 정보를 흡수 → \(V_{\text{robust}}\) 가 일치 추정.

5.3 Sandwich estimator 의 수학적 유도 sketch

왜 식 (8.22) 가 일치 추정인가

표준 ML 이론 (Sandwich 형태):

\[ \sqrt{N} (\hat\beta - \beta_0) \xrightarrow{d} \mathcal{N}(0, M_0^{-1} \Omega_0 M_0^{-1}) \]

여기서:

\(M_0\): 정보 행렬의 기댓값.
\(\Omega_0\): score function 의 분산.

올바르게 명세된 모형에서: \(\Omega_0 = M_0\) (Bartlett 등식) → \(V = M_0^{-1}\) (Naive).

Misspecification 시: \(\Omega_0 \neq M_0\) → Sandwich 가 일반.

GEE 에서:

\(M_0 = E[D^\top V^{-1} D]\)
\(\Omega_0 = E[D^\top V^{-1} (y-\mu)(y-\mu)^\top V^{-1} D]\)
표본 추정: \(\hat M_0\), \(\hat M_1 = \hat\Omega_0\).

→ \(V_{\text{robust}} = \hat M_0^{-1} \hat M_1 \hat M_0^{-1}\) 가 진짜 분산의 일치 추정 (Royall 1986).

5.4 Naive vs Robust 비교 — 진단 도구

두 SE 의 차이가 모형 검진

작동 상관이 정확하면 두 SE 가 비슷, 잘못되면 차이.

두 SE 의 관계	의미
\(\text{SE}_{\text{robust}} \approx \text{SE}_{\text{naive}}\)	작동 상관이 진짜 상관과 비슷 — 모형 OK.
\(\text{SE}_{\text{robust}} > \text{SE}_{\text{naive}}\)	작동 상관이 너무 단순 — 다른 작동 상관 시도.
\(\text{SE}_{\text{robust}} < \text{SE}_{\text{naive}}\)	작동 상관이 너무 복잡 — 데이터에 비해 과적합 가능.

실무 절차:

항상 Robust SE 보고 — 분석의 표준.
Naive 와 비교 — 진단 정보로 사용.
큰 차이면 작동 상관 재검토 — 다른 형태 시도.
Naive·Robust 모두 보고 — 보수적 실무.

5.5 작은 표본의 함정

Sandwich SE 의 작은 표본 편향

Sandwich estimator 는 점근적 (asymptotic) 일치. 표본 작으면 (피험자 수 < 30~50) 하향 편향:

\[ E[V_{\text{robust}}] < V(\hat\beta) \]

→ SE 너무 작음 → 거짓 유의 결과 가능.

보정 방법:

MD (Mancl & DeRouen 2001) 보정: \(M_1\) 의 표본 보정.
KC (Kauermann & Carroll 2001) 보정: 자유도 조정.
FG (Fay & Graubard 2001) 보정: t 분포 사용.

R geepack 의 geeglm 에는 일부 보정 옵션. SAS PROC GENMOD 도 마찬가지.

실무 권고: 피험자 수 < 30 이면 Sandwich SE 만 신뢰하지 말고 보정된 SE 또는 bootstrap 사용.

6 IRLS 와 Sandwich 의 통합 알고리즘

GEE 추정의 전체 흐름

초기값: \(\hat\beta^{(0)}\) — 보통 GLM (Independence) 추정.
반복:
1. 현재 \(\hat\beta^{(t)}\) 에서 \(\hat\mu, \hat V_i\) 계산.
2. Step A — IRLS: WLS 로 \(\hat\beta^{(t+1)}\) 갱신.
3. Step B — Pearson 잔차: \(\hat a^{(t+1)}, \hat\phi^{(t+1)}\) 갱신.
4. 수렴 검사: \(\|\hat\beta^{(t+1)} - \hat\beta^{(t)}\| < \epsilon\).
수렴 후 SE 계산:
- Naive (식 8.21).
- Robust Sandwich (식 8.22).
검정:
- Wald: \(W = \hat\beta_j^2 / \widehat{\text{Var}}(\hat\beta_j)\), robust SE 사용.
- Generalized Wald (다중 모수).
- Score 검정.

7 코드 예시

7.1 Step 1: 5 작동 상관 행렬 직접 구성

import numpy as np


def correlation_independence(n: int) -> np.ndarray:
    """Independence 작동 상관"""
    return np.eye(n)


def correlation_exchangeable(rho: float, n: int) -> np.ndarray:
    """Exchangeable: 모든 lag 동일 ρ"""
    R = np.full((n, n), rho)
    np.fill_diagonal(R, 1.0)
    return R


def correlation_ar1(rho: float, n: int) -> np.ndarray:
    """AR(1): ρ^|j-j'|"""
    idx = np.arange(n)
    lag = np.abs(idx[:, None] - idx[None, :])
    return rho ** lag


def correlation_m_dependent(rho_vec: np.ndarray, m: int, n: int) -> np.ndarray:
    """m-dependent: lag ≤ m 만 비영"""
    R = np.eye(n)
    for k in range(1, min(m + 1, n)):
        for i in range(n - k):
            R[i, i + k] = rho_vec[k - 1]
            R[i + k, i] = rho_vec[k - 1]
    return R


def correlation_unstructured(rho_dict: dict, n: int) -> np.ndarray:
    """Unstructured: 모든 ρ_jj' 자유

    rho_dict: {(j, j'): rho_jj'} for j < j'.
    """
    R = np.eye(n)
    for (j, jp), rho in rho_dict.items():
        R[j, jp] = rho
        R[jp, j] = rho
    return R


# 비교 예시
n = 5
print("Independence:")
print(correlation_independence(n))

print("\nExchangeable (ρ=0.5):")
print(correlation_exchangeable(0.5, n).round(2))

print("\nAR(1) (ρ=0.7):")
print(correlation_ar1(0.7, n).round(3))

print("\n2-dependent (ρ_1=0.6, ρ_2=0.3):")
print(correlation_m_dependent(np.array([0.6, 0.3]), m=2, n=n).round(2))

7.2 Step 2: GEE Sandwich 추정 직접 구현

import numpy as np


def gee_sandwich_se(X: np.ndarray, y: np.ndarray, groups: np.ndarray,
                    beta: np.ndarray, working_corr: callable,
                    family: str = "binomial") -> dict:
    """GEE Sandwich SE 직접 계산

    X:           (N_total, p) 디자인.
    y:           (N_total,) 반응.
    groups:      (N_total,) 군집 ID.
    beta:        (p,) 추정 회귀 계수.
    working_corr: function(n_i) → R_i.
    family:      GLM family.
    """
    p = X.shape[1]
    M0 = np.zeros((p, p))
    M1 = np.zeros((p, p))

    for g in np.unique(groups):
        idx = groups == g
        X_i = X[idx]
        y_i = y[idx]
        n_i = len(y_i)

        # 평균과 분산 함수
        eta = X_i @ beta
        if family == "binomial":
            mu = 1 / (1 + np.exp(-eta))
            v = mu * (1 - mu)
            D_i = X_i * (mu * (1 - mu))[:, None]  # ∂μ/∂β
        elif family == "poisson":
            mu = np.exp(eta)
            v = mu
            D_i = X_i * mu[:, None]
        elif family == "normal":
            mu = eta
            v = np.ones(n_i)
            D_i = X_i

        # 작동 분산-공분산: V_i = A^(1/2) R A^(1/2)
        A_sqrt = np.diag(np.sqrt(v))
        R = working_corr(n_i)
        V_i = A_sqrt @ R @ A_sqrt
        V_i_inv = np.linalg.inv(V_i)

        # M_0 갱신
        M0 += D_i.T @ V_i_inv @ D_i

        # M_1 갱신
        resid = y_i - mu
        M1 += D_i.T @ V_i_inv @ np.outer(resid, resid) @ V_i_inv @ D_i

    M0_inv = np.linalg.inv(M0)

    # Naive vs Robust SE
    naive_cov = M0_inv
    robust_cov = M0_inv @ M1 @ M0_inv

    return {
        "naive_se": np.sqrt(np.diag(naive_cov)),
        "robust_se": np.sqrt(np.diag(robust_cov)),
        "naive_cov": naive_cov,
        "robust_cov": robust_cov,
    }

검증 포인트

위 함수의 출력에서:

작동 상관이 진짜에 가까우면: naive_se ≈ robust_se.
작동 상관 misspecification: robust_se ≠ naive_se.

R geepack::geeglm() 의 vbeta, vbeta.naiv 와 비교 가능.

7.3 Step 3: 작동 상관 비교 진단

import pandas as pd
from statsmodels.genmod.generalized_estimating_equations import GEE
from statsmodels.genmod.cov_struct import (
    Independence, Exchangeable, Autoregressive, Unstructured
)
import statsmodels.api as sm


def compare_working_correlations(df, formula, groups, family):
    """5 가지 작동 상관 결과 + naive/robust SE 비교"""
    structures = {
        "Independence": Independence(),
        "Exchangeable": Exchangeable(),
        "AR(1)": Autoregressive(),
        "Unstructured": Unstructured(),
    }

    rows = []
    for name, cs in structures.items():
        fit = GEE.from_formula(formula, groups=groups, data=df,
                                cov_struct=cs, family=family).fit()

        for param in fit.params.index:
            rows.append({
                "WorkingCorr": name,
                "Parameter": param,
                "Estimate": fit.params[param],
                "Naive SE": np.sqrt(fit.naive_covariance[
                    fit.model.exog_names.index(param)
                ][fit.model.exog_names.index(param)]),
                "Robust SE": fit.bse[param],
                "Ratio": fit.bse[param] / np.sqrt(fit.naive_covariance[
                    fit.model.exog_names.index(param)
                ][fit.model.exog_names.index(param)]),
            })

    return pd.DataFrame(rows)

실용 진단

Ratio = Robust / Naive 가:

약 1.0: 작동 상관 적절.
1.2: 작동 상관 너무 단순 — 다른 형태 시도.
< 0.8: 과적합 가능 (드뭄).

여러 작동 상관에 대해 Estimate 가 비슷 + Robust SE 도 비슷하면 robust 성질 확인.

8 핵심 정리

한 페이지 요약

§ 8.3 — GEE 모형:

4 명세: GLM 의 3 가지 + 작동 상관 행렬 \(R_i(a)\).
Robust 성질: 작동 상관 misspecification 에도 회귀 계수 일치 추정. 효율성 5~20% 손실.
5 작동 상관: Independence (0 모수), Exchangeable (1), AR(1) (1), m-dependent (\(m\)), Unstructured (\(n(n-1)/2\)).
CPM 와 매핑: GEE 의 5 작동 상관 ≈ CPM 의 5 구조 (Independence ↔︎ \(\sigma^2 I\) 추가).

§ 8.4 — GEE 추정:

작동 분산-공분산 (식 8.15): \(V_i = \phi A_i^{1/2} R_i A_i^{1/2}\). 분산 함수 + 작동 상관.
추정 방정식 (식 8.17): GLM 의 (8.10) 의 직접 일반화.
IRLS 두 단계: (a) WLS 로 \(\hat\beta\) 갱신, (b) Pearson 잔차로 \(\hat a, \hat\phi\) 갱신. 수렴까지 반복.
Pearson 잔차 (식 8.20): 분포의 marginal 분산 제거 — 시점 간 비교 가능 단위.
Naive SE (식 8.21): \(V_{\text{naive}} = M_0^{-1}\) (작동 상관 정확 가정).
Sandwich SE (식 8.22): \(V_{\text{robust}} = M_0^{-1} M_1 M_0^{-1}\) (Royall 1986). 작동 상관 misspecification 에도 일치 추정.
빵-속재료 직관: \(M_0^{-1}\) 가 빵, \(M_1\) 이 잔차 기반 empirical 분산.
Naive vs Robust 진단: 두 SE 의 비율로 작동 상관 misspecification 탐지.
작은 표본 함정: Sandwich SE 가 \(N < 30\) 에서 하향 편향 → MD/KC/FG 보정 또는 bootstrap.

§ 8.3 + § 8.4 가 GEE 의 작동 메커니즘 — 모형 정의 (5 작동 상관) + 추정 (IRLS) + SE (Sandwich). § 8.5 Gruder 사례 가 이 절차의 임상 적용.

9 다음 단계

주제	내용	위치
§ 8.5 Gruder 흡연 절제 사례	분석 정밀 재현 + Helmert contrasts 임상 해석	작성 예정 (`08-3-mrm-gee-gruder.qmd`)
Ch.9 GLMM 이항	population-averaged vs subject-specific 비교	미작성 (mm-06 참조)

10 관련 주제

선행 지식

Ch.8 Overview — GEE — 4 절 systematic
§ 8.1-8.2 — GEE 의 토대 — Marginal 정의 + GLM 토대
Ch.6 Overview — CPM — 5 분산-공분산 구조 (GEE 와 매핑)