1 들어가며 — 왜 이항 GLMM 추정이 어려운가

Ch.9 Overview → § 9.4 → § 9.5.1-9.5.5 에서 mixed-effects logistic 모형의 정의를 다뤘다. 본 포스트는 그 모형의 모수를 데이터에서 추정하는 § 9.6 전체 — 본문 + 세 소절 — 를 정리한다.

정규 MRM 추정과 결정적으로 다른 점

§ 4.5 의 정규 MRM 추정 과 비교:

항목	정규 MRM	이항 GLMM
우도 형태	다변량 정규 (closed form)	적분 안에 비선형 함수 (no closed form)
랜덤 효과 적분	해석적으로 가능	수치 근사 필요 (Gauss-Hermite, Laplace, MCMC)
EM 알고리즘	E-step·M-step 모두 닫힌 형태	E-step 도 적분 필요
추정 방법	ML/REML, 행렬 연산만	Fisher scoring + 수치 적분 결합
표준오차	잠근 형태 inverse	Information matrix 의 적분 추정

한 줄 요약

“이항 GLMM 의 우도는 랜덤 효과에 대한 적분 안에 비선형 logistic cdf 가 들어 있어 닫힌 형태로 풀 수 없다 (식 9.40). 그래서 (1) 적분을 Gauss-Hermite quadrature 로 수치 근사하고 (식 9.55-9.59), (2) 이 근사 우도에 Fisher scoring 을 적용해 모수를 update 하며 (식 9.47), (3) 환자별 랜덤 효과는 사후 평균 (empirical Bayes) 으로 별도 추정한다 (식 9.48). 다차원 랜덤 효과에서는 적분 차원이 \(Q^r\) 로 폭발하므로 adaptive quadrature 로 점 수를 줄인다 (식 9.60-9.61).”

2 § 9.6 본문 — Marginal Maximum Likelihood

2.1 기본 설정 — 랜덤 절편 모형 (식 9.36)

출발점 모형

피험자 \(i = 1, \ldots, N\), 각 피험자의 반복 측정 \(j = 1, \ldots, n_i\). 응답 \(Y_{ij} \in \{0, 1\}\), \(\theta_i \sim \mathcal{N}(0, 1)\) (표준화 랜덤 효과).

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = x_{ij}^\top \beta + \sigma_v \theta_i \tag{9.36} \]

추정 대상 모수: \(\Theta = (\beta^\top, \sigma_v)^\top\) — 회귀 계수 벡터와 랜덤 효과 표준편차.

추정 절차의 큰 그림 — 두 단계의 결합

이항 GLMM 추정은 두 단계의 추정 사이를 진동한다:

모수 추정 (\(\beta, \sigma_v\)): 모집단 수준 회귀 계수와 분산 — 모든 환자가 공유.
랜덤 효과 추정 (\(\theta_i\)): 환자별 잠재 효과 — 환자마다 다름.

전자는 marginal maximum likelihood (랜덤 효과를 적분으로 제거) 로, 후자는 empirical Bayes (사후 평균) 로 추정. 두 추정이 분리되는 이유 — 전자는 모집단 추론이라 표본 정보를 모두 합쳐야 하고, 후자는 개인 추론이라 환자 데이터에 더 의존한다.

2.2 Conditional Probability — 식 (9.37-9.38)

조건부 확률

\(z_{ij} = x_{ij}^\top \beta + \sigma_v \theta_i\) 로 두면:

\[ P(Y_{ij} = 1 \mid \theta_i) = \Psi(z_{ij}) \tag{9.37} \]

\[ P(Y_{ij} = 0 \mid \theta_i) = 1 - \Psi(z_{ij}) \tag{9.38} \]

여기서 \(\Psi\) 는 logistic cdf (식 9.33).

직관 — “환자 효과를 알면 시점이 독립”

\(\theta_i\) 가 알려진 값이면 같은 환자의 다른 시점 \(j, j'\) 의 응답은 독립이라고 본다. 이유는 환자의 모든 시점에 공통으로 작용하는 잠재 변수가 \(\theta_i\) 한 개뿐이라 가정했기 때문 — 이 효과를 빼면 남는 것은 시점별 독립 잡음 \(\epsilon_{ij}\) 뿐.

이 가정은 모형 설계에서 매우 중요한 결정이다. 만약 같은 환자의 시점 사이에 (랜덤 효과로 설명되지 않는) 추가 상관이 있으면 이 가정이 깨진다 — 예: 측정 오차에 자기상관이 있는 경우. § 7 의 mixed regression with autocorrelated errors (MRM-AC) 가 이 가정을 완화하는 모형.

2.3 Conditional Independence — 핵심 가정

Conditional Independence Assumption

같은 환자 안의 모든 시점이 \(\theta_i\) 를 조건부로 한 상태에서 독립:

\[ P(Y_{i1}, Y_{i2}, \ldots, Y_{in_i} \mid \theta_i) = \prod_{j=1}^{n_i} P(Y_{ij} \mid \theta_i) \]

이 가정은 mixed-effects 모형의 토대 — 이것 없이는 식 (9.39) 의 곱 형태가 불가능.

직관 — Psychometric literature 의 핵심 개념

심리 측정학 (psychometric literature) 에서 conditional independence 는 잠재 변수 모형의 정의 자체. 잠재 변수 \(\theta_i\) (예: 응시자의 능력) 를 알면 같은 응시자의 모든 문항 응답이 독립이 된다 — 이것이 잠재 변수가 “데이터 안의 모든 종속성을 흡수한다” 는 뜻.

이 발상이 IRT (Item Response Theory) 의 토대이고 (§ 9.5.3 의 2PL/Rasch 와 연결), mixed-effects logistic 도 같은 가정 위에 서 있다.

가정이 깨지는 경우의 결과:

표준오차 과소 추정 (실제 정보보다 더 많은 정보가 있다고 착각).
p-value 가 부풀려지고 type I error 가 증가.
진단으로는 잔차 자기상관, 시간 lag 별 ICC 변화 등을 검사.

2.4 Conditional Likelihood — 식 (9.39)

식 (9.39) — 환자별 조건부 가능도

같은 환자의 \(n_i\) 개 시점 가능도를 곱하면:

\[ \ell(Y_i \mid \theta) = \prod_{j=1}^{n_i} \Psi(z_{ij})^{Y_{ij}} [1 - \Psi(z_{ij})]^{1 - Y_{ij}} \tag{9.39} \]

직관 — Cross-sectional 우도와의 평행

\(Y_{ij} = 1\) 일 때 항이 \(\Psi\), \(Y_{ij} = 0\) 일 때 항이 \(1 - \Psi\). 이는 단일 수준 로지스틱의 Bernoulli 우도 (식 9.5) 와 정확히 같은 형태.

차이는:

단일 수준 (cross-sectional): 각 환자가 한 번만 관측, 환자별 1 개 항.
Mixed-effects (longitudinal): 각 환자가 \(n_i\) 번 관측, 환자별 \(n_i\) 항의 곱.

이 곱이 가능한 이유가 conditional independence — \(\theta_i\) 가 고정되면 시점이 독립.

식 (9.39) 는 아직 conditional 우도. \(\theta_i\) 의 값이 모르므로 그대로는 추정에 못 쓴다. 다음 단계에서 이 \(\theta_i\) 를 적분으로 제거한다.

2.5 Marginal Likelihood — 식 (9.40-9.42)

식 (9.40) — 한 환자의 marginal 가능도

\(\theta_i\) 의 분포 \(g(\theta) = \mathcal{N}(0, 1)\) (또는 \(\mathcal{N}(0, \sigma_v^2)\), § 9.5.4 표기) 로 적분:

\[ h(Y_i) = \int_\theta \ell(Y_i \mid \theta) \, g(\theta) \, d\theta \tag{9.40} \]

식 (9.41-9.42) — 표본 전체의 marginal 가능도

환자가 서로 독립이라 가정 (각 환자가 자기 \(\theta_i\) 를 가짐):

\[ L = \prod_{i=1}^N h(Y_i) \tag{9.41} \]

로그 형태:

\[ \log L = \sum_{i=1}^N \log h(Y_i) \tag{9.42} \]

직관 — “가중 평균 확률” 로 본 식 (9.40)

식 (9.40) 의 적분은 다음과 같이 해석할 수 있다:

\(\theta\) 의 모든 가능한 값에 대해 conditional 우도 \(\ell(Y_i \mid \theta)\) 를 계산.
각 \(\theta\) 값의 확률 \(g(\theta)\) 를 가중치로 평균.
결과: \(\theta\) 가 무엇인지 모를 때의 marginal (평균) 우도.

비유 — 환자의 “잠재 성향” 을 모를 때, 모든 가능한 성향 (\(\theta = -3, -2, \ldots, +3\)) 에 대해 데이터의 가능성을 따져보고, 각 성향의 사전 확률 (정규 분포) 로 가중 평균.

이 적분이 닫힌 형태로 안 풀리는 이유 — \(\Psi(z) = 1/(1 + e^{-z})\) 가 \(\theta\) 에 비선형이고, \(g(\theta) = \exp(-\theta^2/2)/\sqrt{2\pi}\) 와의 곱을 적분할 때 표준 함수로 표현 못 함. 정규 MRM 에서는 우도와 사전이 모두 정규라 곱이 정규로 단힘 — 이항 GLMM 에는 없는 행운.

2.6 First Derivatives — 식 (9.43-9.46)

식 (9.43) — Score 의 일반 형태

\(\eta\) 를 \(\beta\) 또는 \(\sigma_v\) 로 두면:

\[ \frac{\partial \log L}{\partial \eta} = \sum_{i=1}^N h^{-1}(Y_i) \frac{\partial h(Y_i)}{\partial \eta} \tag{9.43} \]

식 (9.44) — \(h(Y_i)\) 의 미분 후 정리

식 (9.40) 에 chain rule 적용 + logistic cdf-pdf 관계 (\(\Psi'(z) = \Psi(z)[1-\Psi(z)]\)) 활용:

\[ \frac{\partial \log L}{\partial \eta} = \sum_{i=1}^N h^{-1}(Y_i) \int_\theta \sum_{j=1}^{n_i} \frac{Y_{ij} - \Psi(z_{ij})}{\Psi(z_{ij})[1 - \Psi(z_{ij})]} \cdot \partial \Psi(z_{ij}) \cdot \frac{\partial z_{ij}}{\partial \eta} \cdot \ell_i \cdot g(\theta) \, d\theta \tag{9.44} \]

여기서 \(\partial \Psi(z_{ij}) = \Psi(z_{ij})[1 - \Psi(z_{ij})]\) (logistic 의 cdf-pdf 관계). 이 항이 분모와 약분되어 다음 형태로 단순화:

\[ \propto \int_\theta \sum_{j=1}^{n_i} [Y_{ij} - \Psi(z_{ij})] \cdot \frac{\partial z_{ij}}{\partial \eta} \cdot \ell_i \cdot g(\theta) \, d\theta \]

\(\partial z_{ij} / \partial \eta\) 는 모수에 따라:

\[ \frac{\partial z_{ij}}{\partial \beta} = x_{ij}, \quad \frac{\partial z_{ij}}{\partial \sigma_v} = \theta_i \tag{9.45-9.46} \]

직관 — Score 의 의미는 “잔차 × 공변량의 가중 평균”

식 (9.44) 의 핵심 항 \(Y_{ij} - \Psi(z_{ij})\) 는 이항 잔차 — 관측치와 예측 확률의 차이. 단일 수준 로지스틱의 score 식 (9.7) 과 정확히 같은 구조.

차이는 가중치:

단일 수준: 잔차 × 공변량 → 직접 합.
Mixed-effects: 잔차 × 공변량 × conditional likelihood × 사전 분포 → \(\theta\) 에 대한 적분.

해석: “환자별 잠재 성향에 대한 가중 잔차의 평균”. \(\ell_i \cdot g(\theta) / h(Y_i)\) 가 사실상 사후 분포 \(f(\theta \mid Y_i)\) — 즉 score 가 사후 기댓값 의 형태. EM 알고리즘의 E-step 과 같은 형태가 자연 등장.

logistic cdf-pdf 관계 \(\Psi' = \Psi(1-\Psi)\) 가 식 (9.44) 의 분자·분모 약분을 만들어 score 형태를 단순화 — 이것이 § 9.5.5 에서 강조한 logistic 의 수학적 매력.

2.7 Probit vs Logit 의 계산 차이

Probit 모형의 score

식 (9.44) 에서 \(\Psi\) 를 \(\Phi\) (정규 cdf) 로, \(\partial \Psi\) 를 \(\phi\) (정규 pdf) 로 대체하면 probit 의 score. 그러나 분자·분모 약분이 일어나지 않아 더 복잡:

\[ \propto \int_\theta \sum_{j=1}^{n_i} \frac{Y_{ij} - \Phi(z_{ij})}{\Phi(z_{ij})[1 - \Phi(z_{ij})]} \cdot \phi(z_{ij}) \cdot \frac{\partial z_{ij}}{\partial \eta} \cdot \ell_i \cdot g(\theta) \, d\theta \]

\(\phi / [\Phi(1-\Phi)]\) 가 단순화되지 않음.

Logistic 이 추정에서 더 단순한 이유 (§ 9.5.5 의 활용)

§ 9.5.5 에서 강조한 \(\psi = \Psi(1 - \Psi)\) 가 추정에서 결정적 역할:

Score 단순화: 식 (9.44) 의 가중치 \(\partial \Psi / [\Psi(1-\Psi)] = 1\) 로 약분 → 잔차 × 공변량 형태.
계산 효율: \(\Psi(z) = 1/(1+e^{-z})\) 만 계산하면 pdf 도 즉시 도출.
수치 안정: erf 같은 외부 함수 호출 불필요.

probit 은 \(\phi\) 와 \(\Phi\) 가 별도 계산 필요 + 수치 적분 (erf) 필요. 알고리즘이 더 복잡하고 약간 더 느리다. 결과는 거의 같지만 (식 9.12 의 \(\beta_L \approx 1.81 \beta_P\) 척도 차이만 있음) 추정 비용이 다르다.

2.8 Newton-Raphson vs Fisher Scoring — 식 (9.47)

식 (9.47) — Fisher Scoring update

\(\Theta_\iota\) 를 iteration \(\iota\) 의 모수 벡터, \(I(\Theta_\iota)\) 를 정보 행렬 (information matrix) 로 두면:

\[ \Theta_{\iota + 1} = \Theta_\iota + I(\Theta_\iota)^{-1} \cdot \frac{\partial \log L}{\partial \Theta_\iota} \tag{9.47} \]

정보 행렬 (Bock & Aitkin 1981) — 음의 2계 도함수의 기댓값:

\[ I(\Theta_\iota) = E\left[ -\frac{\partial^2 \log L}{\partial \Theta_\iota \partial \Theta_\iota^\top} \right] = \sum_{i=1}^N h^{-2}(Y_i) \frac{\partial h(Y_i)}{\partial \Theta_\iota} \left( \frac{\partial h(Y_i)}{\partial \Theta_\iota} \right)^\top \]

직관 — Newton-Raphson 와 Fisher Scoring 의 차이

두 알고리즘 모두 “현재 점에서 곡률 정보로 다음 점을 결정” 하는 2차 방법. 차이는 어떤 곡률을 쓰는가:

Newton-Raphson: 관측 정보 행렬 \(-\partial^2 \log L / \partial \Theta \partial \Theta^\top\) — 현재 데이터에서 직접 계산한 Hessian.
Fisher Scoring: 기대 정보 행렬 \(E[-\partial^2 \log L / \partial \Theta \partial \Theta^\top]\) — Hessian 의 기댓값 (모집단 평균).

왜 Fisher Scoring 이 더 단순한가:

계산 단순성: 기대 Hessian 은 보통 score 의 outer product 합으로 표현 가능 (\(\sum_i s_i s_i^\top\) 형태). 2계 도함수 직접 계산 불필요.
양정치 보장: 기대 Hessian 은 항상 양정치 (정보 행렬의 정의). Newton-Raphson 의 관측 Hessian 은 멀리 있는 점에서 음정치가 될 수 있어 발산 위험.
표준오차 직접 산출: 수렴 시 \(I^{-1}\) 가 그대로 모수 추정의 분산-공분산 → Wald 검정·신뢰구간에 즉시 활용.

단점: 정확한 Hessian 보다 정보가 적어 수렴이 다소 느릴 수 있음. 실무에서는 두 방법을 혼합 — 초기에 Fisher Scoring 으로 안정 수렴, 끝부분에 Newton-Raphson 으로 정밀 마무리 (hybrid 알고리즘).

3 § 9.6.1 — Empirical Bayes 랜덤 효과 추정

3.1 사후 평균 — 식 (9.48-9.49)

식 (9.48) — Empirical Bayes (EAP) 추정

환자별 랜덤 효과의 사후 분포 평균:

\[ \widehat{\theta}_i = E(\theta_i \mid Y_i) = h^{-1}(Y_i) \int_\theta \theta_i \cdot \ell(Y_i \mid \theta) \cdot g(\theta) \, d\theta \tag{9.48} \]

EAP = Expected A Posteriori (Bock & Aitkin 1981 의 명명).

식 (9.49) — 사후 분산

\[ V(\widehat{\theta}_i \mid Y_i) = h^{-1}(Y_i) \int_\theta (\theta_i - \widehat{\theta}_i)^2 \cdot \ell(Y_i \mid Y_i) \cdot g(\theta) \, d\theta \tag{9.49} \]

직관 — Empirical Bayes 는 “데이터 + 사전” 의 절충

\(\theta_i\) 의 사후 분포를 Bayes 정리로 보면:

\[ f(\theta_i \mid Y_i) = \frac{\ell(Y_i \mid \theta_i) \, g(\theta_i)}{h(Y_i)} \]

사후 분포의 평균이 식 (9.48). 이 평균이 Empirical Bayes 추정값.

해석:

\(\widehat{\theta}_i\) 는 환자 \(i\) 의 데이터가 시사하는 잠재 효과와, 모집단 사전 분포 \(g(\theta)\) 사이의 절충.
환자 데이터가 풍부하면 (관측 \(n_i\) 가 큼) → 데이터 우도 \(\ell\) 가 우세 → \(\widehat{\theta}_i\) 가 환자 데이터에 수렴.
환자 데이터가 빈약하면 (관측 \(n_i\) 가 적음) → 사전 \(g\) 가 우세 → \(\widehat{\theta}_i\) 가 0 (사전 평균) 으로 shrink.

이 shrinkage (수축) 효과 가 mixed-effects 모형의 핵심 매력. 환자별로 자유 추정 (fixed effects) 보다 더 안정한 추정이 자동 — 데이터가 적은 환자가 모집단으로 끌어당겨진다.

“Empirical” Bayes 의 의미: 사전 분포의 모수 (\(\sigma_v\)) 를 데이터에서 추정 (식 9.40-9.42 의 marginal MLE) 하고, 그 추정값을 사용해 사후 평균을 계산. 완전한 Bayes (사전을 미리 지정) 와의 차이.

3.2 Subject-Specific vs Marginal Probabilities

두 종류의 추정 확률

Subject-specific 확률 \(\widehat{P}_{ss}\): 특정 환자의 랜덤 효과 \(\theta^*\) 에서의 응답 확률.

\[ \widehat{P}_{ss}(\theta^*) = \Psi(x^\top \widehat{\beta} + \widehat{\sigma}_v \theta^*) \]

Marginal 확률 \(\widehat{P}_m\): 환자 모집단 전체 평균 응답 확률.

\[ \widehat{P}_m = \int_\theta \widehat{P}_{ss}(\theta) \, g(\theta) \, d\theta \]

직관 — 같은 모형, 두 가지 답

질문이 다르다:

“평균 환자 (\(\theta = 0\)) 의 응답 확률은?” → subject-specific \(\widehat{P}_{ss}(0)\).
“전체 모집단의 평균 응답률은?” → marginal \(\widehat{P}_m\).

비선형 link 때문에 두 답이 다르다 — § 9.4 의 핵심 통찰. Logit 모형에서 \(\Psi\) 가 비선형이라 적분이 평균과 교환 불가:

\[ \Psi(E[z]) \neq E[\Psi(z)] \]

(좌변이 subject-specific, 우변이 marginal.)

수치 예시 — 만약 모집단의 절반은 \(\theta = -2\), 절반은 \(\theta = +2\) 라면:

subject-specific at \(\theta = 0\): \(\Psi(0) = 0.5\).
marginal: \(0.5 \cdot \Psi(-2) + 0.5 \cdot \Psi(+2) \approx 0.5 \cdot 0.119 + 0.5 \cdot 0.881 = 0.5\).

이 경우 우연히 같지만, 비대칭 분포나 공변량이 추가되면 차이가 명확해진다. 실무에서는 목적에 맞춰 선택:

의사가 환자에게 설명: subject-specific (“당신 같은 환자라면”).
정책 결정·전체 효과: marginal (“전체 환자에서 평균”).

3.3 Zeger Marginalization — 단순한 근사

Logistic 모형의 marginalization 근사

랜덤 절편 모형에 대해, 추정된 회귀 계수를 다음으로 나누면 marginal 계수 근사:

\[ \widehat{\beta}_M \approx \widehat{\beta}_{SS} \, / \sqrt{k^2 \widehat{\sigma}_v^2 + 1} \]

여기서 \(k = \dfrac{16 \sqrt{3}}{15 \pi}\).

직관 — 식 (9.16) 의 역방향 + 보정

§ 9.5 에서 본 식 (9.16) — subject-specific 과 marginal 의 척도 차이:

\[ \beta_M = \frac{\beta_{SS}}{\sqrt{(\sigma_v^2 + \sigma_\epsilon^2) / \sigma_\epsilon^2}} \]

logit 의 \(\sigma_\epsilon^2 = \pi^2/3\) 그대로 쓰면 부정확. Zeger et al. (1988) 가 보정 인자 \(k = 16\sqrt{3}/(15\pi) \approx 0.588\) 을 도입:

\[ \sqrt{(\sigma_v^2 + \pi^2/3) / (\pi^2/3)} \approx \sqrt{k^2 \sigma_v^2 + 1} \]

근사의 의미: \(k^2 \sigma_v^2 = (16\sqrt{3} / 15\pi)^2 \sigma_v^2\) 가 \(\sigma_v^2 / (\pi^2/3)\) 와 거의 같음 (\(k^2 \approx 0.346\), \(1/(\pi^2/3) \approx 0.304\) — 약간 다르지만 가까움). 차이를 정확히 보정하는 인자가 \(k\).

왜 적분 없이 marginal 을 얻을 수 있는가: 일반적으로 식 (9.40) 의 적분으로 marginal 을 계산해야 하지만, 랜덤 절편 모형 + logistic link 의 특수 경우에서는 잠근 근사식이 존재. 실무에서 빠른 marginal 환산이 필요할 때 유용.

4 § 9.6.2 — Multiple Random Effects 의 추정

4.1 Cholesky 미분 — 식 (9.50-9.54)

다중 랜덤 효과 모형 (식 9.50)

§ 9.5.2 의 식 (9.24) 와 동일:

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = x_{ij}^\top \beta + z_{ij}^\top T \theta_i \tag{9.50} \]

추정 대상: 스칼라 \(\sigma_v\) 가 \(r \times r\) Cholesky 인자 \(T\) (하삼각, \(r^* = r(r+1)/2\) 자유 모수) 로 대체. 식 (9.46) 의 \(\partial z_{ij}/\partial \sigma_v = \theta_i\) 자리에 \(\partial z_{ij}/\partial T\) 가 들어간다.

왜 행렬 미분이 필요한가

랜덤 효과가 1 개일 때 (\(r = 1\)): \(T = \sigma_v\) 스칼라, 미분도 스칼라. 단순.

랜덤 효과가 \(r\) 개일 때: \(T\) 가 \(r \times r\) 하삼각 행렬, 자유 모수 \(r^* = r(r+1)/2\) 개. 각 모수에 대한 score 와 정보 행렬을 계산하려면 행렬 미분 도구가 필요.

핵심 도구 (Magnus 1988):

vec 연산자: 행렬을 벡터로 변환 (열을 차례로 쌓음).
Kronecker 곱 (\(\otimes\)): 행렬 미분에서 곱셈 규칙을 처리.
Elimination matrix \(J_r\): 하삼각 행렬의 자유 모수만 추출 (위쪽 0 부분 제거).

식 (9.51-9.54) — 도함수 도출

vectorize:

\[ \text{vec}(z_{ij}^\top T \theta_i) = (\theta_i^\top \otimes z_{ij}^\top) \cdot \text{vec}(T) \tag{9.51} \]

elimination matrix 로 자유 모수만 추출:

\[ \text{vec}(z_{ij}^\top T \theta_i) = (\theta_i^\top \otimes z_{ij}^\top) \, J_r^\top \, v(T) \tag{9.53} \]

도함수:

\[ \frac{\partial z_{ij}}{\partial v(T)} = (\theta_i^\top \otimes z_{ij}^\top) \, J_r^\top = [J_r (\theta_i \otimes z_{ij})]^\top \tag{9.54} \]

여기서 \(v(T)\) 는 \(T\) 의 하삼각 자유 모수 (\(r^*\) 개) 를 모은 벡터.

직관 — 행렬을 벡터로 펴고 자유 모수만 추출

복잡해 보이지만 핵심은 두 단계:

Vectorize: \(z_{ij}^\top T \theta_i\) 가 스칼라이지만 \(T\) 의 모수에 대한 미분을 계산하려면 \(T\) 를 벡터로 펴야 한다 (벡터 미적분 표준 표기). vec 연산자가 이 역할.
Elimination: \(T\) 가 하삼각이라 위쪽 (\(r(r-1)/2\) 개) 원소는 항상 0. 자유 모수는 하삼각의 \(r^*\) 개. Elimination matrix \(J_r\) 가 vec(\(T\)) 에서 자유 모수만 추출.

식 (9.54) 의 결과: 하나의 도함수 벡터가 \(r^* \times 1\) 차원 (자유 모수와 같은 수). 이 벡터를 식 (9.46) 의 \(\partial z / \partial \sigma_v = \theta_i\) 자리에 대입하면 multiple random effects 의 score 와 정보 행렬이 1 차원 경우와 같은 형태로 계산된다.

핵심 메시지 — 단일 랜덤 효과의 알고리즘이 행렬 미분 도구를 통해 다중 랜덤 효과로 곧바로 일반화된다. 새로운 알고리즘을 짜는 것이 아니라 기존 score 식 (9.44) 의 한 항만 교체.

5 § 9.6.3 — 랜덤 효과 분포에 대한 적분

5.1 적분 근사 방법의 풍경

식 (9.40) 의 적분을 어떻게 풀 것인가

\[ h(Y_i) = \int_\theta \ell(Y_i \mid \theta) \, g(\theta) \, d\theta \]

이 적분이 닫힌 형태로 안 풀리므로 다양한 근사 방법이 제안되었다.

4 가지 접근법 비교

1. Marginal Quasi-Likelihood (MQL) [Goldstein-Rasbash]

Taylor 전개를 모형의 fixed part 주변에서 1차 또는 2차로 펼침.
빠르지만 1차 전개는 분산 성분 추정에 큰 음의 편향 (Breslow-Lin 1995).
MLwiN 소프트웨어에서 사용.

2. Penalized/Predictive Quasi-Likelihood (PQL) [Goldstein-Rasbash]

Taylor 전개에 random part 도 포함.
MQL 보다 정확하지만 여전히 편향 가능.

3. Laplace Approximation [Raudenbush et al. 2000a]

적분을 정규 분포로 근사 (사후 분포의 mode 주변 Gaussian).
다변량 Taylor + Laplace 결합으로 정확도와 속도 균형.
HLM 소프트웨어에서 사용.
우도 비교 검정에 deviance 직접 사용 가능 (PQL 은 불가).

4. Gauss-Hermite Quadrature [Stroud-Sechrest 1966]

정규 분포 가정 하에 적분을 합으로 근사.
점 수 \(Q\) 를 늘리면 임의 정확도 가능.
차원 폭발 (\(Q^r\)) 이 단점 — adaptive quadrature 로 완화.
MIXOR, MIXNO, NLMIXED 소프트웨어에서 사용.

5. MCMC (Gibbs Sampling) [Geman-Geman 1984]

사후 분포에서 표본 추출.
가장 유연 (어떤 분포든 처리), 가장 느림.
BUGS, Stan, brms 등에서 사용.
다차원 랜덤 효과 (\(r > 5\)) 에 유리.

본 § 9.6.3 의 초점은 Gauss-Hermite quadrature.

5.2 Gauss-Hermite Quadrature — 식 (9.55-9.59)

핵심 아이디어

표준 정규 \(g(\theta) = \mathcal{N}(0, 1)\) 에 대한 적분을 \(Q\) 개 점·가중치의 합으로 정확히 근사:

\[ \int_{-\infty}^{\infty} f(\theta) \, e^{-\theta^2/2} \, d\theta \approx \sum_{q=1}^Q A_q \, f(B_q) \]

\(B_q\): 미리 계산된 노드 (Hermite 다항식의 근), \(A_q\): 미리 계산된 가중치. Stroud-Sechrest (1966) 표 참조.

이 방법은 다항식 적분에 대해 정확 (\(Q\) 점이면 차수 \(2Q-1\) 까지 정확). 일반 함수에 대해서도 매끄러운 함수면 빠르게 수렴.

다변량 확장 (식 9.55)

\(r\) 차원 랜덤 효과의 경우, 노드는 \(r\) 차원 벡터 \(B_q^\top = (B_{q1}, B_{q2}, \ldots, B_{qr})\), 가중치는 일변량 가중치의 곱:

\[ A(B_q) = \prod_{h=1}^r A(B_{qh}) \tag{9.55} \]

총 적분 점 수: \(Q^r\).

Quadrature 를 사용한 모형 (식 9.56-9.59)

선형 예측자에 quadrature 노드 대입:

\[ z_{ijq} = x_{ij}^\top \beta + z_{ij}^\top T B_q \tag{9.56} \]

각 노드에서 conditional likelihood:

\[ \ell(Y_i \mid B_q) = \prod_{j=1}^{n_i} \Psi(z_{ijq})^{Y_{ij}} [1 - \Psi(z_{ijq})]^{1 - Y_{ij}} \tag{9.57} \]

근사 marginal likelihood (식 9.40 의 합 근사):

\[ h(Y_i) \approx \sum_{q=1}^{Q^r} \ell(Y_i \mid B_q) \, A(B_q) \tag{9.58} \]

근사 score (식 9.44 의 합 근사):

\[ \frac{\partial \log L}{\partial \eta} \approx \sum_{i=1}^N h^{-1}(Y_i) \sum_{q=1}^{Q^r} \sum_{j=1}^{n_i} \frac{Y_{ij} - \Psi(z_{ijq})}{\Psi(z_{ijq})[1 - \Psi(z_{ijq})]} \, \partial \Psi(z_{ijq}) \, \frac{\partial z_{ijq}}{\partial \eta} \, \ell(Y_i \mid B_q) \, A(B_q) \tag{9.59} \]

직관 — 적분이 합으로 바뀌는 마법

식 (9.40) 의 적분 \(\int \cdots g(\theta) d\theta\) 가 식 (9.58) 의 합 \(\sum_q \cdots A(B_q)\) 으로 바뀌었다. 핵심은:

적분 변수 \(\theta\) → 미리 정해진 \(Q\) 개 노드 \(B_q\) 로 대체.
분포 \(g(\theta)\) → 미리 정해진 가중치 \(A(B_q)\) 로 대체.
모든 함수 호출은 노드 위치에서만 평가.

비유 — 곡선 아래 면적을 trapezoidal rule 로 근사할 때 \(x\) 축의 미세한 점들의 함수값을 더하는 것과 같다. 차이는 점들의 위치 (Hermite 다항식의 근, 정규 분포의 무게가 큰 곳에 집중) 와 가중치가 최적화되어 있다는 것.

왜 정규 가정이 필수인가: Gauss-Hermite 의 노드·가중치는 \(e^{-\theta^2/2}\) 형태의 분포를 가정해 미리 계산. 다른 분포 (uniform, Cauchy 등) 도 가능하지만 (Bock-Aitkin 1981 의 비모수 분포 등), 정규가 가장 효율.

이 합이 정확해지려면 \(Q\) 가 충분히 커야 한다. Longford (1993) 와 Jansen (1990): 일변량 모형은 \(Q = 5\) 이상이면 거의 정확. 다차원은 차원당 점 수를 줄여도 됨 — Bock et al. (1988) 의 5 차원 인수분석에서 차원당 3 점만으로 충분.

LIMDEP·EGRET 의 기본값은 일변량 20 점 — 안전 마진. Lesaffre-Spiessens (2001) 는 매우 높은 점 수가 필요한 사례를 보여 — \(Q\) 에 대한 결과의 의존성을 항상 점검 권고.

5.3 차원 폭발 (\(Q^r\) 문제)

차원의 저주 — Curse of Dimensionality

총 적분 점 수 \(Q^r\):

\(r\) (랜덤 효과 수)	\(Q = 5\)	\(Q = 10\)	\(Q = 20\)
1	5	10	20
2	25	100	400
3	125	1,000	8,000
4	625	10,000	160,000
5	3,125	100,000	3,200,000
6	15,625	1,000,000	64,000,000

각 점에서 score·정보 행렬을 계산하므로 \(r > 5\) 면 실용적 한계.

직관 — 차원이 늘면 적분 영역이 폭발적으로 커진다

1 차원에서는 \(\theta\) 가 실수 한 개라 \(Q\) 점이면 충분. 2 차원에서는 \(\theta = (\theta_1, \theta_2)\) 평면이라 격자 점 수가 \(Q^2\). \(r\) 차원에서는 \(r\) 차원 hypercube 의 격자 점이 \(Q^r\).

비유 — 1 차원 도로는 100 m 마다 점 찍으면 1 km 가 10 점. 2 차원 평면은 같은 밀도면 100 점. 3 차원 공간은 1000 점. 차원이 늘수록 같은 정확도를 위해 폭발적으로 더 많은 점 필요.

해결책 두 가지:

차원당 점 수 줄이기: Bock et al. (1988) 처럼 \(r = 5\) 이면 \(Q = 3\) → \(3^5 = 243\) 점 (적당). 정확도 약간 손실.
Adaptive quadrature: 사후 분포의 위치·분산에 맞춰 점을 재배치 → 적은 점으로 같은 정확도. 다음에 설명.

5.4 Adaptive Gauss-Hermite Quadrature — 식 (9.60-9.61)

식 (9.60-9.61) — Adaptive 노드와 가중치

각 환자별로 사후 평균 \(\bar{\theta}_i\) (즉 식 9.48 의 \(\widehat{\theta}_i\)) 와 사후 표준편차 \(s_i\) 를 추정 (각 iteration 마다 갱신).

Adaptive 노드:

\[ B_{iq} = \bar{\theta}_i + s_i \cdot B_q \tag{9.60} \]

Adaptive 가중치:

\[ A_{iq} = \sqrt{2\pi} \cdot s_i \cdot \exp(B_q^2 / 2) \cdot \phi(B_{iq}) \cdot A_q \tag{9.61} \]

여기서 \(\phi\) 는 표준 정규 pdf, \(B_q, A_q\) 는 표준 Gauss-Hermite 노드·가중치.

직관 — “각 환자의 사후 분포에 맞춰 점 재배치”

표준 Gauss-Hermite 의 노드는 표준 정규 \(\mathcal{N}(0, 1)\) 에 최적화. 그러나 실제로 적분해야 하는 분포는 \(f(\theta \mid Y_i) \propto \ell(Y_i \mid \theta) g(\theta)\) — 환자별 사후 분포. 이 분포가 \(\mathcal{N}(\bar{\theta}_i, s_i^2)\) 에 가까운 경우가 많다.

발상: 각 환자의 사후 분포 위치 (\(\bar{\theta}_i\)) 와 폭 (\(s_i\)) 에 맞춰 노드를 재배치 + 가중치를 보정. 식 (9.60) 이 노드 재배치 (\(\bar{\theta}_i\) 만큼 이동, \(s_i\) 만큼 폭 조정), 식 (9.61) 이 가중치 보정 (Jacobian + 분포 변경에 따른 재가중).

결과: 차원당 3-5 점만으로도 정확한 적분 (Pinheiro-Bates 1995, Bock-Shilling 1997). \(Q^r\) 폭발이 크게 완화.

비유 — 표준 quadrature 가 “지도 상의 격자점에서 측정” 이라면, adaptive quadrature 는 “관심 지역 (사후 분포 봉우리) 주변에 점을 집중 배치”. 같은 점 수로 더 정확한 결과.

매 iteration 마다 \(\bar{\theta}_i, s_i\) 를 갱신해야 하므로 한 iteration 의 비용은 더 크지만, 점 수가 훨씬 적어 전체적으로 빠름. 다차원 랜덤 효과 (\(r \geq 3\)) 에서 표준보다 압도적으로 유리.

현재 표준 도구:

R lme4::glmer(..., nAGQ = N): \(N \geq 1\) 이면 adaptive Gauss-Hermite quadrature (\(N\) 점). 단 \(r = 1\) 일 때만 \(N > 1\) 가능.
Stata meqrlogit, meglm: 다차원 adaptive quadrature.
SAS PROC NLMIXED, PROC GLIMMIX (METHOD=QUAD): 비슷.

5.5 소프트웨어 비교

주요 GLMM 패키지

소프트웨어	적분 방법	다중 랜덤 효과	분포 옵션	비고
EGRET	Gauss-Hermite	단일만	정규	이항만 지원
LIMDEP	Gauss-Hermite	단일만	정규	기본 20 점
MIXOR	Gauss-Hermite	다중	정규 또는 균일	Hedeker-Gibbons 1996a
MIXNO	Gauss-Hermite	다중	정규 또는 균일	Hedeker 1999, 명목형
SAS PROC NLMIXED	(Adaptive) Gauss-Hermite	다중	정규	일반 비선형 모형
SAS PROC GLIMMIX	PQL 기본, QUAD 옵션	다중	정규	Adaptive quadrature 옵션
HLM	Laplace	다중 (2-3 수준)	정규	LR test 가능
MLwiN	MQL/PQL	다중 (3 수준 이상)	정규	빠르나 편향 가능
BUGS / Stan	Gibbs / HMC	다중 (제한 없음)	자유	가장 유연, 가장 느림
R lme4 (glmer)	Laplace 또는 AGQ	다중 (단 AGQ 는 \(r = 1\) 만)	정규	가장 보편
R glmmTMB	Laplace	다중	정규	빠른 대안
Python statsmodels	VB / Laplace	제한적	정규	발전 중

실무 권고 — 어느 도구를 언제

랜덤 효과 1-2 개: R lme4::glmer(..., nAGQ = 10 ~ 25) 또는 Stata meqrlogit — adaptive quadrature 로 정확.

랜덤 효과 3-5 개: SAS PROC NLMIXED (METHOD=ISAMP) 또는 R glmmTMB — Laplace + Importance sampling.

랜덤 효과 6 개 이상: R brms (Stan) 또는 BUGS — MCMC 만 실용적.

빠른 prototype: R lme4::glmer(..., nAGQ = 1) (Laplace) — 정확도는 낮지만 빠른 탐색용.

\(Q\) 에 대한 민감도 점검 (Lesaffre-Spiessens 권고): 같은 모형을 \(Q = 5, 10, 20, 30\) 으로 적합하고 추정값이 안정될 때까지 늘림. 안정되지 않으면 모형 식별성·데이터 한계 의심.

6 응용 분야

분야	활용	추정 도구 권고
임상시험 종단 (이항 endpoint)	약물 효과 + 환자 이질성 추정	R `lme4::glmer` (단순 모형) 또는 SAS `NLMIXED` (복잡)
다기관 임상시험 (3 수준)	환자·기관·시점의 중첩 구조	HLM (Laplace) 또는 SAS `GLIMMIX`
IRT (교육 측정)	응시자 능력 + 문항 난이도	R `mirt` 또는 SAS `IRT`
행동 유전학 (쌍둥이)	MZ/DZ tetrachoric correlation	OpenMx 또는 R `umx` (probit)
생태학 종단	동물/구역 반복 관측, 환경 효과	R `lme4` (논문 표준) 또는 `glmmTMB`
사회 조사 (다수준)	학교·학생 중첩	HLM, MLwiN

7 코드 예시

7.1 Step 1: Marginal Likelihood 의 적분 직접 계산 (numpy)

import numpy as np
from scipy.special import expit  # logistic cdf
from scipy.integrate import quad


def conditional_likelihood(theta: float, y: np.ndarray, x: np.ndarray,
                            beta: np.ndarray, sigma_v: float) -> float:
    """식 (9.39): 한 환자의 conditional likelihood (theta 주어진 상태).

    y: 한 환자의 시점별 응답 (n_i 차원).
    x: 한 환자의 시점별 공변량 (n_i × p).
    """
    z = x @ beta + sigma_v * theta
    p = expit(z)
    # 곱 형태 (수치 안정 위해 log 후 exp)
    log_lik = np.sum(y * np.log(p + 1e-300) + (1 - y) * np.log(1 - p + 1e-300))
    return np.exp(log_lik)


def marginal_likelihood_quad(y: np.ndarray, x: np.ndarray,
                              beta: np.ndarray, sigma_v: float) -> float:
    """식 (9.40): scipy.integrate.quad 로 적분."""
    integrand = lambda theta: conditional_likelihood(theta, y, x, beta, sigma_v) \
                              * np.exp(-theta**2 / 2) / np.sqrt(2 * np.pi)
    result, _ = quad(integrand, -8, 8)  # ±8σ 면 충분
    return result


# 더미 데이터 (한 환자, 5 시점)
y_demo = np.array([1, 1, 0, 1, 0])
x_demo = np.column_stack([np.ones(5), np.arange(5)])  # 절편 + 시간
beta_demo = np.array([-0.5, 0.3])
sigma_v_demo = 1.2

h_y = marginal_likelihood_quad(y_demo, x_demo, beta_demo, sigma_v_demo)
print(f"Marginal likelihood h(Y_i) = {h_y:.6e}")
print(f"Log marginal = {np.log(h_y):.4f}")

검증 포인트

scipy.integrate.quad 는 적응적 적분 — Gauss-Kronrod 알고리즘. 이 결과가 Gauss-Hermite quadrature 의 정답.

실제 추정에서는 환자 수 \(N \times\) iteration 수 만큼 이 적분을 반복해야 하므로 quad 는 너무 느리다. Gauss-Hermite quadrature (다음 단계) 가 같은 결과를 훨씬 빠르게 산출.

7.2 Step 2: Gauss-Hermite Quadrature 직접 구현 (numpy)

import numpy as np
from numpy.polynomial.hermite_e import hermegauss
from scipy.special import expit


def marginal_likelihood_ghq(y: np.ndarray, x: np.ndarray,
                             beta: np.ndarray, sigma_v: float,
                             Q: int = 20) -> float:
    """식 (9.58): Gauss-Hermite quadrature 로 식 (9.40) 근사.

    표준 정규 N(0,1) 에 대해 hermegauss(Q) 로 노드·가중치 산출.
    """
    # 표준 정규에 대한 Hermite_e 노드·가중치
    nodes, weights = hermegauss(Q)
    # 가중치는 e^(-x²/2)/√(2π) 와 일치하도록 정규화
    weights = weights / np.sqrt(2 * np.pi)

    h = 0.0
    for q in range(Q):
        z = x @ beta + sigma_v * nodes[q]  # 식 (9.56)
        p = expit(z)
        log_lik = np.sum(y * np.log(p + 1e-300) + (1 - y) * np.log(1 - p + 1e-300))
        h += weights[q] * np.exp(log_lik)
    return h


# 같은 데이터에 적용
y_demo = np.array([1, 1, 0, 1, 0])
x_demo = np.column_stack([np.ones(5), np.arange(5)])
beta_demo = np.array([-0.5, 0.3])
sigma_v_demo = 1.2

# Q 별 결과 — 수렴 확인
for Q in [3, 5, 10, 20, 40]:
    h = marginal_likelihood_ghq(y_demo, x_demo, beta_demo, sigma_v_demo, Q)
    print(f"Q = {Q:3d}: h(Y_i) = {h:.6e}")

검증 포인트 — Q 수렴

\(Q\) 가 5 이상이면 결과가 안정되어야 한다 (Longford 1993, Jansen 1990 의 권고와 일치). \(Q = 20\) 결과가 quad 결과와 거의 동일.

Step 1 의 quad 와 Step 2 의 GHQ (\(Q = 20\)) 가 일치하면 quadrature 가 정확. 이 코드를 모형 적합 코드의 inner loop 에 두면 marginal MLE 가 만들어진다.

7.3 Step 3: Empirical Bayes 추정 (numpy)

import numpy as np
from numpy.polynomial.hermite_e import hermegauss
from scipy.special import expit


def empirical_bayes_theta(y: np.ndarray, x: np.ndarray,
                           beta: np.ndarray, sigma_v: float,
                           Q: int = 30) -> tuple:
    """식 (9.48-9.49): EAP 추정과 사후 분산.

    Returns: (theta_hat, var_theta_hat).
    """
    nodes, weights = hermegauss(Q)
    weights = weights / np.sqrt(2 * np.pi)

    # 각 노드에서 conditional likelihood
    h = 0.0  # marginal likelihood (식 9.40)
    numerator = 0.0  # ∫ θ · ℓ(Y|θ) g(θ) dθ
    for q in range(Q):
        z = x @ beta + sigma_v * nodes[q]
        p = expit(z)
        log_lik = np.sum(y * np.log(p + 1e-300) + (1 - y) * np.log(1 - p + 1e-300))
        ell = np.exp(log_lik)
        h += weights[q] * ell
        numerator += weights[q] * nodes[q] * ell

    theta_hat = numerator / h  # 식 (9.48)

    # 사후 분산 (식 9.49)
    var_numerator = 0.0
    for q in range(Q):
        z = x @ beta + sigma_v * nodes[q]
        p = expit(z)
        log_lik = np.sum(y * np.log(p + 1e-300) + (1 - y) * np.log(1 - p + 1e-300))
        ell = np.exp(log_lik)
        var_numerator += weights[q] * (nodes[q] - theta_hat) ** 2 * ell

    var_theta_hat = var_numerator / h
    return theta_hat, var_theta_hat


# 예: 응답이 모두 1 인 환자 → 양수 theta 예상
y_high = np.array([1, 1, 1, 1, 1])
y_low = np.array([0, 0, 0, 0, 0])
y_mixed = np.array([1, 0, 1, 0, 1])

x_demo = np.column_stack([np.ones(5), np.arange(5)])
beta_demo = np.array([0.0, 0.0])  # 절편 0, 시간 효과 0 (단순화)
sigma_v_demo = 1.5

for label, y in [("high", y_high), ("low", y_low), ("mixed", y_mixed)]:
    theta_hat, var_theta = empirical_bayes_theta(y, x_demo, beta_demo, sigma_v_demo)
    print(f"y = {label:5s}: theta_hat = {theta_hat:+.3f}, "
          f"SD(theta_hat) = {np.sqrt(var_theta):.3f}")

결과 해석 — Shrinkage 의 작동

y_high (모두 1): \(\widehat\theta > 0\), 양수 효과 추정.
y_low (모두 0): \(\widehat\theta < 0\), 음수 효과 추정.
y_mixed (반반): \(\widehat\theta \approx 0\), 사전 분포로 수축.

같은 환자에 시점이 더 많으면 (예: \(n_i = 20\)) → 사후 분산 감소, 추정이 데이터에 가까워짐. 같은 환자에 시점이 적으면 (예: \(n_i = 1\)) → 사후 분산 큼, 추정이 사전 (0) 으로 수축.

이 shrinkage 가 mixed-effects 의 핵심 — fixed-effects 라면 환자별 자유 추정으로 변동이 폭발하지만, mixed-effects 는 사전 분포가 자동으로 안정화.

7.4 Step 4: 실무 적합 — R `lme4::glmer` 비교 (Python에서 호출)

import numpy as np
import pandas as pd
import statsmodels.api as sm


# 종단 이항 데이터 시뮬레이션 (앞 09-3 의 함수 재사용)
def simulate_longitudinal_binary(n_subjects: int, n_times: int,
                                  beta0: float, beta_treat: float,
                                  beta_time: float, sigma_v: float,
                                  seed: int = 2026) -> pd.DataFrame:
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n_subjects):
        treat = rng.binomial(1, 0.5)
        upsilon = rng.normal(0, sigma_v)
        for j in range(n_times):
            t = j
            eta = beta0 + beta_treat * treat + beta_time * t + upsilon
            p = 1.0 / (1.0 + np.exp(-eta))
            y = rng.binomial(1, p)
            rows.append({"subject": i, "time": t, "treat": treat, "y": y})
    return pd.DataFrame(rows)


df = simulate_longitudinal_binary(n_subjects=300, n_times=5,
                                   beta0=-1.0, beta_treat=0.8,
                                   beta_time=0.3, sigma_v=1.2)

# Python statsmodels (variational Bayes)
model = sm.BinomialBayesMixedGLM.from_formula(
    "y ~ treat + time", random={"a": "0 + C(subject)"}, data=df)
result = model.fit_vb()

print("True parameters: beta0=-1.0, beta_treat=0.8, beta_time=0.3, sigma_v=1.2")
print("\nstatsmodels VB estimates:")
print(f"  beta:    {result.fe_mean.round(3).tolist()}")
print(f"  sigma_v: {np.sqrt(result.vc_mean[0]):.3f}  (estimated)")

실무 비교 — Python vs R

R 등가 코드:

library(lme4)
fit <- glmer(y ~ treat + time + (1 | subject),
             data = df, family = binomial,
             nAGQ = 20)  # adaptive Gauss-Hermite quadrature, 20 점
summary(fit)

# Empirical Bayes 랜덤 효과
ranef_estimates <- ranef(fit)$subject

Python statsmodels.BinomialBayesMixedGLM 은 변분 근사 (VB) 라 정확도가 quadrature 보다 낮다. R lme4::glmer(nAGQ = 20) 또는 Stata meqrlogit 가 정확도 측면에서 표준.

빠른 결과가 필요하면 R lme4 Laplace (nAGQ = 1) 도 가능 — 작은 \(\sigma_v\) 에서는 정확하나 클 때 편향.

8 관련 주제

선행 지식

§ 9.5.1-9.5.3 — Mixed-effects logistic 모형 정의와 ICC
§ 9.5.4-9.5.5 — Multilevel form 과 response functions
§ 4.5 MRM 추정 — 정규 MRM 의 ML/REML/EM/Fisher Scoring (구조적 비교)
§ 9.2-9.3 — Single-level Bernoulli MLE 의 score (식 9.7-9.9), Newton-Raphson

후속 주제

§ 9.7 — Illustration: 정신과 임상시험 데이터에 mixed-effects logistic 적용
§ 9.8-9.10 — 추정 세부 (Hessian, optimizer, 표준오차)
§ 9.11 — Subject-specific vs population-averaged 효과의 정량 비교 (Zeger marginalization 의 활용)
Ch.10 GLMM 순서형 — 누적 logit 의 추정 (적분 차원 동일, link 만 다름)
Ch.11 GLMM 명목 — 다항 logit 의 추정 (적분 차원 = 범주 - 1)
Ch.12 GLMM 카운트 — Poisson/NB 의 추정

관련 개념

Bock & Aitkin (1981) — Marginal MLE + Gauss-Hermite quadrature 의 원전
Pinheiro & Bates (1995) — Adaptive quadrature 의 발전
Raudenbush et al. (2000a, 2000b) — Laplace approximation 의 정확도
Breslow & Lin (1995) — PQL 의 편향 (Bock-Aitkin quadrature 정당화)
Liu & Pierce (1994), Rabe-Hesketh et al. (2002) — Adaptive quadrature 이론
GEE (Ch.8) — Marginal 우도를 회피하는 대안 (작동 상관 + sandwich)