1 서론 — LRT 의 \(\chi^2\) 근사에는 오차가 있다
Ch.1-14 전체에서 반복적으로 등장한 도구 — 우도비 검정 (LRT):
\[ \Lambda = 2\{l(\widehat\theta) - l(\theta_0)\}. \]
귀무 하 점근 분포: \(\Lambda \overset{d}{\to} \chi_p^2\). 이 \(\chi^2\) 근사 덕분에 표준 가설 검정이 작동한다.
그러나 유한 표본에서 오차:
\[ E(\Lambda) = p + \epsilon_p + O(n^{-2}), \qquad \epsilon_p = O(n^{-1}). \]
즉 \(\chi^2\) 근사는 평균부터 틀리기 시작. \(p\) 대신 \(p(1 + b_p)\).
Bartlett 의 아이디어 (1937): 보정 통계량
\[ \Lambda' = \frac{\Lambda}{1 + b_p}. \]
\(\Lambda'\) 의 모든 cumulant 가 \(\chi_p^2\) 의 cumulant 와 \(O(n^{-2})\) 까지 일치. 즉 한 차수 더 정확한 \(\chi^2\) 근사.
이 포스트는 overview (14-1) 에서 개요만 다룬 §15.3 을 완전 유도 + 컴퓨터 구현 + Feigl-Zelen 수치 재현 으로 심화한다. 핵심은 (a) 단일 스칼라 \(b_p\) 가 왜 모든 cumulant 를 동시 보정하는지, (b) GLM 의 6 개 불변 스칼라가 어떤 수학적 구조를 공유하는지, (c) 실무에서 Bartlett 가 언제 큰 차이를 만드는지.
2 모든 cumulant 동시 보정 — 놀라운 수학
2.1 기본 결과 (15.6)
McCullagh-Nelder 가 기록하는 식 (15.6):
\[ \kappa_r(\Lambda) = (r-1)! \cdot 2^{r-1} \cdot p \cdot (1 + b_p)^r + O(n^{-2}), \quad r \geq 1. \]
\(r\) = cumulant 차수. 우변의 \((r-1)! \cdot 2^{r-1} \cdot p\) 는 \(\chi_p^2\) 의 \(r\)-차 cumulant.
2.2 왜 이것이 놀라운가
일반적으로 분포의 모양을 고칠 때 여러 모수 조정이 필요하다: - 평균 맞추기 → 1 모수 - 분산 맞추기 → 1 모수 더 - 왜도 맞추기 → 1 모수 더 - 첨도 맞추기 → 1 모수 더 …
그런데 (15.6) 은 단 하나 의 스칼라 \(b_p\) 만 바꾸면 모든 cumulant 가 맞춰진다고 말한다.
2.3 비결 — 곱셈적 구조
(15.6) 의 핵심 구조는 우변의 \((1 + b_p)^r\). \(\Lambda' = \Lambda/(1 + b_p)\) 의 \(r\)-차 cumulant:
\[ \kappa_r(\Lambda') = \frac{\kappa_r(\Lambda)}{(1 + b_p)^r} = (r-1)! \cdot 2^{r-1} \cdot p + O(n^{-2}). \]
분모의 \((1+b_p)^r\) 이 분자의 \((1+b_p)^r\) 과 정확히 상쇄. 모든 차수에서 동시 일치.
2.4 확률 변수의 관점
\(\Lambda \approx (1 + b_p) \cdot \chi_p^2\) 가 스케일 변환 관계라는 뜻. LRT 통계량이 점근적으로 “\(\chi^2\) 의 배수 \(1+b_p\)” 이므로, 그 배수로 나누면 순수 \(\chi^2\).
점근 이론은 \(\Lambda \to \chi^2_p\) 을 보장. 유한 표본에서 \(\Lambda\) 는 \((1 + b_p)\) 배 부풀려져 있다. 오직 scale 에만 미세한 편향 — shape 은 이미 \(\chi^2\) 이다.
이 scale bias 를 고치는 것이 Bartlett. 다른 분포 교정 방법 (예: Edgeworth 전개) 은 shape 도 수정하지만 Bartlett 은 scale 만. 단순성 이 매력.
2.5 수학적 증명 스케치
Appendix C (McCullagh 1987, Ch.7) 의 텐서 미적분으로 증명. 핵심 도구:
- 정보 부등식 — MLE 의 공분산 = \(I^{-1}\).
- 3차·4차 cumulant — \(\kappa_{rst}, \kappa_{rs,tu}, \kappa_{rstu}\) 등 Bartlett 식별자.
- Edgeworth 전개 — \(\Lambda\) 의 밀도를 고차 항으로 확장.
- 미적분 상쇄 — 특정 조합이 \((1 + b_p)^r\) 구조로 정확히 인수분해.
이 증명은 텐서 수학의 걸작. 직관적으로는 “단순” 해 보이는 보정이 뒷단에서 거대한 수학적 구조 위에 서 있다.
3 복합 귀무가설 (15.7-15.8)
3.1 일반적 상황
실무 대부분의 가설 검정은 중첩 모형 비교: - \(H_0\): \(\theta\) 가 \(q\)-차원 부분공간. - \(H_1\): \(\theta\) 가 \(p\)-차원 전체공간 (\(p > q\)).
우도비 (15.7):
\[ \Lambda(\widehat\theta, \widehat\theta_0) = 2\{l(\widehat\theta) - l(\widehat\theta_0)\} = D(Y; \widehat\theta_0) - D(Y; \widehat\theta). \]
여기서 \(\widehat\theta_0\) 는 \(H_0\) 하의 MLE. 이탈도 차이로 표현.
3.2 기대값
두 개별 보정을 조합:
\[ E(\Lambda_{H_0,H_1}) = \{p + \epsilon_p\} - \{q + \epsilon_q\} + O(n^{-2}) = (p-q)\{1 + b_{pq}\} + O(n^{-2}). \]
복합 Bartlett 인수 (15.8):
\[ \boxed{\; b_{pq} = \frac{p b_p - q b_q}{p - q} = \frac{\epsilon_p - \epsilon_q}{p - q} \; } \]
두 전체 보정의 가중 차이. 가중치는 각 공간의 차원.
3.3 조정된 통계량
\[ \Lambda' = \frac{\Lambda(\widehat\theta, \widehat\theta_0)}{1 + b_{pq}} \sim \chi^2_{p-q} + O(n^{-2}). \]
자유도 = \(p - q\) (두 모형의 차원 차이).
3.4 계산 흐름
실무에서 \(b_{pq}\) 계산:
Step 1: H_1 (p-모형) 에서 ε_p 계산 — 본문 (15.9)
Step 2: H_0 (q-모형) 에서 ε_q 계산 — 동일 공식 적용
Step 3: b_{pq} = (ε_p - ε_q) / (p - q)
Step 4: Λ' = Λ / (1 + b_{pq})
Step 5: Λ' 을 χ²_{p-q} 임계값과 비교각 모형에서 \(\epsilon\) 계산이 개별 로 이뤄진다. 두 결과의 차이로 \(b_{pq}\).
4 §15.3.2 — \(\epsilon_p\) 의 6 개 불변 스칼라 구조
4.1 도입 행렬
GLM 의 \(\epsilon_p\) 계산에 사용되는 대각 행렬:
\[ D^{(1)} = \text{diag}\{\mu'_i\}, \quad \mu'_i = d\mu_i/d\eta_i. \]
\[ D^{(2)} = \text{diag}\{\mu''_i - \mu'^2_i \, d\log V_i/d\mu_i\}. \]
\[ W = \text{diag}\{\mu'^2_i / V_i\} = D^{(1)} V^{-1} D^{(1)}. \]
\(W\) = 표준 IRLS 가중치. \(D^{(1)}\) = 역링크의 1차 미분. \(D^{(2)}\) = 비대칭 / 곡률 보정.
4.2 두 사영 행렬
\[ Q = X(X^T W X)^{-1} X^T, \quad P = D^{(1)} Q D^{(1)} V^{-1}. \]
\(Q\): 대칭 · \(\widehat\eta\) 의 점근 공분산 (스케일 \(\sigma^2\)) 까지 결정.
\(P\): 비대칭 사영 · \(\widehat\mu\) 의 점근 공분산을 원 척도에서 표현. \(D^{(1)} V^{-1}\) 가중이 비대칭 유발.
4.3 Mixed Cumulants 의 특별한 대칭성
GLM 의 “놀라운 성질” (McCullagh-Nelder):
\[ \kappa_{rs,t} = \kappa_{rt,s} = \kappa_{st,r}, \quad \ldots \]
모든 지수 순열에 대해 대칭. 일반 모형에서는 성립하지 않는 GLM 특유의 구조.
일반 모형의 mixed cumulants 는 다른 종류: - \(\kappa_{rs,t}\): 2차 미분과 1차 미분의 공분산. - \(\kappa_{rt,s}\): 2차 미분과 1차 미분의 공분산 (지수 순서 다름). - 이 둘은 일반적으로 다를 수 있다.
GLM 의 경우 지수족 + 선형 예측자 의 특수 구조로 인해 두 값이 같아진다. \(x_r^i x_s^i x_t^i\) 의 곱 은 지수 순서에 무관하기 때문.
이 성질이 (15.9) 의 6 스칼라 식을 닫힌 매트릭스 형태 로 쓸 수 있게 해 준다. GLM 밖에서는 훨씬 복잡한 텐서 계산 필요.
4.4 Residual Second Derivative Matrix
\(V_{rs}\) = “\(U_{rs}\) 를 \(U_r\) 로 회귀한 잔차”:
\[ V_{rs} = \sum_i x_r^i x_s^i D_i^{(2)} V_i^{-1} \{\delta_{ij} - P_{ij}\} (Y_j - \mu_j). \]
실무 계산에서 근사:
\[ V_{rs} \simeq \sum_i x_r^i x_s^i D_i^{(2)} V_i^{-1} (Y_i - \widehat\mu_i). \]
이 양이 \(\epsilon_p\) 의 “잔차 Hessian 성분” 을 담는다.
4.5 6 개 불변 스칼라
McCullagh-Cox (1986) 의 공식:
(a) \(\rho_4\) — 4차 cumulant 기여:
\[ \rho_4 = \sum_i P_{ii}^2 \rho_{4i}, \qquad \rho_{4i} = \kappa_{4i} / \kappa_{2i}^2. \]
\(\rho_{4i}\) = 관측치 \(i\) 의 표준화 첨도 (정규에서 0).
(b) \(\rho_{13}^2\) — 대각 집중 왜도:
\[ \rho_{13}^2 = \sum_{ij} (P_{ii} \kappa_{3i}/\kappa_{2i}) V_i^{-1} P_{ij} (P_{jj} \kappa_{3j}/\kappa_{2j}). \]
대각 원소 \(P_{ii}\) 의 3차 왜도 상관 구조.
(c) \(\rho_{23}^2\) — 쌍별 왜도:
\[ \rho_{23}^2 = \sum_{ij} (V_i^{-1} P_{ij})^3 \kappa_{3i} \kappa_{3j}. \]
\(P_{ij}\) 의 3 제곱에 비례.
(d) 잔차 Hessian 분산 I:
\[ \nu_{rs,tu} \kappa^{r,s} \kappa^{t,u} = q^T D^{(2)} V^{-1} (I - P) D^{(2)} q, \]
\(q\) = \(Q_{ii}\) 벡터.
(e) 잔차 Hessian 분산 II:
\[ \nu_{rs,tu} \kappa^{r,t} \kappa^{s,u} = \sum_{ij} Q_{ij}^2 [D^{(2)} V^{-1} (I - P) D^{(2)}]_{ij}. \]
(f) 혼합 잔차:
\[ \nu_{r,s,tu} \kappa^{r,s} \kappa^{t,u} = q^T D^{(2)} V^{-1} (I - P) q^*, \]
\(q^*_j = q_j W_j \kappa_{3j}/\kappa_{2j}\).
4.6 최종 공식
Bartlett 식 (15.9):
\[ \boxed{\; \epsilon_p = -\frac{1}{4}(a) + \frac{1}{4}(b) + \frac{1}{6}(c) - \frac{1}{4}(d) + \frac{1}{2}(e) - \frac{1}{2}(f). \;} \]
6 스칼라의 선형 결합. 각 계수는 Cox-McCullagh 유도의 구체적 상수.
4.7 Canonical Link 의 단순화
McCullagh-Nelder: canonical link 에서는 \(D^{(2)} = 0\).
유도: canonical link 에서 \(g'(\mu) V = 1\), 즉 \(V_i = 1/g'(\mu_i) = \mu'_i\). 따라서 \(d \log V_i/d\mu_i = \mu''_i/\mu'^2_i\). \(D^{(2)}\) 정의 대입:
\[ D_i^{(2)} = \mu''_i - \mu'^2_i \cdot \frac{\mu''_i}{\mu'^2_i} = 0. \]
결과: (d), (e), (f) 항 모두 0. \(\epsilon_p\) 가 \(-\frac{1}{4}(a) + \frac{1}{4}(b) + \frac{1}{6}(c)\) 의 3 항 선형 결합.
Canonical link 의 정의적 특성: 로그우도가 정준 모수 \(\theta\) 에 대해 2차에서 정확. 3차 이상 미분 항이 “순수 데이터” 에 대한 선형 함수.
이 성질이 Bartlett 계산에서 \(D^{(2)}\) 항의 사라짐을 유발. 계산 간단화 + 근사 정확도 상승.
실무 권고: 가능하면 canonical link 사용. Logit (이항), log (포아송), inverse (gamma) 등.
5 §15.3.3 — Exponential Regression Example
5.1 모형 (15.11)
\(Y_1, \ldots, Y_n\) 독립 지수분포, \(\eta_i = \log \mu_i = \beta_0 + \beta_1 x_i\).
지수분포: \(f(y; \mu) = \mu^{-1} e^{-y/\mu}\). \(V(\mu) = \mu^2\) (constant CV family, Ch.8 Gamma 의 \(\nu = 1\) 특수 사례).
5.2 매트릭스 계산
지수 회귀에서:
\[D^{(1)} = \text{diag}\{\mu_i\}, \quad D^{(2)} = \text{diag}\{-\mu_i\}, \quad V_i = \mu_i^2.\]
\(W = D^{(1)} V^{-1} D^{(1)} = \text{diag}\{\mu_i \cdot \mu_i^{-2} \cdot \mu_i\} = I\) — 단위 가중치.
\(Q = X (X^T X)^{-1} X^T\) — 표준 선형 회귀 햇 행렬.
\(P = D^{(1)} Q D^{(1)} V^{-1} = \text{diag}\{\mu_i\} Q \text{diag}\{\mu_i^{-1}\}\). 비대각 · 비대칭 이지만 \(P_{ii} = Q_{ii}\) (대각 원소는 scaling 으로 변하지 않음).
5.3 특수 구조의 활용
\(\kappa_{2i} = \mu_i^2, \kappa_{3i} = 2\mu_i^3, \kappa_{4i} = 6\mu_i^4\) (지수 분포 cumulants).
\(\rho_{4i} = 6\mu_i^4 / \mu_i^4 = 6\), \(\kappa_{3i}/\kappa_{2i} = 2\mu_i\).
(a) \(\rho_4 = \sum P_{ii}^2 \cdot 6 = 6 \sum Q_{ii}^2\).
(b), (c) 도 \(P_{ij}\) 와 \(\mu\) 의 조합으로 단순화.
최종 결과 (McCullagh-Nelder 의 (15.10)):
\[ \boxed{\; \epsilon_p = \frac{1}{6} \sum_{ij} Q_{ij}^3 - \frac{1}{4} q^T (I - Q) q. \;} \]
5.4 Translation Invariance
지수 회귀는 translation type — 모수 변환 \((\beta_0, \beta_1) \to (\beta_0 + c, \beta_1)\) 에 대해 \(\mu_i \to e^c \mu_i\) 로 이동. \(Q\) 는 \(\mu\) 무관 이므로 \(\epsilon_p\) 도 \(\beta\) 무관. 데이터만의 함수.
5.5 Feigl-Zelen 백혈병 생존 데이터
17 명 AML (급성 골수성 백혈병) 환자. \(Y\) = 생존 시간 (주), \(x\) = log(초기 백혈구 수). 모형 (15.11):
\[\log \mu_i = \beta_0 + \beta_1 x_i.\]
5.6 단순 검정: \(H_0: \beta_1 = 0\)
원 LRT: \(\Lambda = 6.826\) on 1 df. \(p\)-value = \(0.89\%\) (매우 유의).
Bartlett 계산:
\[ Q_{ij} = \frac{1}{n} + \frac{(x_i - \bar x)(x_j - \bar x)}{\sum(x_i - \bar x)^2}. \]
\[ \sum Q_{ij}^3 = \frac{4}{n} + \frac{(\sum(x_i - \bar x)^3)^2}{(\sum(x_i - \bar x)^2)^3}. \]
\(\epsilon_2\) (2 모수 전체 모형):
\[ \epsilon_2 = (4/17 + 0.0000385)/6 - 0.0688/4 = 0.0220. \]
\(\epsilon_1\) (1 모수 귀무 모형): \(\epsilon_1 = 1/(6n) = 1/102 = 0.00980\).
\(b_{21}\) (식 15.8):
\[ b_{21} = \frac{2 \cdot (0.0220/2) - 1 \cdot (0.00980/1)}{2 - 1} = 0.0220 - 0.00980 = 0.0122. \]
조정된 통계량: \(\Lambda' = 6.826 / 1.0122 = 6.744\).
새 \(p\)-value: \(P(\chi_1^2 > 6.744) = 0.94\%\).
5.7 해석 — 단순 검정에서의 미미한 영향
0.89% → 0.94% — 사실상 같은 결론. Bartlett 보정이 필요 없음을 확증.
이 결과는 \(n = 17\) 에도 \(\chi^2\) 근사가 이미 정확함 을 보여 준다. 실무 함의: \(n\) 이 적당하면 Bartlett 생략 가능.
5.8 큰 Impact 의 예 — 모형 적합도 검정
관심을 바꿔 “모형 (15.11) 자체가 적합한가?” 질문. \(H_0\): 모든 \(\mu_i\) 가 자유 (\(p = 17\)). \(H_1\): 모형 (15.11) 로 제약 (\(q = 2\)).
LRT: \(\Lambda = 19.46\) on \(15 = 17 - 2\) df.
\(\epsilon_{17}\) = \(17/6 = 2.833\) (포화 모형).
\(b_{17,2}\):
\[ b_{17,2} = \frac{17 \cdot (2.833/17) - 2 \cdot (0.0220/2)}{17 - 2} = \frac{2.833 - 0.0220}{15} = 0.187. \]
조정된 통계량: \(\Lambda' = 19.46 / 1.187 = 16.39\).
\(p\)-value 변화: \(P(\chi_{15}^2 > 19.46) = 19.4\%\) → \(P(\chi_{15}^2 > 16.39) = 35.6\%\). 19% 감소.
5.9 해석 — 적합도 검정의 극적 변화
\(p\)-value 가 1.8 배 증가. 원 통계량이 한계 유의였다면 Bartlett 조정으로 무의 결론 가능.
영향 작음: - 가설 검정 자유도 \(p - q\) 작음 (1-2). - \(n\) 이 충분히 큼 (\(\geq 50\)). - Canonical link.
영향 큼: - \(p - q\) 큼 (적합도 검정, 여러 모수 동시 검정). - \(n\) 작음 또는 \(p \sim n\). - Non-canonical link. - 비정규 반응 (Gamma, 역 가우시안 등 긴 꼬리).
Feigl-Zelen 예제가 두 양상을 한 데이터에서 보여 준다. 실무자는 검정 종류에 따라 선별 적용.
6 Lattice Case — 이산 반응의 한계
6.1 이론적 한계
McCullagh-Nelder 경고: 이산 (lattice) 분포 (이항, 포아송 등) 의 Bartlett 보정은 연속 분포만큼 도움이 되지 않는다.
이유: lattice distribution 의 \(\chi^2\) 근사 오차는 \(O(n^{-1/2})\) 이 본질적 한계. Edgeworth 전개의 continuity correction 오차. Bartlett 의 \(O(n^{-2})\) 향상이 이 한계 아래로 내려갈 수 없다.
6.2 실무 함의
- 연속 반응 (정규, Gamma, 역 가우시안): Bartlett 효과 완전 실현. \(O(n^{-2})\) 정확도.
- 이항 · 포아송: Bartlett 해도 \(O(n^{-1/2})\) 한계. 개선 있지만 제한적.
6.3 이산 반응의 대안
- Edgeworth 전개: Bartlett 보다 구체적이지만 복잡.
- Saddlepoint 근사: 꼬리 정확도 뛰어남. Durbin (1980).
- Exact 검정: Fisher 정확 검정, 순열 검정.
- Bootstrap calibration: 시뮬레이션으로 실제 \(p\)-value 분포 추정.
7 현대적 대안 — Bootstrap-Bartlett
7.1 아이디어
이론적 \(\epsilon_p\) 계산이 복잡하면 bootstrap 으로 추정:
- \(\widehat\theta_0\) 에서 \(B\) 개 가상 데이터 생성.
- 각 데이터에서 \(\Lambda\) 계산.
- 평균 \(\bar\Lambda^*\).
- \(\widehat b_p = \bar\Lambda^* / p - 1\).
- \(\Lambda' = \Lambda / (1 + \widehat b_p)\).
7.2 장점
- 이론 유도 불필요 — 어떤 GLM 에도 적용.
- 모형 복잡도 무관 — 랜덤효과, nuisance, 혼합 모형에서도 가능.
- 직관적 — Monte Carlo 에 익숙한 분석가에게 접근성.
7.3 단점
- 계산 비용 — \(B = 1000\) 정도 필요.
- \(B\) 가 작으면 \(\widehat b_p\) 자체에 noise.
- 이론 공식이 존재하면 정확한 이론 \(b_p\) 가 더 정밀.
8 Python 실전 — Feigl-Zelen 데이터
8.1 데이터 입력
import numpy as np
import statsmodels.api as sm
import pandas as pd
from scipy import stats
# Feigl & Zelen (1965) AML 생존 데이터
# 각 튜플: (survival_weeks, log_WBC_count)
# 실제 데이터 (Cox-Snell 1981, Table 6.2)
data_fz = np.array([
(65, 3.36), (156, 2.88), (100, 3.63), (134, 3.41), (16, 3.78),
(108, 4.02), (121, 4.00), (4, 4.23), (39, 3.73), (143, 3.85),
(56, 3.97), (26, 4.51), (22, 4.54), (1, 5.00), (1, 5.00),
(5, 4.72), (65, 5.00),
])
df = pd.DataFrame(data_fz, columns=['weeks', 'log_wbc'])
print(df)8.2 지수 회귀 적합
# Gamma (shape=1 = exponential) + log link
X = sm.add_constant(df['log_wbc'])
m_full = sm.GLM(df['weeks'], X,
family=sm.families.Gamma(link=sm.families.links.log())).fit(scale=1)
m_null = sm.GLM(df['weeks'], np.ones(len(df)),
family=sm.families.Gamma(link=sm.families.links.log())).fit(scale=1)
LRT = 2 * (m_full.llf - m_null.llf)
print(f"LRT (β₁=0 vs full): {LRT:.3f}")
print(f"p-value: {1 - stats.chi2.cdf(LRT, 1):.4f}")기대: \(\Lambda \approx 6.83\), \(p \approx 0.009\).
8.3 Bartlett 인수 계산
def bartlett_epsilon_exp(X_design):
"""지수 회귀에서 (15.10) 으로 ε_p 계산."""
X = np.asarray(X_design)
n = X.shape[0]
# Q = X (X'X)^-1 X'
Q = X @ np.linalg.inv(X.T @ X) @ X.T
q = np.diag(Q) # Q_ii 벡터
# (15.10): ε_p = (1/6) Σ Q_ij^3 - (1/4) q^T (I - Q) q
term1 = np.sum(Q ** 3) / 6
term2 = q @ (np.eye(n) - Q) @ q / 4
return term1 - term2
# 전체 모형 (p=2)
epsilon_2 = bartlett_epsilon_exp(X)
print(f"ε_2 (full model, p=2): {epsilon_2:.4f}")
# 귀무 모형 (q=1, intercept only)
X0 = np.ones((len(df), 1))
epsilon_1 = bartlett_epsilon_exp(X0)
print(f"ε_1 (null, q=1): {epsilon_1:.4f}")
print(f"이론값 1/(6n) = {1/(6*len(df)):.4f}")
# Bartlett 인수
b_21 = (epsilon_2 - epsilon_1) / (2 - 1)
print(f"b_21 = (ε_2 - ε_1)/(p-q) = {b_21:.4f}")
# 조정된 LRT
LRT_adj = LRT / (1 + b_21)
print(f"\n원 LRT = {LRT:.3f}, 조정 LRT = {LRT_adj:.3f}")
print(f"원 p = {1-stats.chi2.cdf(LRT, 1):.4f}, "
f"조정 p = {1-stats.chi2.cdf(LRT_adj, 1):.4f}")기대: 원 LRT ≈ 6.83, 조정 LRT ≈ 6.74, \(p\) 변화 미미.
8.4 모형 적합도 검정의 큰 Impact
# 포화 모형 (p=17) — 각 관측치에 자유 μ
# 이탈도 = 2*(log-likelihood of saturated - full model)
# Gamma exp 에서 saturated μ_i = Y_i
mu_hat_full = m_full.fittedvalues
log_lik_full = -np.sum(df['weeks']/mu_hat_full + np.log(mu_hat_full))
log_lik_sat = -np.sum(df['weeks']/df['weeks'] + np.log(df['weeks'])) # = -n - Σ log(Y)
LRT_gof = 2 * (log_lik_sat - log_lik_full)
print(f"\n모형 적합도 LRT (15 df): {LRT_gof:.3f}")
# ε_17 — 포화 모형 (Q = I 에 해당, 단 X 가 n×n 동정 행렬일 경우)
# 실제로는 ε_n = n/6 (McCullagh-Nelder)
epsilon_17 = 17 / 6
print(f"ε_17 = n/6 = {epsilon_17:.3f}")
# b_{17,2}
b_17_2 = (epsilon_17 - epsilon_2) / (17 - 2)
print(f"b_{{17,2}} = {b_17_2:.4f}")
# 조정된 통계량
LRT_gof_adj = LRT_gof / (1 + b_17_2)
print(f"\n원 적합도 LRT = {LRT_gof:.2f}, 조정 = {LRT_gof_adj:.2f}")
print(f"원 p = {1-stats.chi2.cdf(LRT_gof, 15):.3f}, "
f"조정 p = {1-stats.chi2.cdf(LRT_gof_adj, 15):.3f}")기대: 원 LRT ≈ 19.46, 조정 ≈ 16.4, \(p\) 19% → 36% (대략 2 배 증가).
8.5 Bootstrap-Bartlett 검증
B = 500
LRT_boot = []
np.random.seed(42)
mu_null = m_null.fittedvalues
for b in range(B):
# 지수 분포 시뮬레이션
y_boot = np.random.exponential(scale=mu_null, size=len(df))
try:
m_f_boot = sm.GLM(y_boot, X,
family=sm.families.Gamma(link=sm.families.links.log())).fit(scale=1, disp=0)
m_n_boot = sm.GLM(y_boot, np.ones(len(df)),
family=sm.families.Gamma(link=sm.families.links.log())).fit(scale=1, disp=0)
LRT_boot.append(2*(m_f_boot.llf - m_n_boot.llf))
except:
continue
mean_LRT_boot = np.mean(LRT_boot)
b_boot = mean_LRT_boot / 1 - 1 # p - q = 1
print(f"\nBootstrap-Bartlett:")
print(f"Mean LRT in bootstrap: {mean_LRT_boot:.3f}")
print(f"b_p (bootstrap): {b_boot:.4f}")
print(f"b_21 (theoretical): {b_21:.4f}")기대: Bootstrap \(b\) 와 이론 \(b\) 가 비슷. 차이가 크면 점근 근사 문제 있음.
9 요약 — §15.3 의 네 가지 교훈
9.1 교훈 1 — 단일 스칼라로 모든 cumulant 보정
Bartlett 의 천재성: \((1 + b_p)^r\) 의 곱셈적 구조로 한 번의 scale 교정이 모든 차수 동시 교정. Edgeworth 와 달리 shape 수정 불필요.
9.2 교훈 2 — GLM 의 Mixed Cumulant 대칭성
GLM 의 특수 성질 덕분에 6 개 불변 스칼라가 매트릭스 연산 으로 계산됨. 일반 모형에서는 훨씬 복잡. Canonical link 에서 3 스칼라 로 추가 단순화.
9.3 교훈 3 — 보정 효과는 자유도에 비례
Feigl-Zelen 예제가 보여 주듯: - 1 df 검정: Bartlett 영향 미미. - 15 df 적합도 검정: \(p\)-value 2 배 변화.
실무 권고: 자유도 큰 검정에서 Bartlett 필수. 단일 계수 검정은 생략 가능.
9.4 교훈 4 — Lattice Case 의 한계
이산 반응에서 Bartlett 은 \(O(n^{-1/2})\) 이하로 개선 못 함. Continuous 분포에서만 full \(O(n^{-2})\) 효과.
9.5 한 줄 정리
“LRT 의 \(\chi^2\) 근사에는 \(O(n^{-1})\) scale bias 가 있고, Bartlett 은 단일 스칼라 로 이것을 \(O(n^{-2})\) 로 교정 한다. GLM 은 이 계산을 매트릭스 연산 으로 감당한다.”
10 관련 주제
선행 지식
- Further Topics — Ch.15 개관
- Bias Adjustment — §15.2
- Score Tests for Extra Parameters (McCullagh §12.3) — LRT·Score·Wald 비교
- GLM 적합도 측정 — Deviance (McCullagh §2.3)
- Exponential · Gamma 분포 (McCullagh §8.2)
관련 개념
- Model Checking — 개관 (McCullagh Ch.12) — 이탈도 검정의 진단 맥락
- 점근 검정 이론 (Casella-Berger Ch.10)
참고 문헌
- Bartlett, M. S. (1937). “Properties of sufficiency and statistical tests.” Proc. R. Soc. Lond. A 160: 268-282. 원 논문.
- Bartlett, M. S. (1954). “A note on the multiplying factors for various \(\chi^2\) approximations.” JRSS B 16: 296-298.
- McCullagh, P. & Cox, D. R. (1986). “Invariants and likelihood ratio statistics.” Ann. Statist. 14: 1419-1430. 식 (15.9) 의 출처.
- Cordeiro, G. M. (1983, 1987). Bartlett corrections in GLM 논문 시리즈.
- Lawley, D. N. (1956). “A general method for approximating to the distribution of likelihood ratio criteria.” Biometrika 43: 295-303. 다변량 분석으로 확장.
현대적 대안
- Durbin, J. (1980). “Approximations for densities of sufficient estimators.” Biometrika 67: 311-333. Saddlepoint.
- Shao, J. (2003). Mathematical Statistics. 점근 이론의 현대적 정리.
후속 주제