GLMM Parameter Estimation — 준-우도 방정식과 이차 형식 교대 추정 (McCullagh §14.4)

1 서론 — 왜 완전 우도 대신 준-우도인가

§14.3 에서 GLMM 의 구조 — 조건부 GLM + 정규 랜덤효과 — 를 세웠다. 추정 은 어떻게 하는가?

가장 당연한 답: 완전 최대 가능도 (Full MLE). 랜덤효과 \(\gamma\) 를 적분하여 주변 우도를 구성하고 최대화.

\[ L(\beta, \sigma^2, \sigma_b^2) = \prod_j \int f(y_{\cdot j} \mid \gamma_j) \cdot \phi(\gamma_j; 0, \sigma_b^2) \, d\gamma_j. \]

문제: 이 적분은 닫힌 해가 거의 없다. Gauss-Hermite 구적 · Laplace 근사 · MCMC 같은 수치 방법이 필요. 1989 년 당시에는 매우 부담스러웠고, 2026 년에도 고차원·비선형 문제에서 여전히 비싸다.

McCullagh-Nelder 의 철학적 선택: “확신할 수 있는 최소 가정” 으로 시작. 즉 완전 분포 가 아니라 주변 평균 · 공분산 행렬 만 명세하고 그 위에서 준-우도 점수 방정식을 쓴다.

이번 글은 §14.4 의 세 가지 핵심을 심화한다.

1. \(\beta\) 의 추정 방정식 (14.6) — Ch.9 의 준-우도 구조 재활용.
2. 공분산 \(V(\mu, \sigma^2)\) 의 가법 분해 (14.7) — 각 분산 성분을 식별 가능 모집단에 대응.
3. 이차 형식 \(Q_r\) 과 그 기대값 (14.8-9) — 분산 성분 추정을 선형 연립방정식 으로 환원.

2 준-우도 점수 방정식 (14.6)

2.1 설정

주변 모형의 평균과 공분산:

\[ E(Y) = \mu(\beta), \qquad \text{cov}(Y) = V(\mu, \sigma^2). \]

여기서 \(\sigma^2 = (\sigma_1^2, \ldots, \sigma_k^2)\) 는 분산 성분 벡터 (Ch.14 에서 \(\sigma^2, \sigma_b^2\) 등 여러 성분).

2.2 점수 방정식

Ch.9 의 준-우도 점수 방정식 (§9.2 식 9.5) 을 그대로 적용:

\[ \boxed{\;U(\widehat\beta, \sigma^2) = D^T V^{-1} (y - \mu(\beta)) = 0\;} \tag{14.6} \]

여기서 \(D = \partial \mu / \partial \beta^T\) 는 평균의 \(\beta\) 에 대한 Jacobian.

2.3 직관

이 식의 각 요소 해석:

• \((y - \mu)\): 잔차 — 관측과 모형의 차이.
• \(V^{-1}\): 정보 가중 — 분산이 작은 관측치에 더 큰 가중치.
• \(D^T\): \(\beta\) 에 대한 민감도 — 각 모수가 어떤 평균 성분에 영향을 주는지.

요약하면: “\(\beta\) 의 각 방향으로 가중 잔차 합이 0 이 되도록 조정”.

이것이 표준 GLM 의 IRLS 와 구조적으로 동일. 차이는 \(V\) 가 대각이 아니라 cluster 내 공분산 을 담은 밀도 행렬이라는 점.

2.4 \(\sigma^2\) 에 대한 의존성

식 (14.6) 의 해는 \(\sigma^2\) 에 의존. 그러나 McCullagh-Nelder 의 관찰:

“Ordinarily, but with some important exceptions, the solution… depends on ratios of dispersion components.”

즉 \(\widehat\beta\) 가 \(\sigma^2, \sigma_b^2\) 의 절대값이 아니라 비율 에 의존하는 경우가 많다. 이 성질이 추정의 구조적 안정성 을 낳는다 — 분산 성분이 다소 부정확해도 \(\widehat\beta\) 는 크게 흔들리지 않는다.

2.5 \(\widehat\beta\) 의 점근 공분산

\[ \text{cov}(\widehat\beta) \simeq (D^T V^{-1} D)^{-1}. \]

이 식이 \(\sigma^2\) 에 민감. 분산 성분이 편향되면 SE 도 편향된다. 점 추정은 견고하지만 구간 추정 · 가설 검정 에서는 \(\sigma^2\) 추정의 정확도가 중요.

두 층의 정밀도 차이

양	\(\sigma^2\) 정확도에 대한 민감도
\(\widehat\beta\) 점 추정	낮음 — 비율만 대략 맞으면 됨
\(\text{SE}(\widehat\beta)\)	높음 — 절대값 정확히 필요
가설 검정 p-value	매우 높음 — SE × t 통계량 조합

실무 함의: 분산 성분 추정의 정확도 가 GLMM 분석의 병목. 이 때문에 §14.4 가 \(\sigma^2\) 추정을 위한 별도 방정식 체계를 도입한다.

3 공분산의 가법 분해 (14.7)

3.1 일반 형태

McCullagh-Nelder 의 관찰: 실무에서 등장하는 대부분의 GLMM 에서 공분산 행렬이 가법적으로 분해 된다.

\[ V(\mu, \sigma^2) = \sigma_1^2 V_1(\mu) + \sigma_2^2 V_2(\mu) + \cdots + \sigma_k^2 V_k(\mu). \tag{14.7} \]

각 \(V_j(\mu)\) 는 고정 구조 행렬 (평균에만 의존), \(\sigma_j^2\) 는 스칼라 크기 모수.

3.2 Rank 의 의미

McCullagh-Nelder 의 관찰: “각 \(V_j\) 의 rank 는 모집단 \(j\) 에서 추출된 표본 원소 수와 같다.”

이 뜻은: - \(V_1\) (within-cluster, 오차 분산) — rank = \(n\) (모든 관측치가 오차 변동 받음). - \(V_2\) (between-cluster) — rank = cluster 수 \(J\) (각 cluster 에만 독립 변동). - \(V_3\) (더 깊은 계층) — rank 가 더 작음.

rank 계층 이 자유도 계층 과 대응 — 각 분산 성분이 식별되는 관측 차원이 다르다.

3.3 왜 가법 분해가 중요한가

계산 · 해석 양면에서 이점:

계산: - 각 \(V_j\) 를 별도로 구성 가능. 코드 모듈화 쉬움. - \(\sigma_j^2\) 갱신 시 \(V_j\) 곱만 재계산. - 대규모 sparse 구조에서 메모리 절약.

해석: - 각 \(\sigma_j^2\) 이 물리적 모집단 에 대응. “between-school variance” 같은 레이블링. - 분산의 출처 를 분리해 보고 가능.

4 §14.3 예제 — 분해의 구체화

4.1 주변 공분산 유도

§14.3 (13-3 포스트) 의 주변 분산 공식 (14.5) 재활용:

\[ \text{var}(Y_{ij(k)}) = \sigma^2 E[V(M_{ij(k)})] + \text{var}(M_{ij(k)}). \]

\(V(M) = M\) (Poisson 조건부) 이면:

\[ \text{var}(Y_{ij(k)}) = \sigma^2 \mu_{ij(k)} + \{e^{\sigma_b^2} - 1\} \mu_{ij(k)}^2. \]

4.2 공분산 (교차 관측)

같은 cow class 내 의 서로 다른 관측치 (다른 site 또는 treatment) 간 공분산은 0 이 아님:

\[ \text{cov}(Y_{ij(k)}, Y_{i'j(k')}) = \text{cov}(M_{ij(k)}, M_{i'j(k')}). \]

조건부로는 독립이지만, \(\gamma_j\) 를 공유하므로 주변으로는 양의 공분산.

계산:

\[ M_{ij(k)} \cdot M_{i'j(k')} = e^{\alpha_i + \gamma_j + \tau_k} \cdot e^{\alpha_{i'} + \gamma_j + \tau_{k'}} = e^{\alpha_i + \alpha_{i'} + 2\gamma_j + \tau_k + \tau_{k'}}. \]

\(E[e^{2\gamma_j}] = e^{2\sigma_b^2}\) (로그-정규 2차 적률). 그리고 \(E[M_{ij(k)}] E[M_{i'j(k')}] = e^{\alpha_i + \alpha_{i'} + \tau_k + \tau_{k'}} \cdot e^{\sigma_b^2}\). 차이:

\[ \text{cov}(M, M') = e^{\alpha_i + \alpha_{i'} + \tau_k + \tau_{k'}} \cdot (e^{2\sigma_b^2} - e^{\sigma_b^2}) = \mu_{ij(k)} \mu_{i'j(k')} \cdot (e^{\sigma_b^2} - 1) / e^{\sigma_b^2} \cdot e^{\sigma_b^2} = \mu_{ij(k)} \mu_{i'j(k')} (e^{\sigma_b^2} - 1). \]

(위 계산은 \(\mu = e^{\sigma_b^2/2} \cdot e^{\alpha + \tau}\) 관계를 사용.)

다른 cow class 간: \(\gamma_j, \gamma_{j'}\) 가 독립이므로 공분산 0.

4.3 행렬 형태 — \(V_1, V_2\) 의 식별

전체 공분산 행렬을 보면:

대각: \(\sigma^2 \mu + (e^{\sigma_b^2} - 1) \mu^2\).

같은 cow class 내 교차: \((e^{\sigma_b^2} - 1) \mu \mu'\).

다른 cow class: 0.

이를 (14.7) 형태로 쓰면:

\[ V = \underbrace{\sigma^2 \cdot \text{diag}(\mu)}_{V_1 \text{ 곱 } \sigma^2_1} + \underbrace{(e^{\sigma_b^2} - 1) \cdot \mu J \mu^T}_{V_2 \text{ 곱 } \sigma_2^2}. \]

여기서:

• \(\sigma_1^2 = \sigma^2\), \(V_1(\mu) = \text{diag}(\mu)\) — 대각 행렬, within-class 변동.
• \(\sigma_2^2 = e^{\sigma_b^2} - 1\), \(V_2(\mu) = \mu J \mu^T\) — rank-low 행렬, between-class 변동.
  • \(J\) 는 블록 대각 단위 행렬: \(J_{rs} = 1\) if \(r, s\) 가 같은 cow class, 0 otherwise.
  • \(\mu J \mu^T\) 는 요소별로 \(\mu_r \cdot J_{rs} \cdot \mu_s\).

4.4 \(J\) 의 기하적 의미

\(J\) 는 “같은 cluster 에 속하는가” 의 지시 행렬. \(n \times n\) 이지만 rank = 클러스터 수 \(J\) (각 클러스터가 길이 \(n_j\) 의 all-ones 벡터를 생성).

\(\mu J \mu^T\) 는 각 클러스터 내에서 \(\mu\) 벡터의 outer product 를 합한 것. 이 행렬의 rank 역시 클러스터 수.

직관: \(V_2 = \mu J \mu^T\) 의 해석

\(\gamma_j\) 가 cow class \(j\) 의 “멀티플라이어” 역할을 하는 랜덤효과. 같은 class 의 모든 관측치가 같은 \(e^{\gamma_j}\) 로 scaling 된다.

따라서 관측치 \((r, s)\) 가 같은 class 이면 둘 다 같은 \(e^{\gamma_j}\) 배율로 흔들림 → 양의 공분산. 공분산 크기는 두 기대값의 곱 \(\mu_r \mu_s\) 에 비례 — 큰 관측치일수록 더 많이 흔들리기 때문.

\(\mu J \mu^T\) 의 \((r, s)\) 원소가 \(\mu_r J_{rs} \mu_s\) 로 이 구조를 정확히 담는다.

5 이차 형식 \(Q_r\) 과 기대값 공식 (14.8-9)

5.1 동기 — \(\sigma^2\) 를 어떻게 추정할 것인가

\(\beta\) 추정에는 (14.6) 이 있지만, \(\sigma^2_j\) 성분 추정은 별도 방정식이 필요. McCullagh-Nelder 의 해법: 이차 형식의 기대값 등치 로 MoM 일반화.

5.2 이차 형식 정의

\(k\) 개 적절한 사영 행렬 \(P_r\) (\(r = 1, \ldots, k\)) 을 잡고

\[ Q_r = (Y - \mu)^T P_r (Y - \mu). \tag{14.8} \]

\(P_r\) 선택: \(r\) 번째 랜덤효과에 대응하는 주변 합계 사영 이 자연스러움. 예:

• \(P_1\) = 개별 관측치 사영 = \(I\) (또는 $I - $ (고정효과 공간 사영)).
• \(P_2\) = cow class 평균 사영 (각 class 의 관측치를 평균내어 4 또는 8 값으로 축소).

일반적으로 \(P_r\) 은 대칭 · 멱등 (사영) 이어야 한다.

5.3 기대값 공식

\(\text{cov}(Y) = V = \sum_j \sigma_j^2 V_j\) 를 대입:

\[ E(Q_r) = E[(Y-\mu)^T P_r (Y-\mu)] = \text{tr}(P_r V) = \sum_{j=1}^{k} \text{tr}(P_r V_j) \cdot \sigma_j^2. \tag{14.9} \]

결과: \(E(Q_r)\) 이 \(\sigma_1^2, \ldots, \sigma_k^2\) 의 선형 함수.

5.4 MoM 연립방정식

\(k\) 개 이차 형식에 대해 \(k\) 개 방정식:

\[ \begin{pmatrix} Q_1 \\ Q_2 \\ \vdots \\ Q_k \end{pmatrix} = \begin{pmatrix} \text{tr}(P_1 V_1) & \cdots & \text{tr}(P_1 V_k) \\ \vdots & & \vdots \\ \text{tr}(P_k V_1) & \cdots & \text{tr}(P_k V_k) \end{pmatrix} \begin{pmatrix} \sigma_1^2 \\ \vdots \\ \sigma_k^2 \end{pmatrix}. \]

\(k \times k\) 선형계 — \(\sigma^2\) 벡터를 직접 풀 수 있다.

5.5 음수 추정치 가능성

“As usual with variance components, the estimates obtained need not be positive.”

선형계의 해가 음수 가 나올 수 있다. 물리적으로는 불가능 (\(\sigma^2 \geq 0\)) 하지만, 표본 변동 때문에 발생.

실무 대응: 1. 0 으로 절단 (truncation): 음수 추정치를 0 으로. 편향 도입. 2. 제약 최적화 (constrained optimization): \(\sigma^2 \geq 0\) 제약 하 우도 최대화. 3. REML: 이 문제에서 더 효율적 (아래).

5.6 구체적 계산 — §14.3 예제

\(V_1 = \text{diag}(\mu)\), \(V_2 = \mu J \mu^T\) 에 대해:

\(P_1 = I - H\) (고정효과 공간에 직교 사영, \(H\) 는 햇 행렬):

\[\text{tr}(P_1 V_1) = \text{tr}((I-H)\text{diag}(\mu)) = \sum (1 - h_i) \mu_i,\]

\[\text{tr}(P_1 V_2) = \text{tr}((I-H) \mu J \mu^T) = \text{지역 대칭 항}.\]

\(P_2\) = cow class 평균 구성 사영:

두 행렬 \(V_1, V_2\) 에 대한 대각합 계산. 구체적 공식은 설계의 세부에 의존.

실무 팁 — ANOVA 분해와의 연결

§14.2 의 ANOVA 표 (Table 14.1) 는 이 공식의 구체적 응용. 각 “Source” 줄이 특정 \(P_r\) 에 대응:

• Error → \(P_1 = I - H\): \(E(MS_E) = \sigma^2\) (balanced 정규 경우).
• Cow class → $P_2 = $ between-class 사영: \(E(MS_{\text{cow}}) = \sigma^2 + n_j \sigma_b^2\).

선형 혼합 모형의 MoM 추정이 ANOVA 표의 일반화. Table 14.1 의 \(E(MS)\) 열이 (14.9) 의 특수 사례다.

6 교대 반복 알고리즘

6.1 상호 의존성

(14.6) 과 (14.8-9) 는 서로 의존:

• (14.6): \(\widehat\beta\) 추정에 \(V\) 가 필요, \(V\) 는 \(\sigma^2\) 에 의존.
• (14.8): \(Q_r\) 계산에 \(\widehat\mu\) 가 필요, \(\widehat\mu = \mu(\widehat\beta)\).

6.2 알고리즘

McCullagh-Nelder 의 실무 방법:

초기화: σ² = 초기값 (예: 순수 GLM 잔차로부터 단일 분산)
1. 주어진 σ² 로 V 구성 → (14.6) 풀어 β̂ 얻음
2. β̂ 에서 μ̂ 계산
3. μ̂ 를 (14.8) 에 대입해 Q_r 계산
4. (14.9) 연립방정식 풀어 σ̂² 갱신
5. 수렴 확인; 미수렴이면 1 로

수렴 성질: 실무에서 보통 5-10 반복 에 수렴. 각 단계의 로그 우도 감소가 보장되는 것은 아니지만 (준-우도이므로) 실제로는 잘 작동.

6.3 Ch.13 Weibull 과의 유비

§13.3 Weibull 의 \((\beta, \alpha)\) 교대 추정과 구조적으로 같다: - Ch.13: offset 에 \(\alpha \log t\) → \(\beta\) 와 \(\alpha\) 교대. - Ch.14: 분산 \(V = V(\sigma^2)\) → \(\beta\) 와 \(\sigma^2\) 교대.

두 경우 모두 모수 공간을 분할 하고 각 부분공간에서 부분 최대화를 반복. Gauss-Seidel 스타일 block coordinate descent 의 통계 응용.

7 REML 대안

7.1 간략 소개

Restricted Maximum Likelihood (REML) (§7.2 참조) 은 분산 성분 추정의 표준 대안.

차이: - MoM (14.8-9): 이차 형식의 기대값 등치. 간단, 선형계. - REML: 고정효과를 적분 제거한 뒤 분산 성분 우도 최대화. 더 효율적이지만 \(V(\mu, \sigma^2)\) 의 반복 역행렬 계산 필요.

7.2 언제 REML 이 나은가

• 정규성 가정 이 성립할 때 점근 효율 더 높음.
• 불균형 설계 에서 MoM 편향이 큼 → REML 개선.
• 음수 추정치 문제에서 제약 더 자연스러움.

7.3 언제 MoM 이 나은가

• 계산 비용 제약.
• 정규성 가정 부적절 — REML 도 정규성 없이 쓸 수 있지만 효율 주장이 약해짐.
• GLMM 의 비정규 반응: 정규 아니면 REML 의 이론적 기반이 약해지므로 MoM 로 충분.

McCullagh-Nelder 의 선택: 간단한 MoM — “less efficient but simpler”.

8 Python 실전 — §14.3 예제 구체화

8.1 \(V\) 행렬 명시적 구성

import numpy as np
import pandas as pd

np.random.seed(42)
n_classes = 8
n_per_class = 4  # 각 class 에 4 사이트

# 참 모수
sigma2 = 0.25
sigma_b2 = 0.1
mu = np.random.uniform(2, 10, n_classes * n_per_class)  # 예시 평균

# V_1 = diag(μ), V_2 = μ J μ^T
V1 = np.diag(mu)

# J: 같은 cow class 인지 표시
cow_class = np.repeat(np.arange(n_classes), n_per_class)
J = (cow_class[:, None] == cow_class[None, :]).astype(float)

V2 = np.outer(mu, mu) * J

# 전체 공분산
V = sigma2 * V1 + (np.exp(sigma_b2) - 1) * V2

print(f"V 차원: {V.shape}")
print(f"rank(V_1) = {np.linalg.matrix_rank(V1)} (= n)")
print(f"rank(V_2) = {np.linalg.matrix_rank(V2)} (= 클러스터 수)")

8.2 \(Q_r\) 계산과 \(\sigma^2\) 복원

# 가상 데이터 생성 (V 로부터 다변량 정규)
y = np.random.multivariate_normal(mu, V)

# 사영 행렬들
# P_1 = I (단순화, 고정효과 공간 무시)
# P_2 = 클러스터 평균 - 전체 평균 의 사영
n = len(y)

# 클러스터별 평균 연산자
M = np.zeros((n_classes, n))
for j in range(n_classes):
    M[j, cow_class == j] = 1.0 / n_per_class

P1_proj = np.eye(n)  # 단순화: (I - H) 대신 I
# P_2: 클러스터 평균을 복제한 뒤 전체 평균을 뺌
ones_n = np.ones((n, n)) / n
cluster_mean_expanded = M.T @ M * n_per_class  # n x n projection
P2_proj = cluster_mean_expanded / n_per_class - ones_n

# Q_r 계산
r = y - mu
Q1 = r @ P1_proj @ r
Q2 = r @ P2_proj @ r

# 기대값 계수
t_P1V1 = np.trace(P1_proj @ V1)
t_P1V2 = np.trace(P1_proj @ V2)
t_P2V1 = np.trace(P2_proj @ V1)
t_P2V2 = np.trace(P2_proj @ V2)

print(f"\nE(Q_1) = {t_P1V1:.2f}·σ² + {t_P1V2:.2f}·σ_2²")
print(f"E(Q_2) = {t_P2V1:.4f}·σ² + {t_P2V2:.2f}·σ_2²")

# 관측 Q 와 기대값 등치
A = np.array([[t_P1V1, t_P1V2], [t_P2V1, t_P2V2]])
b = np.array([Q1, Q2])
sigma2_est = np.linalg.solve(A, b)

print(f"\n관측 Q_1 = {Q1:.2f}, Q_2 = {Q2:.4f}")
print(f"σ² 추정: {sigma2_est[0]:.3f} (참: {sigma2:.3f})")
print(f"σ_2² 추정 (=exp(σ_b²)-1): {sigma2_est[1]:.4f} (참: {np.exp(sigma_b2)-1:.4f})")
print(f"역산 σ_b² ≈ log(1+σ_2²): {np.log(1+sigma2_est[1]):.4f} (참: {sigma_b2:.4f})")

기대: 한 번의 시뮬레이션에서 단일 추정치가 참값 근처. 반복 시뮬로 표본 분포를 확인하면 분산 성분 추정의 변동을 체감 가능.

8.3 전체 교대 반복

def glmm_alternating(y, X, design_info, max_iter=10):
    """
    X: 고정효과 설계 행렬
    design_info: cluster 구조
    """
    n = len(y)

    # 초기화
    beta = np.linalg.lstsq(X, y, rcond=None)[0]
    sigma2 = np.var(y - X @ beta)
    sigma_b2 = 0.01  # 작은 초기값

    for it in range(max_iter):
        # 1. 현재 σ² 로 V 구성
        mu_hat = np.exp(X @ beta)  # log link 가정
        V1 = np.diag(mu_hat)
        V2 = np.outer(mu_hat, mu_hat) * J
        V = sigma2 * V1 + (np.exp(sigma_b2) - 1) * V2

        # 2. β̂ 갱신 (14.6)
        # D = ∂μ/∂β^T = diag(μ) X (로그 링크)
        D = mu_hat[:, None] * X
        beta_new = beta + np.linalg.solve(D.T @ np.linalg.solve(V, D),
                                           D.T @ np.linalg.solve(V, y - mu_hat))

        # 3-4. σ² 갱신 (14.8-9)
        mu_hat = np.exp(X @ beta_new)
        V1 = np.diag(mu_hat)
        V2 = np.outer(mu_hat, mu_hat) * J
        r = y - mu_hat

        Q1 = r @ P1_proj @ r
        Q2 = r @ P2_proj @ r

        A = np.array([[np.trace(P1_proj @ V1), np.trace(P1_proj @ V2)],
                      [np.trace(P2_proj @ V1), np.trace(P2_proj @ V2)]])
        sigma_new = np.linalg.solve(A, [Q1, Q2])

        # σ² ≥ 0 으로 절단
        sigma2 = max(sigma_new[0], 1e-6)
        sigma_b2 = max(np.log(1 + sigma_new[1]), 1e-6)

        # 수렴 확인
        change = np.max(np.abs(beta - beta_new))
        beta = beta_new
        print(f"iter {it}: β change = {change:.5f}, σ² = {sigma2:.3f}, σ_b² = {sigma_b2:.3f}")
        if change < 1e-4:
            break

    return beta, sigma2, sigma_b2

이 코드는 개념 구현. 실제 현장에서는 statsmodels.MixedLM 또는 R lme4::glmer 사용.

9 요약 — §14.4 의 네 가지 교훈

9.1 교훈 1 — 준-우도 점수 방정식의 재활용

(14.6) 은 Ch.9 에서 이미 본 구조. GLMM 으로 확장할 때 추가 복잡도는 \(V\) 가 클러스터 공분산 을 담는다는 점뿐. 이것이 McCullagh-Nelder 의 철학 — “새 문제를 기존 도구로 공격”.

9.2 교훈 2 — 가법 분해가 계산 · 해석을 명확화

(14.7) 의 \(V = \sum \sigma_j^2 V_j\) 구조는: - 각 \(\sigma_j^2\) 을 물리적 모집단에 대응시켜 해석 용이. - 고정 \(V_j\) 행렬 로 계산 모듈화. - Rank 계층 으로 식별 가능성 점검.

9.3 교훈 3 — MoM 의 이차 형식 구조

\(Q_r\) 이 “사영 행렬 × 잔차 내적” 구조라는 점이 중요. ANOVA 표의 SS 줄 이 이것의 구체적 사례. “분산 분석 → 혼합 모형 → GLMM” 으로 이어지는 추정 철학의 연속성.

9.4 교훈 4 — 교대 반복 = 모수 공간 분할

\(\beta\) 와 \(\sigma^2\) 를 분리해서 추정하는 것이 완전 동시 최대화보다 계산 안정. Weibull (§13.3), 음이항 (§11.2), Box-Tidwell (§11.4) 등 McCullagh-Nelder 의 여러 교대 알고리즘의 공통 패턴.

9.5 한 줄 정리

“준-우도 + 가법 분해 + 이차 형식 MoM + 교대 반복” 네 요소가 GLMM 추정 파이프라인의 McCullagh-Nelder 판본 을 구성한다. 이후 30 년간 PQL · Laplace · MCMC 가 이 틀 위에 정교한 대안을 쌓아 올렸다.

10 관련 주제

선행 지식

• Components of Dispersion — 개관 (McCullagh Ch.14)
• Linear Mixed Models — 라틴 정방 (McCullagh §14.2)
• Non-Linear Mixed Models (GLMM) — Table 14.2 aliasing (McCullagh §14.3)
• Quasi-likelihood Functions (McCullagh Ch.9) — 준-우도 점수 방정식의 기원
• Optimal Estimating Functions (McCullagh §9.4)

관련 개념

• Non-Linear Parameters in the Variance Function (McCullagh §11.2) — 교대 추정의 다른 사례
• Non-Linear Parameters in the Covariates (McCullagh §11.4) — Box-Tidwell 교대
• Parametric Survival — Weibull 교대 (McCullagh §13.3)
• Extended Quasi-Likelihood (McCullagh §9.6)

실무 참고 — 현대 GLMM 구현

• R lme4::glmer — Laplace/AGQ 기반
• R nlme::lme — REML 기반 LMM
• Python statsmodels.MixedLM — REML 기반 LMM
• Python pymc, Stan — 베이지안 GLMM
• SAS PROC GLIMMIX, PROC NLMIXED

후속 주제

• Salamander Mating Experiment (McCullagh §14.5) — 이 알고리즘의 실전 적용
• Ch.14 Exercises (McCullagh §14.7)