Likelihood Functions for Log-linear Models

포아송 로그우도·이탈도·과산포 메커니즘·점근이론 (McCullagh & Nelder §6.2)

로그선형 모형의 우도 구조를 정리한다. 포아송 로그우도 \(\sum(y\log\mu - \mu)\) 의 기하, 이탈도 \(D = 2\sum y\log(y/\hat\mu)\) 의 Kullback–Leibler 해석과 \(\mu^{1/3}\) 근사, 과산포의 세 메커니즘(랜덤 구간, 클러스터 포아송, 감마 혼합 → 음이항), quasi-Poisson 와 NB 의 실무 선택, IRLS 점근 이론을 수식과 직관으로 함께 전개한다.

Statistics
GLM
저자

Kwangmin Kim

공개

2026년 04월 15일

1 왜 우도를 따로 다루는가

§6.1 의 개관에서 “\(\log \mu = \boldsymbol{\beta}^\top \mathbf{x}\)” 가 로그선형 모형의 뼈대임을 봤다. 이제 이 뼈대를 데이터로 붙잡아 추정 하는 기계 — 우도함수 — 를 전개할 차례이다.

§6.2 의 중심 질문들:

  • 포아송 로그우도는 어떻게 생겼으며 왜 그 모양인가
  • 이탈도 \(D = 2\sum y\log(y/\hat\mu)\) 는 무엇을 측정하는가
  • 과산포는 어떤 메커니즘에서 발생하며, 그에 따라 어떤 분포를 쓰는가
  • 모수가 많은 희소 분할표에서도 점근 이론이 유효한가

이 포스트는 §6.2.2 (로그우도), §6.2.3 (과산포의 세 메커니즘과 음이항), §6.2.4 (점근 이론) 를 순서대로 따라가며 수식 뒤의 직관을 정리한다.

2 §6.2.2 포아송 로그우도

2.1 단일 관측의 기여

한 관측 \(y\) 의 포아송 로그우도(상수 제외):

\[ \ell(\mu; y) = y \log \mu - \mu. \]

모수 무관 상수 \(-\log y!\) 는 생략 — 모수 추정·검정에 영향을 주지 않는다.

2.2 세 가지 척도에서의 형태

\(\ell(\mu)\) 를 세 가지 축에서 그려보면 이 함수의 본질을 알 수 있다 (교재 Fig. 6.2, \(y = 1\)).

축 특징
\(\mu\) (원래 척도) 강한 비대칭 — \(\mu\) 가 \(y\) 보다 작을 때보다 클 때 더 빠르게 떨어짐
\(\log \mu\) (정준 모수) 거의 대칭이지만 여전히 약간 왜곡
\(\mu^{1/3}\) (큐브 루트) 거의 완벽한 이차함수

왜 \(\mu^{1/3}\) 가 특별한가: Anscombe (1953) 변환 \(Y^{2/3}\) 이 포아송을 가장 대칭화하는 것과 같은 이유이다. 가까운 최댓값 주변에서

\[ y \log \mu - \mu \approx y \log y - y - \frac{9 y^{1/3}(\mu^{1/3} - y^{1/3})^2}{2} \]

로 근사된다 (Exercise 6.2). 이것은 포아송 로그우도가 \(\mu^{1/3}\) 척도에서 이차함수로 잘 근사된다는 강력한 주장.

실무적 함의: 수치 최적화에서 \(\mu^{1/3}\) 파라미터화를 쓰면 뉴턴-라프슨이 빠르게 수렴한다. 비대칭이 줄어 신뢰구간도 \(\mu^{1/3}\) 척도에서 구하면 정규 근사가 훨씬 좋다.

2.3 신호편차 근사

위 근사의 신호편차(signed deviance) 제곱근

\[ g(y) = 3 y^{1/6}(y^{1/3} - \mu^{1/3}) \]

가 \(\mu \to \infty\) 에서 표준 정규에 근사. §6.1 의 tail 확률 근사 공식의 주요 항이 이것이다.

3 전체 로그우도와 이탈도

독립 관측 벡터 \(\mathbf{y} = (y_1, \ldots, y_n)\) 의 로그우도:

\[ \ell(\boldsymbol{\mu}; \mathbf{y}) = \sum_i (y_i \log \mu_i - \mu_i). \tag{6.3} \]

3.1 이탈도 — 세 가지 표현

포화 모형의 최대 로그우도는 \(\hat\mu_i = y_i\) 에서 달성. 이탈도는

\[ \begin{aligned} D(\mathbf{y}; \hat{\boldsymbol{\mu}}) &= 2\ell(\mathbf{y}; \mathbf{y}) - 2\ell(\hat{\boldsymbol{\mu}}; \mathbf{y}) \\ &= 2\sum_i \{y_i \log(y_i / \hat\mu_i) - (y_i - \hat\mu_i)\}. \end{aligned} \]

상수항 포함 모형에서의 단순화: 로그선형 모형에 절편이 있으면 \(\sum(y_i - \hat\mu_i) = 0\) (스코어 방정식의 한 귀결). 따라서

\[ D(\mathbf{y}; \hat{\boldsymbol{\mu}}) = 2\sum_i y_i \log(y_i / \hat\mu_i). \]

이것이 가장 자주 보는 형태.

3.2 \(\mu^{1/3}\) 근사

포아송 로그우도의 큐브 루트 이차 근사로부터

\[ D(\mathbf{y}; \boldsymbol{\mu}) \approx 9\sum_i y_i^{1/3}(y_i^{1/3} - \mu_i^{1/3})^2. \]

의미: 이탈도가 \(\mu^{1/3}\) 척도에서 가중 제곱합으로 근사된다. 가중치는 \(y_i^{1/3}\). 포아송 GLM 의 이탈도 잔차가 \(\mu^{1/3}\) 변환 공간에서 대략 정규처럼 행동하는 이론적 근거.

3.3 Pearson \(X^2\) 근사

테일러 전개 \((y - \mu)/\mu\) 로부터

\[ D(\mathbf{y}; \boldsymbol{\mu}) \approx \sum_i \frac{(y_i - \mu_i)^2}{\mu_i} = X^2. \]

이것이 Pearson (1900) 의 \(X^2\) 통계량. 큐브 루트 근사보다 정확도가 낮지만 해석이 직관적.

세 가지 이탈도 표현의 계층
표현 정확도 용도
\(2\sum y_i \log(y_i/\hat\mu_i)\) 정확 (정의) 표준 이탈도, LR 검정
\(9\sum y_i^{1/3}(y_i^{1/3} - \hat\mu_i^{1/3})^2\) \(O(\mu^{-3/2})\) 잔차 분석·시각화
\(\sum (y_i - \hat\mu_i)^2 / \hat\mu_i\) (Pearson \(X^2\)) \(O(\mu^{-1/2})\) 과산포 진단

\(\mu\) 가 작을수록 이탈도와 Pearson 이 괴리된다 — 희소 데이터에서는 선택이 중요.

3.4 이탈도의 KL 해석

\(D/2\) 는 경험 분포 \(\tilde{\pi}_i = y_i/y_\bullet\) 와 적합 분포 \(\hat\pi_i\) 의 KL 발산. 포아송 우도 관점에서도 이 해석이 유지된다 — 포아송의 총합이 고정되면 자동으로 다항이 되고, 조건부 KL 발산이 이탈도와 일치한다.

4 §6.2.3 과산포 — 세 가지 발생 메커니즘

4.1 메커니즘 1: 랜덤 구간 길이

관측 구간 길이 \(T\) 가 확정이 아니라 랜덤이면, 총 관측 수 \(Y\) 는 조건부 포아송이지만 주변 분포에서는 포아송이 아니다.

\[ Y \mid T \sim \text{Poisson}(\lambda T), \qquad T \text{ 랜덤}. \]

전체 분산 분해:

\[ \mathrm{Var}(Y) = \mathrm{E}[\mathrm{Var}(Y \mid T)] + \mathrm{Var}[\mathrm{E}(Y \mid T)] = \lambda \mathrm{E}(T) + \lambda^2 \mathrm{Var}(T). \]

\(\mathrm{Var}(T) > 0\) 이면 \(\mathrm{Var}(Y) > \mathrm{E}(Y)\) — 과산포.

실제 예: 동물 추적에서 배터리 수명이 다르거나, 카운터 작동 시간이 관측마다 미묘하게 다른 경우.

4.2 메커니즘 2: 클러스터 포아송

한 사건이 여러 개의 하위 관측을 만드는 경우.

\[ Y = Z_1 + Z_2 + \cdots + Z_N, \qquad N \sim \text{Poisson}, \; Z_j \text{ i.i.d., } N \perp Z. \]

모멘트 계산:

\[ \mathrm{E}(Y) = \mathrm{E}(N)\mathrm{E}(Z), \qquad \mathrm{Var}(Y) = \mathrm{E}(N)\mathrm{E}(Z^2) = \mathrm{E}(N)\{\mathrm{Var}(Z) + (\mathrm{E}Z)^2\}. \]

따라서 과산포 조건 은 \(\mathrm{E}(Z^2) > \mathrm{E}(Z)\).

실제 예:

  • 교통사고: 한 사고가 여러 차량을 포함 (\(N\) = 사고 수, \(Z\) = 사고당 피해 차량 수)
  • 동물 번식: 한 어미가 여러 새끼를 낳음 (\(N\) = 번식 수, \(Z\) = 한 번에 낳는 새끼 수)
  • 광고 클릭: 한 봇이 여러 번 클릭 (\(N\) = 봇 수, \(Z\) = 봇당 클릭)

4.3 메커니즘 3: 감마 혼합 → 음이항

개체마다 발생률이 다른 경우. 관측 수 \(Y\) 가 조건부 포아송이고 조건부 평균 \(Z\) 가 감마분포:

\[ Y \mid Z \sim \text{Poisson}(Z), \qquad Z \sim \text{Gamma}(\text{mean} = \mu,\; \text{index} = \phi\mu). \]

즉 \(\mathrm{E}(Z) = \mu\), \(\mathrm{Var}(Z) = \mu/\phi\).

주변 분포 는 음이항:

\[ \Pr(Y = y;\; \mu, \phi) = \frac{\Gamma(y + \phi\mu)\, \phi^{\phi\mu}}{y!\, \Gamma(\phi\mu) (1+\phi)^{y+\phi\mu}}, \quad y = 0, 1, 2, \ldots \]

모멘트:

\[ \mathrm{E}(Y) = \mu, \qquad \mathrm{Var}(Y) = \mu(1 + \phi)/\phi = \mu \cdot \sigma^2 \]

여기서 \(\sigma^2 = (1 + \phi)/\phi = 1 + 1/\phi\) 는 상수 산포 인자. 이것이 식 (6.1) 의 선형 분산 \(\mathrm{Var}(Y) = \sigma^2 \mu\) 과 정확히 일치.

직관 — 왜 감마 혼합인가:

감마는 포아송의 켤레 사전분포. 이것이 수학적으로 편리한 이유도 있지만, 더 중요한 것은 감마의 “\(\mathrm{Var}(Z) = \mu/\phi\)” 구조가 포아송 자체와 같은 평균-분산 관계 를 재현해 전체 모형이 여전히 linear variance 관계를 유지한다는 점이다.

4.4 대안 감마 혼합 — 이차 분산 음이항

지수를 \(\nu\) 로 고정(평균 \(\mu\) 에 무관)하면:

\[ \mathrm{Var}(Y) = \mu + \mu^2/\nu. \]

분산이 평균에 이차 의존. 이것이 실무에서 “NB2” 모형으로 불리는 버전. statsmodels 의 NegativeBinomial 기본값.

혼합 방식 분산 구조 이름
감마 index = \(\phi\mu\) (mean-dependent) \(\sigma^2 \mu\) (선형) NB1
감마 index = \(\nu\) (고정) \(\mu + \mu^2/\nu\) (이차) NB2

4.5 NB vs Quasi-Poisson — McCullagh 의 실무 권고

두 접근의 차이:

항목 NB Quasi-Poisson
분포 가정 완전 지정 평균-분산 관계만
추정 MLE (분포 명시) 추정함수 (pseudo-score)
점추정 \(\hat{\boldsymbol{\beta}}\) NB MLE Poisson MLE 와 동일
표준오차 NB 정보행렬 \(\sigma^2 \times\) Poisson 정보행렬
분산 가정 틀릴 때 편향 가능 대체로 강건

McCullagh 의 결론: “두 접근의 \(\hat{\boldsymbol{\beta}}\) 차이는 \(O_p(\phi^{-2})\) 로, 과산포가 심하지 않으면 무시 가능. 따라서 분포를 특정하지 않는 quasi-Poisson 이 실무적으로 더 안전” (§6.2.3).

4.6 Quasi-Poisson 보정의 실행

포아송 MLE 로 \(\hat{\boldsymbol{\beta}}\) 를 구한 뒤:

\[ \hat\sigma^2 = X^2 / (n - p) \]

로 산포를 추정하고 공분산을 \(\hat\sigma^2 \times (\mathbf{X}^\top \mathbf{W} \mathbf{X})^{-1}\) 로 부풀린다. 점추정은 불변. §4.5·§5.5 에서 본 패턴이 포아송에도 그대로 반복.

5 §6.2.4 점근 이론

5.1 일치성과 정규성의 조건

\(\hat{\boldsymbol{\beta}}\) 의 일치성과 점근 정규성은 정보행렬 \(\mathbf{i}_\beta\) 의 고유값이 무한대로 증가 하면 성립. 실무에서 이 조건은 다음 두 경우에 만족된다.

  1. \(p\) 고정, \(n \to \infty\) — 관측 수가 증가
  2. \(n, p\) 고정, 각 \(\mu_i \to \infty\) — 셀마다 기대도수가 커짐

실무적 의미: (a) 관측 많으면 OK, (b) 관측 적어도 셀 도수가 크면 OK. (a)·(b) 모두 실패하는 희소 분할표(작은 \(n\) + 작은 \(\mu_i\)) 에서만 점근이 나빠진다.

5.2 점근 공분산

\[ \mathrm{Cov}(\hat{\boldsymbol{\beta}}) \approx \sigma^2\, \mathbf{i}_\beta^{-1}, \quad \mathbf{i}_\beta = \mathbf{X}^\top \mathrm{diag}(\boldsymbol{\mu}) \mathbf{X}. \]

가중치가 \(\mu_i\) 자체라는 점이 특징이다 — 이항의 \(\mu_i(1 - \pi_i)\) 와 비교. 포아송은 분산 = 평균이라 가중치가 단순해진다.

5.3 산포 추정 — 식 (6.4)

\[ \hat\sigma^2 = \frac{X^2}{n - p} = \frac{1}{n-p}\sum_i \frac{(y_i - \hat\mu_i)^2}{\hat\mu_i}. \tag{6.4} \]

왜 이탈도 \(D/(n-p)\) 가 아닌가: §5.5 에서와 같은 이유. 희소 데이터에서 \(D\) 기반 추정량은 편향되며 Pearson 기반이 더 견고하다 (희소·비희소 모두 일치 추정량).

5.4 유효 자유도 — t 근사

산포를 추정하면 \(\hat{\boldsymbol{\beta}}\) 의 분포가 정확한 정규가 아니라 \(t\) 분포 근사가 된다.

\[ f = \frac{n - p}{1 + \frac{1}{2}\tilde\gamma_2} \]

\(\tilde\gamma_2\) 는 표준화 첨도(kurtosis). 이것은 (4.14)·(5.13) 의 로그선형 버전이다. 대부분의 실무 상황(평균 \(\mu_i \ge 1\) 정도)에서 \(f \approx n - p\) 로 두고 \(t_{n-p}\) 또는 정규 근사를 써도 무방.

작은 평균(\(\mu_i < 1\)) 영역: 이 보정이 실제로 필요해진다. 희소 분할표에서 보수적 추론을 하려면 \(t_f\) 분위수를 쓴다.

6 응용 — 어떻게 우도가 작동하는지

6.1 선형 대조의 t 검정

Fisher (1949) 의 결핵균 생물학적 검정(§6.3.1) 에서 4처리의 라틴 정방 설계. 로그 척도에서 고용량 vs 저용량 대조:

\[ \hat{\beta} = 0.2095, \quad \text{SE} = 0.0124, \quad \hat\sigma^2 = 0.201. \]

해석: “용량 2배 → 반응 \(\exp(0.2095) = 1.233\) 배 증가”. 23.3% 증가.

Weybridge vs Standard 대조:

\[ \hat{\beta} = 0.0026, \quad \text{SE} = 0.0123. \]

해석: “같은 용량일 때 Weybridge 는 Standard 대비 거의 차이 없음 (\(e^{0.003} \approx 1.003\)). 통계적으로 유의하지 않다.”

상대 효능: \(2^{0.2121/0.2095} = 2.017\). “Weybridge 는 Standard 의 약 두 배 효능” (동일 반응을 위해 절반 용량 필요).

7 코드 예시

7.1 Step 1: 순수 Python — 포아송 로그우도의 세 척도 시각화

import numpy as np
import matplotlib.pyplot as plt

y = 1
mu = np.linspace(0.2, 4, 500)
loglik = y * np.log(mu) - mu
logmu = np.log(mu)
mu13 = mu ** (1/3)

# 이차 근사 on mu^(1/3)
approx = y * np.log(y) - y - 9 / 2 * y ** (1/3) * (mu13 - y ** (1/3)) ** 2

fig, axes = plt.subplots(1, 3, figsize=(12, 3.5))
axes[0].plot(mu, loglik); axes[0].set(title="l(mu) vs mu", xlabel="mu")
axes[1].plot(logmu, loglik); axes[1].set(title="l(mu) vs log mu", xlabel="log mu")
axes[2].plot(mu13, loglik, label="exact")
axes[2].plot(mu13, approx, "--", label="quadratic on mu^(1/3)")
axes[2].set(title="l(mu) vs mu^(1/3)", xlabel="mu^(1/3)")
axes[2].legend()
plt.tight_layout()

세 번째 그림에서 점선(이차 근사) 과 실선(정확) 이 거의 겹침. 큐브 루트 척도에서 이차 근사가 매우 정확함이 시각적으로 확인된다.

7.2 Step 2: 세 가지 이탈도 표현의 수치 비교

import numpy as np

rng = np.random.default_rng(0)
n = 200
x = rng.normal(size=n)
mu_true = np.exp(1.0 + 0.5 * x)
y = rng.poisson(mu_true)

# 모형 적합 (포아송 MLE — 간단한 직접 최적화)
import statsmodels.api as sm
X = sm.add_constant(x)
fit = sm.GLM(y, X, family=sm.families.Poisson()).fit()
mu_hat = fit.predict(X)

def deviance_exact(y, mu_hat):
    mask = y > 0
    term = np.zeros_like(y, dtype=float)
    term[mask] = y[mask] * np.log(y[mask] / mu_hat[mask]) - (y[mask] - mu_hat[mask])
    term[~mask] = -(y[~mask] - mu_hat[~mask])  # y log y = 0 at y=0
    return 2 * np.sum(term)

def deviance_cuberoot(y, mu_hat):
    return 9 * np.sum(y ** (1/3) * (y ** (1/3) - mu_hat ** (1/3)) ** 2)

def pearson_chi2(y, mu_hat):
    return np.sum((y - mu_hat) ** 2 / mu_hat)

D_exact = deviance_exact(y, mu_hat)
D_cube = deviance_cuberoot(y, mu_hat)
X2 = pearson_chi2(y, mu_hat)

print(f"Exact deviance      = {D_exact:.3f}  (statsmodels = {fit.deviance:.3f})")
print(f"Cube-root approx    = {D_cube:.3f}  (오차 {abs(D_cube-D_exact)/D_exact*100:.1f}%)")
print(f"Pearson X^2         = {X2:.3f}  (오차 {abs(X2-D_exact)/D_exact*100:.1f}%)")

큐브 루트 근사는 정확한 이탈도와 1% 이내로 일치하지만, Pearson \(X^2\) 는 몇 % 차이를 보인다. \(\mu\) 가 작은 영역에서 이 차이가 커진다.

7.3 Step 3: Poisson vs Quasi-Poisson vs NB 비교

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

rng = np.random.default_rng(1)
n = 300
x = rng.normal(size=n)

# 감마 혼합으로 과산포 데이터 생성 → NB2
phi = 2.0  # var(Z) = mu/phi → sigma^2 = 1 + 1/phi = 1.5
mu_true = np.exp(1.0 + 0.4 * x)
Z = rng.gamma(shape=phi * mu_true, scale=1 / phi)  # E(Z)=mu, Var(Z)=mu/phi
y = rng.poisson(Z)

df = pd.DataFrame({"x": x, "y": y})

fit_p = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit()
fit_q = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit(scale="X2")
fit_nb = smf.glm("y ~ x", data=df, family=sm.families.NegativeBinomial(alpha=1/phi)).fit()

for name, f in [("Poisson", fit_p), ("Quasi-Poisson", fit_q), ("NegBin", fit_nb)]:
    print(f"=== {name} ===")
    print(f"  beta = {f.params.values.round(3)}")
    print(f"  SE   = {f.bse.values.round(3)}")
print(f"\nQuasi scale = {fit_q.scale:.3f} (이론: 1 + 1/phi = {1 + 1/phi})")

관찰: 세 모형 모두 점추정 \(\hat{\boldsymbol{\beta}}\) 은 거의 동일. SE 만 다름 — Poisson 이 가장 작고 (과산포 무시), Quasi와 NB 는 거의 같은 부풀려진 값을 준다. 이것이 §6.2.3 의 “두 추정이 \(O_p(\phi^{-2})\) 로 차이난다” 의 수치적 확인.

7.4 R 대응

# 포아송, 과산포 보정, 음이항
fit_p <- glm(y ~ x, data = df, family = poisson())
fit_q <- glm(y ~ x, data = df, family = quasipoisson())
library(MASS)
fit_nb <- glm.nb(y ~ x, data = df)

lapply(list(poisson = fit_p, quasi = fit_q, negbin = fit_nb), summary)

8 자주 걸리는 함정

함정 증상 처방
\(\mu^{1/3}\) 근사가 모든 영역에서 정확하다고 가정 극히 작은 \(\mu\) 에서 오차 정확 이탈도 직접 계산
이탈도 \(D/(n-p)\) 로 \(\sigma^2\) 추정 희소 데이터에서 편향 Pearson \(X^2/(n-p)\) 사용
과산포를 NB 로만 대응 분포 오지정 위험 Quasi-Poisson 먼저
가우시안 근사로 모든 \(\beta\) CI 구성 희소·작은 \(\mu\) 에서 비대칭 \(t_f\) 근사, \(\mu^{1/3}\) 척도 CI
큐브 루트 이차 근사의 단조성 깨짐 \(\mu > 38\) 영역 해당 영역 확률 무시할 수 있음 (교재 지적)
이탈도를 절대 적합도 통계로 사용 희소에서 \(\chi^2\) 근사 나쁨 중첩 모형 비교(LR)에만 사용
NB2 와 NB1 혼동 분산 구조 \(\mu + \mu^2/\nu\) vs \(\sigma^2 \mu\) 데이터의 \(X^2/\text{df}\) 선형성 확인
offset 누락으로 “rate → count” 오류 계수에 노출 효과 섞임 \(\log(\text{exposure})\) offset 필수

9 관련 주제

선행 지식

후속 주제 (placeholder)

관련 개념

10 참고문헌

  • McCullagh, P. & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.), §6.2. Chapman & Hall.
  • Fisher, R. A. (1949). A biological assay of tuberculins. Biometrics, 5, 300–316.
  • Anscombe, F. J. (1953). Contribution to the discussion of H. Hotelling’s paper. JRSS B, 15, 229–230.
  • Plackett, R. L. (1981). The Analysis of Categorical Data (2nd ed.). Griffin.
  • Pearson, K. (1900). On the criterion that a given system of deviations… Philosophical Magazine, 50, 157–175.
  • Cameron, A. C. & Trivedi, P. K. (2013). Regression Analysis of Count Data (2nd ed.). Cambridge.
  • Hilbe, J. M. (2011). Negative Binomial Regression (2nd ed.). Cambridge.

Subscribe

Enjoy this blog? Get notified of new posts by email: