Kwangmin Kim - Parameters in the Variance Function — 음이항의 \(k\) 와 반올림 오차 분산

1 서론 — 분산함수 속 모수가 왜 까다로운가

표준 GLM 의 분산 구조는

\[ \mathrm{var}(Y_i) = \phi \, V(\mu_i) \]

에서 두 역할 을 엄격히 분리한다.

\(\phi\) — 산포 모수(dispersion): 전역 상수. \(\hat\beta\) 의 정규 방정식에서 소거 되므로 \(\phi\) 를 몰라도 \(\hat\beta\) 를 구할 수 있다. 추정은 사후에 Pearson \(X^2/(n-p)\) 등으로 한다
\(V(\mu)\) — 분산함수 형태: 분포족이 정하는 고정 함수. 포아송은 \(V = \mu\), 감마는 \(V = \mu^2\), 이항은 \(V = \pi(1-\pi)/m\)

이 관습이 깨지는 두 상황이 있다.

음이항(Negative Binomial): \(V(\mu; k) = \mu + \mu^2/k\). 여기서 \(k\) 는 \(\phi\) 처럼 \(\hat\beta\) 에서 분리되지 않는다. \(k\) 가 바뀌면 가중치 \(W = \text{diag}(1/V(\mu_i))\) 가 바뀌어 모든 \(\hat\beta_j\) 가 움직인다.
반올림 오차: \(V(\mu; \tau^2, \sigma^2) = \tau^2 + \sigma^2 \mu^2\). 두 분산 성분의 비율 \(\sigma^2/\tau^2\) 이 미지. 이 비율이 바뀌면 작은 \(\mu\) 와 큰 \(\mu\) 의 상대 가중치 자체가 변한다.

두 경우 모두 산포가 아닌 구조 모수 가 분산함수에 들어간다. McCullagh & Nelder (1989, §11.2) 는 이들을 GLM 프레임워크 안에서 처리하는 체계를 제시한다.

1.1 직관: \(\phi\) 와 \(k\) 의 기하학적 차이

\(\mathrm{var}(Y_i) = \phi \cdot V(\mu_i)\) 를 그래프로 그리면, \(\phi\) 는 분산 곡선 전체를 위아래로 평행이동 시킨다 (수직 스케일). 모든 관측치의 상대적 가중치 \(1/V(\mu_i)\) 는 불변이다. 따라서 \(\phi\) 는 \(\hat\beta\) 에 영향을 주지 않는다.

반면 \(V(\mu; k) = \mu + \mu^2/k\) 의 \(k\) 는 분산 곡선의 모양 자체를 바꾼다. \(k\) 가 작으면 큰 \(\mu\) 의 분산이 이차적으로 폭발 — 큰 관측치의 가중치가 크게 줄어든다. \(k\) 가 크면 포아송에 가까워져 \(V \propto \mu\). 상대 가중치가 \(k\) 에 의존 하므로 \(\hat\beta\) 가 \(k\) 와 얽힌다.

2 음이항 — Poisson-Gamma 혼합으로부터의 유도

2.1 생성 메커니즘

포아송 모형의 한계: 실제 count 데이터는 종종 \(\mathrm{var}(Y) > E(Y)\) — 과산포. 포아송 평균 \(\Lambda_i\) 자체가 개체별로 무작위로 변동한다고 가정하자.

\[ \Lambda_i \sim \text{Gamma}(k, \alpha), \qquad Y_i \mid \Lambda_i \sim \text{Poisson}(\Lambda_i) \]

감마 분포의 파라미터화: 형태 \(k\), 척도 \(\alpha\), 따라서 \(E(\Lambda_i) = k\alpha\), \(\mathrm{var}(\Lambda_i) = k\alpha^2\).

주변분포를 구하면:

\[ \Pr(Y = y) = \int_0^\infty \frac{e^{-\lambda} \lambda^y}{y!} \cdot \frac{\lambda^{k-1} e^{-\lambda/\alpha}}{\alpha^k \Gamma(k)} d\lambda \]

이 적분은 감마 함수의 정의로 정리되어

\[ \Pr(Y = y; \alpha, k) = \frac{(y+k-1)!}{y!(k-1)!} \frac{\alpha^y}{(1+\alpha)^{y+k}}, \quad y = 0, 1, 2, \ldots \tag{11.2} \]

2.2 평균과 분산 — 조건부 기댓값 공식

\[ E(Y) = E\!\left[ E(Y | \Lambda) \right] = E(\Lambda) = k\alpha = \mu \]

\[ \mathrm{var}(Y) = E[\mathrm{var}(Y | \Lambda)] + \mathrm{var}[E(Y | \Lambda)] = E(\Lambda) + \mathrm{var}(\Lambda) = k\alpha + k\alpha^2 \]

\(\mu = k\alpha\) 로 치환하면 \(\alpha = \mu/k\), 따라서

\[ \mathrm{var}(Y) = \mu + \mu^2/k \tag{11.3} \]

2.3 직관: 과산포의 두 원천 분해

\(V(\mu) = \mu + \mu^2/k\) 의 두 항은 서로 다른 출처를 가진다.

\(\mu\) 항 — 포아송 자체의 분산. \(\Lambda\) 가 고정되었을 때의 conditional variance
\(\mu^2/k\) 항 — 개체 간 \(\Lambda\) 의 변동성. \(k\) 가 혼합의 “집중도” 를 결정

\(k \to \infty\) 이면 감마 분포가 \(\mu\) 에 집중 (분산 \(\to 0\)) → 순수 포아송. \(k\) 가 작으면 \(\Lambda\) 가 크게 들쭉날쭉해 마지널 분산이 이차적으로 커진다. \(k\) 는 “과산포 강도의 역수” 로 읽으면 직관이 맞다.

2.4 로그우도와 canonical link 의 문제

(11.2) 를 다시 쓰면:

\[ \ell = y \log\!\left(\frac{\alpha}{1+\alpha}\right) - k \log(1 + \alpha) + c(y, k) \]

여기서 \(c(y, k)\) 는 \(\alpha\) 에 의존하지 않는 항. 고정 \(k\) 에서 이 로그우도는 canonical link

\[ \eta = \log\!\left(\frac{\alpha}{1+\alpha}\right) = \log\!\left(\frac{\mu}{\mu+k}\right) \tag{11.4} \]

와 함께 지수족 꼴이다. 그러나 (11.4) 의 링크 자체가 \(k\) 에 의존한다. \(k\) 를 고정하지 않으면 링크 함수가 변동하여 IRLS 의 설계행렬이 매 반복마다 달라진다.

2.5 실무적 선택: 로그 링크 + 음이항 분산함수

실전에서는 canonical link 를 버리고 로그 링크 를 쓴다.

\[ \log \mu_i = \mathbf x_i^T \boldsymbol\beta, \qquad V(\mu_i; k) = \mu_i + \mu_i^2/k \]

장점 1: 링크가 \(k\) 와 분리되어 해석·알고리즘이 단순
장점 2: 포아송 회귀와 동일한 링크라 계수 해석 일관 — \(\beta_j\) 는 로그 rate ratio
장점 3: 로그 링크는 모든 계수 조합에서 \(\mu > 0\) 보장

이 결합이 glm.nb (R), statsmodels.NegativeBinomial (Python) 의 기본 동작이다.

3 \(k\) 의 추정 — 네 가지 경로

3.1 ML 추정: digamma 방정식

\(\boldsymbol\beta\) 를 고정한 상태에서 \(k\) 의 스코어:

\[ \frac{\partial \ell}{\partial k} = \sum_i \left[ \psi(y_i + k) - \psi(k) + \log\!\left(\frac{k}{\mu_i + k}\right) + \frac{y_i - \mu_i}{\mu_i + k} \right] \]

여기서 \(\psi(\cdot)\) 는 digamma 함수. 이 식 \(= 0\) 은 비선형 방정식 — 해석해가 없다. Newton-Raphson 또는 Brent 법으로 푼다.

문제점: digamma 함수 호출로 계산이 느리고, \(\beta\) 와 \(k\) 가 얽혀 이중 반복(외부 \(k\), 내부 IRLS)이 필요.

3.2 Moment matching — Anscombe (1949)

\(\mathrm{var}(Y) = \mu + \mu^2/k\) 를 표본분산과 맞춘다. 단일 표본(공변량 없음)에서

\[ \hat k_{\text{mom}} = \frac{\bar y^2}{s^2 - \bar y} \]

이 단순하고 설명가능하지만, 구조화된 공변량 데이터에서는 일반화가 애매하다.

3.3 Mean deviance = 1

이항·포아송 GLM 에서 \(\phi\) 추정에 쓰는 관습을 \(k\) 에 적용:

\[ \hat k \text{ s.t. } \frac{D(\hat\beta(k); k)}{n - p} = 1 \]

\(k\) 를 격자로 움직이며 각 \(k\) 에서 GLM 을 적합해 이탈도를 계산, 1 에 가장 가까운 \(k\) 선택. 음이항의 saddlepoint 근사에서는 이 조건이 ML 과 유사한 추정을 준다.

3.4 Pearson \(X^2 = n - p\)

\[ \hat k \text{ s.t. } \sum_i \frac{(y_i - \hat\mu_i)^2}{V(\hat\mu_i; k)} = n - p \]

이쪽이 계산적으로 더 안정적이다. Pearson 통계량은 이탈도보다 \(k\) 에 대한 감도가 크고 단조 — 격자 탐색에서 \(k\) 를 유일하게 결정한다.

3.5 직관: 네 추정량이 왜 다른가

같은 \(V\) 를 다른 목적함수 로 접근한다.

방법	무엇을 최적화	언제 쓰나
ML	정확 로그우도	원리적 최적, 표본이 크면 SE 도 최적
Moment	처음 두 모멘트	빠른 초기값, 단일 표본
Deviance = 1	분포의 “올바름” 기준	과산포 정도의 보정용
Pearson = \(n-p\)	분산함수의 경험적 적합	수렴 안정, \(V\) 가 맞는지 감각적 판단

일반적으로 ML 이 효율적이지만, 표본 크기가 작거나 분포 가정이 불확실할 때는 Pearson 매칭이 더 robust 하다.

4 음이항 적합 알고리즘 — Profile + IRLS

4.1 이중 반복 구조

외부 루프: k 갱신
  ├── 내부 루프: k 고정, 로그 링크로 GLM 적합 (IRLS)
  │       → hat_beta, hat_mu_i
  └── k 업데이트: 스코어 방정식 or moment 매칭
수렴: |k_new - k_old| < tolerance

4.2 왜 이 구조가 필요한가

\(V(\mu; k)\) 가 \(k\) 에 의존하므로 가중치 \(W = \text{diag}(1/V(\mu_i; k))\) 도 \(k\) 와 얽힌다. \(\boldsymbol\beta\) 와 \(k\) 를 동시에 풀면 연립 비선형 방정식이 되지만 — 이를 변수 분리(variable splitting) 로 접근하면 각 하위 문제가 단순해진다.

\(k\) 고정 시 \(\boldsymbol\beta\): 표준 IRLS (닫힌형 반복)
\(\boldsymbol\beta\) 고정 시 \(k\): 일변수 비선형 방정식 (Brent or Newton)

이 전략은 coordinate descent 의 예시로, mixed model 의 \((\boldsymbol\beta, \boldsymbol\sigma)\) 분리 적합과 같은 철학이다.

4.3 Profile likelihood 로 신뢰구간 구성

\(k\) 의 프로파일 로그우도:

\[ \ell_p(k) = \max_{\boldsymbol\beta} \ell(\boldsymbol\beta, k) \]

95% 신뢰구간:

\[ \left\{ k : 2(\ell_p(\hat k) - \ell_p(k)) < \chi^2_{1, 0.05} \right\} \]

Wald 구간보다 profile 이 나은 이유: \(k\) 의 로그우도는 보통 비대칭(오른쪽 꼬리가 길다). \(\hat k\) 의 표준오차로 만든 대칭 구간은 상한을 과소, 하한을 음수 까지 보낼 위험이 있다 (\(k > 0\) 이어야 함에도). Profile 은 이 경계를 자연스럽게 존중한다.

5 Manton et al. (1981) — \(\alpha\) 와 \(k\) 모두 공변량 의존

기본 음이항 GLM: \(\alpha\) 만 공변량에 따라 변하고 \(k\) 는 전역 상수. 이는 “과산포 강도가 모든 관측치에서 동일” 이라는 제약.

Manton et al. 은 둘 다 공변량에 의존:

\[ \log \alpha_i = \mathbf x_i^T \boldsymbol\beta, \qquad \log k_i = \mathbf u_i^T \boldsymbol\gamma \]

언제 필요한가: 집단별로 과산포 정도가 다른 것이 실질적 관심사일 때.

보험: 운전자 유형별 claim 평균과 과산포 모두 다를 수 있음
질병 지도: 지역별 발생 rate 와 개인 간 변동 모두 추정 대상
생태: 서식지별 개체수 평균과 집락(clustering) 강도 분리

이 모형은 McCullagh-Nelder 의 표준 GLM 프레임워크를 벗어난다 — Ch.10 의 joint mean-dispersion 과 구조적으로 유사하나, \(k\) 가 dispersion parameter 가 아니라는 점에서 다르다.

5.1 왜 Ch.10 의 이중 GLM 으로 부족한가

Ch.10 의 이중 GLM 은 평균·산포 두 선형 예측자 를 병렬로 두고 두 루프를 교대로 적합한다. 이것이 가능한 이유는 Fisher 정보 행렬이 \((\boldsymbol\beta, \boldsymbol\gamma)\) 에 대해 블록 대각 이기 때문 — \(\phi\) 변화가 \(\hat\beta\) 를 “평행 이동” 시키지 않는다 (§09-1 의 기하학적 직관).

그러나 음이항의 \(k\) 는 산포 모수가 아니라 분산함수 \(V(\mu; k) = \mu + \mu^2/k\) 의 형태 모수 다. \(\log k_i = \mathbf u_i^T \boldsymbol\gamma\) 로 놓으면 \(V\) 자체가 관측치마다 곡선 모양을 바꾸고, \(\hat\beta\) 의 정규 방정식에서 \(\gamma\) 를 깨끗이 분리할 수 없다. 블록 대각 구조가 깨지므로 이중 IRLS 대신 \((\beta, \gamma)\) 결합 Newton-Raphson 또는 Laplace/AGQ 근사가 필요하다.

5.2 현대적 구현

NBMM (Negative Binomial Mixed Model): 개체별 무작위 효과 \(b_i\) 를 도입해 \(Y_i \mid b_i \sim \text{NB}(\mu_i, k)\) 로 조건부 음이항. R glmmTMB::glmmTMB(..., family = nbinom2) 가 표준이며, Python 은 statsmodels 의 GLMM 확장이 제한적으로 지원한다
GAMLSS (Generalized Additive Models for Location, Scale and Shape): 분포의 모든 모수 (\(\mu\), \(\sigma\), \(\nu\), \(\tau\)) 를 각각 독립된 공변량 선형 예측자로 설명하는 일반 프레임워크. R gamlss 패키지가 NB 의 \(k\) 공변량 의존을 표준 지원하며, Box-Cox-t·Box-Cox-Cole-Green 등 다중 모수 분포까지 확장 가능

실무 지침: \(k\) 가 공변량에 의존한다는 구체적 근거 (예: 집단별 \(\hat k\) 추정치의 체계적 차이) 가 있을 때만 이 확장을 고려한다. 근거 없이 도입하면 식별성 붕괴로 \(\hat\gamma\) 의 표준오차가 폭발하고 수렴이 실패하는 경우가 흔하다.

6 반올림 오차 분산함수 — \(V = \tau^2 + \sigma^2 \mu^2\)

6.1 동기 — 기록 관행이 만드는 이분산성

감마 오차를 가진 양의 연속 측정값 (예: 반응 시간, 수명, 소득)을 고정 소수점 으로 기록하면 분산 구조가 어떻게 바뀌는가.

순수 감마: \(\mathrm{var}(Z) = \sigma^2 \mu^2\) (\(\sigma = 1/\sqrt{\text{shape}}\), 비례 오차) 기록 단위 \(\epsilon = 10^{-d}\): \(Y = \epsilon \lfloor Z/\epsilon \rceil\) (가장 가까운 정수 배로 반올림)

6.2 Exercise 11.1 의 핵심 결과

Kolassa & McCullagh (1987) 증명: \(\epsilon\) 이 작을 때 \(R = (Z - Y)/\epsilon\) 은 \([-1/2, 1/2]\) 의 균등분포 에 점근적으로 독립이다. 정확하게는

\[ \kappa_r(Y) \simeq \begin{cases} \kappa_r(Z) & r \text{ 홀수} \\ \kappa_r(Z) + \epsilon^r \kappa_r(R) & r \text{ 짝수} \end{cases} \]

\(R\) 의 누율: \(\kappa_2(R) = 1/12\), \(\kappa_4(R) = -1/120\), 홀수 누율은 0.

6.3 직관: 왜 짝수 누율만 증가하는가

\(R\) 은 대칭 분포 (\([-1/2, 1/2]\) 의 uniform). 따라서 홀수 모멘트와 누율이 0. 반올림이 가하는 “무작위 노이즈” 는 평균과 비대칭을 건드리지 않고 분산·첨도만 부풀린다.

\(Z\) 와 \(Y\) 의 첫째 누율(평균)은 같다 — 반올림이 편향을 만들지 않는다 (대칭 반올림이므로). 둘째 누율(분산)은

\[ \mathrm{var}(Y) \simeq \mathrm{var}(Z) + \epsilon^2 \kappa_2(R) = \mathrm{var}(Z) + \epsilon^2/12 \tag{11.5} \]

상수 \(\epsilon^2/12\) 가 분산에 더해진다. 이것이 \(V(\mu) = \sigma^2 \mu^2 + \epsilon^2/12\) 로 나타나는 이유다. 통상 \(\tau^2 \equiv \epsilon^2/12\) 로 쓴다.

6.4 Exercise 11.2 — 감마 분포 + \(d\) 소수점 반올림

\(Z \sim \text{Gamma}(\text{shape} = \nu, \text{mean} = \mu_Z)\), \(\mathrm{var}(Z) = \mu_Z^2/\nu\). 반올림 후

\[ \mathrm{var}(Y) \simeq \frac{\mu_Z^2}{\nu} + \frac{\epsilon^2}{12} \tag{11.6} \]

첫 항은 비례 오차, 둘째 항은 상수 반올림 오차. 유효 숫자 반올림 (scientific notation 의 유효 자릿수 \(d\)) 이면 \(\epsilon \propto \mu\) 이므로 둘째 항도 \(\mu^2\) 에 비례 — 분산 구조가 다시 \(\propto \mu^2\) 로 돌아온다.

6.5 분산 구조가 주는 실무적 함정

순수 감마 GLM 은 \(V(\mu) \propto \mu^2\) 를 가정하여 가중치 \(w_i \propto 1/\mu_i^2\) 을 쓴다. 작은 \(\mu\) 의 관측치가 큰 가중치를 받는다.

문제: 작은 \(\mu\) 에서는 반올림 오차 \(\epsilon^2/12\) 이 비례 오차를 지배 하므로 실제 분산이 모형보다 크다. 가중치가 과대평가 → 계수 추정이 작은 \(\mu\) 데이터에 과적합.

수정된 가중치:

\[ w_i = \frac{1}{\tau^2 + \sigma^2 \mu_i^2} \]

작은 \(\mu_i\) 에서 \(\tau^2\) 이 지배, 큰 \(\mu_i\) 에서 \(\sigma^2 \mu_i^2\) 이 지배. 상대 가중치가 부드럽게 전환된다.

6.6 \(\sigma^2/\tau^2\) 의 추정

음이항의 \(k\) 와 같은 방식.

\(\tau^2, \sigma^2\) 초기값 설정 (예: \(\tau^2 = \epsilon^2/12\) 알려진 값, \(\sigma^2 = 1\))
로그 링크 + \(V(\mu; \tau^2, \sigma^2)\) 로 quasi-likelihood GLM 적합 → \(\hat\mu_i\)
Pearson 통계량이 \(n - p\) 가 되도록 \(\sigma^2/\tau^2\) 조정
수렴까지 반복

Quasi-likelihood 체계(McCullagh §9.2) 에서 분포 가정 없이 \(V(\mu; k)\) 만 올바르면 \(\hat\beta\) 의 consistency 가 보장된다. \(\tau^2, \sigma^2\) 자체가 실무적 관심사가 아닐 때는 Pearson 매칭으로 충분.

7 응용 분야

분야	활용	구체적 예시
공정 품질	결함 수 count	반도체 wafer 의 defect count, 과산포 음이항
보건 역학	질병 발생 rate	지역별 희귀 암 count, Poisson 대신 NB 로 overdispersion 흡수
마케팅	클릭·구매 수	광고 노출당 클릭 분포, NB regression 의 \(k\) 가 고객 이질성
운영 리서치	콜센터 호출 수	시간대·팀별 호출 count, 과산포가 작업 로드 변동을 반영
반응 시간 측정	ms 단위 기록	심리학 실험의 RT 데이터, 반올림 분산함수
금융 데이터	가격·거래량	정수 tick size 반올림, \(\tau^2\) 성분이 microstructure 노이즈
천문학	사진 측광	디지털 CCD count 의 Poisson 노이즈 + 읽기 잡음 (상수 \(\tau^2\))
생태 개체수	Count 데이터	종 분포의 과산포, NB regression 기본 도구

8 예시 — 자동차 보험 claim 수 (음이항 적용)

운전자 \(i\) 의 claim 수 \(Y_i\), 연령·경험·지역 공변량 \(\mathbf x_i\). 포아송 모형을 적합했더니 Pearson \(X^2/(n-p) = 2.4\) → 과산포.

음이항 모형:

\[ \log \mu_i = \mathbf x_i^T \boldsymbol\beta, \qquad \mathrm{var}(Y_i) = \mu_i + \mu_i^2/k \]

\(\hat k = 1.8\) 추정. 해석: \(\mu = 0.1\) 에서 포아송 분산은 0.1, 음이항 분산은 \(0.1 + 0.01/1.8 = 0.106\) (거의 포아송과 같음). \(\mu = 2\) 에서 포아송 2, 음이항 \(2 + 4/1.8 = 4.22\) (두 배 이상). 과산포가 고위험 운전자에서 집중적으로 나타난다 는 실무적 통찰을 수식이 그대로 보여준다.

9 코드 예시 — 음이항 + 반올림 분산

9.1 Step 1: Python — NB 의 \(k\) 직접 추정

import numpy as np
from scipy.special import digamma, gammaln
from scipy.optimize import brentq

rng = np.random.default_rng(7)
n = 300
x = rng.normal(size=n)
mu = np.exp(0.3 + 0.5 * x)
k_true = 2.0
# Negative binomial 샘플링: Poisson-Gamma 혼합
lam = rng.gamma(k_true, mu / k_true)
y = rng.poisson(lam)

def nb_glm_inner(x, y, k, n_iter=50, tol=1e-8):
    """k 고정하고 NB GLM (로그 링크) 적합 — IRLS."""
    X = np.column_stack([np.ones(n), x])
    beta = np.array([np.log(y.mean() + 1e-6), 0.0])
    for _ in range(n_iter):
        eta = X @ beta
        mu_hat = np.exp(eta)
        V = mu_hat + mu_hat ** 2 / k
        w = mu_hat ** 2 / V                        # 로그 링크: dmu/deta = mu
        z = eta + (y - mu_hat) / mu_hat
        WX = X * w[:, None]
        beta_new = np.linalg.solve(X.T @ WX, X.T @ (w * z))
        if np.max(np.abs(beta_new - beta)) < tol:
            beta = beta_new
            break
        beta = beta_new
    return beta, mu_hat

def nb_loglik(k, x, y):
    beta, mu_hat = nb_glm_inner(x, y, k)
    ll = np.sum(
        gammaln(y + k) - gammaln(k) - gammaln(y + 1)
        + k * np.log(k / (mu_hat + k))
        + y * np.log(mu_hat / (mu_hat + k))
    )
    return ll

# Profile likelihood 로 k 탐색
k_grid = np.linspace(0.3, 8.0, 40)
lls = np.array([nb_loglik(k, x, y) for k in k_grid])
k_hat_profile = k_grid[np.argmax(lls)]

# Brent 법으로 스코어 = 0 찾기 (정밀화)
def score_k(k):
    beta, mu_hat = nb_glm_inner(x, y, k)
    return np.sum(
        digamma(y + k) - digamma(k)
        + np.log(k / (mu_hat + k))
        + (y - mu_hat) / (mu_hat + k)
    )

k_hat = brentq(score_k, 0.3, 8.0)
beta_hat, _ = nb_glm_inner(x, y, k_hat)

print(f"k_hat (profile grid) = {k_hat_profile:.3f}")
print(f"k_hat (Brent)        = {k_hat:.3f}   (true = {k_true})")
print(f"beta_hat             = {beta_hat}  (true = [0.3, 0.5])")

gammaln 을 이용해 factorial 을 로그 스케일에서 안정적으로 계산한다. 내부 IRLS 는 mu_hat.var() 이 큰 값에서도 overflow 없이 동작.

9.2 Step 2: R — MASS::glm.nb 표준 적합

library(MASS)
set.seed(7)
n <- 300
x <- rnorm(n)
mu <- exp(0.3 + 0.5 * x)
k_true <- 2.0
lam <- rgamma(n, shape = k_true, scale = mu / k_true)
y <- rpois(n, lam)

fit <- glm.nb(y ~ x, link = log)
summary(fit)
cat("theta (k) =", fit$theta, "  true =", k_true, "\n")

# Profile likelihood 로 k 의 95% 신뢰구간
confint(profile(fit))               # 각 계수
# k 의 profile CI 는 fit$SE.theta 와 fit$theta 로부터 비대칭 구성
k_profile <- MASS:::profile.glm(fit, which = "theta")

glm.nb 는 nbinomial family 의 특수화로, \(k\) (R 에서는 theta) 의 ML 추정을 digamma 방정식으로 내부 처리한다.

9.3 Step 3: 반올림 분산함수 시뮬레이션

import numpy as np
import statsmodels.api as sm

rng = np.random.default_rng(42)
n = 500
x = rng.uniform(0.5, 5.0, n)
mu_true = np.exp(1.0 + 0.3 * x)
shape = 10.0                                    # sigma^2 = 1/shape = 0.1
z = rng.gamma(shape, mu_true / shape)           # 순수 감마

eps = 0.1                                        # 소수점 1자리 반올림
y = eps * np.round(z / eps)

# 순수 감마 GLM (잘못된 모형)
X = sm.add_constant(x)
fit_gamma = sm.GLM(y, X, family=sm.families.Gamma(link=sm.families.links.log())).fit()

# Quasi-likelihood with V(mu) = tau^2 + sigma^2 mu^2
tau2 = eps ** 2 / 12
sigma2 = 1.0 / shape

def irls_with_rounding(y, X, tau2, sigma2_guess, n_iter=30):
    beta = np.array([np.log(y.mean()), 0.0])
    sigma2 = sigma2_guess
    for _ in range(n_iter):
        eta = X @ beta
        mu = np.exp(eta)
        V = tau2 + sigma2 * mu ** 2
        w = mu ** 2 / V
        z = eta + (y - mu) / mu
        WX = X * w[:, None]
        beta = np.linalg.solve(X.T @ WX, X.T @ (w * z))
        # Pearson 매칭으로 sigma2 갱신
        resid = y - np.exp(X @ beta)
        pearson = np.sum(resid ** 2 / (tau2 + sigma2 * np.exp(X @ beta) ** 2))
        sigma2 *= pearson / (n - X.shape[1])
    return beta, sigma2

beta_rob, sigma2_hat = irls_with_rounding(y, X, tau2, sigma2)
print(f"순수 감마 GLM beta  = {fit_gamma.params}")
print(f"반올림 보정 beta    = {beta_rob}  (true = [1.0, 0.3])")
print(f"sigma^2 estimate   = {sigma2_hat:.4f}  (true = {1/shape})")

반올림 보정 모형은 작은 \(\mu\) 에서 \(\tau^2\) 성분이 가중치를 안정시켜 절편 추정이 덜 왜곡된다.

10 진단과 실무 가이드

진단	어떤 문제를 잡는가
Pearson \(X^2/(n-p)\)	포아송 적합 후 1 초과 → 과산포, 음이항 고려
\(\hat k\) 의 profile CI	Wald 구간이 0 을 포함하면 profile 로 확인 — \(k > 0\) 경계
\(\log V(\hat\mu_i)\) vs \(\log \hat\mu_i\)	기울기 2 → \(V \propto \mu^2\) (감마), 기울기 1 → 포아송, 그 사이 → NB
작은 \(\hat\mu\) 잔차의 이상 분산	반올림·검출 한계 등 상수 분산 성분 의심 → \(V = \tau^2 + \sigma^2 \mu^2\)
NB vs Poisson LRT	\(H_0: 1/k = 0\) 의 경계 검정 — \(\chi^2_1\) 의 \(1/2\) 혼합

실무 체크리스트:

포아송 GLM 부터 적합. Pearson \(X^2/(n-p) \gg 1\) 이면 과산포 진단
과산포 원인 고려: 진짜 혼합 (NB), 영 초과 (ZINB), 무작위 효과 (mixed model)? NB 는 가장 간단
\(\hat k\) 가 매우 작으면 (< 0.5) NB 의 이차항이 과대 — 데이터에 극단치나 혼합 구조 재점검
\(\hat k\) 가 매우 크면 (> 50) 포아송과 구분 안 됨 — NB 사용의 근거 약함
측정 반올림이 심한 continuous outcome 은 \(V = \tau^2 + \sigma^2 \mu^2\) quasi-likelihood 로 가중치 보정

11 Bibliographic Notes

Anscombe (1949): 음이항 \(k\) 의 모멘트·최대우도 추정. 단일 표본·다수 표본 처리
Cameron & Trivedi (1998): Count regression 의 현대 표준 교과서. NB1, NB2 구분 체계
Manton et al. (1981): \(\alpha\) 와 \(k\) 모두 공변량 의존시키는 일반화. 지역별 사망률 데이터
Kolassa & McCullagh (1987): Rounding error 의 점근 누율 유도 (Exercise 11.1)
Hilbe (2011): Negative Binomial Regression 전용 단행본. 수렴 문제·profile likelihood 상세
Lawless (1987): NB regression 의 ML vs moment 추정 비교
Breslow (1984): quasi-likelihood 관점에서의 과산포 처리

12 정리

분산함수 속 구조 모수는 산포 모수 \(\phi\) 와 근본적으로 다르다. \(\hat\beta\) 의 정규 방정식에서 분리되지 않아 이중 반복 알고리즘이 필요하며, 분포 가정의 강도가 링크·공분산 해석에까지 영향을 미친다.

음이항의 \(k\) 는 Poisson-Gamma 혼합의 집중도로 해석되며, \(V(\mu) = \mu + \mu^2/k\) 의 이차항을 통해 과산포를 흡수한다
Canonical link 는 \(k\) 에 의존하므로 실무에서는 로그 링크 + NB 분산함수 조합이 표준
\(k\) 의 추정은 ML (digamma 방정식), moment, deviance=1, Pearson=\(n-p\) 의 네 경로 — 표본 크기·안정성에 따라 선택
Profile likelihood 는 Wald 보다 비대칭 CI 를 정확히 반영하여 \(k\) 의 경계(\(k > 0\))를 존중
반올림 오차는 상수 분산 성분 \(\tau^2 = \epsilon^2/12\) 를 추가하여 \(V(\mu) = \tau^2 + \sigma^2 \mu^2\) 꼴로 나타나며, 작은 \(\mu\) 에서 가중치를 안정화

다음 절(§11.3) 은 연결함수 속 비선형 모수 — Box-Cox 멱 링크와 Pregibon 선형화 — 를 다룬다.

13 관련 주제

선행 지식

Models with Additional Non-Linear Parameters — 개관 (McCullagh Ch.11) — 세 출구 전체 지도
Likelihood Functions for Log-linear Models — 포아송 로그우도·NB 등장 (McCullagh §6.2) — NB 의 GLM 맥락
Quasi-likelihood Functions 개관 (McCullagh Ch.9) — 분산함수만 쓰는 추론 체계
이항 자료의 과산포 — Over-dispersion (McCullagh §4.5) — 과산포 처리 3 메커니즘

후속 주제

Parameters in the Link Function — Box-Cox \(\lambda\)·Pregibon 선형화 (McCullagh §11.3)
Non-linear Parameters in the Covariates — Box-Tidwell (McCullagh §11.4)
Zero-Inflated Count Models — NB 를 넘어선 과산포 처리 (현대 확장)
Negative Binomial Mixed Model — 다수준 count 데이터 (Manton 확장의 현대 버전)

관련 개념

Extended Quasi-likelihood — Q⁺ 정의 (McCullagh §9.6) — 산포 추정의 일반 체계
Adjustments of the Estimating Equations — 첨도 보정 (McCullagh §10.5) — \(\rho_4\) 가 NB 에서 어떻게 계산되는가
Gamma GLM — 분산 함수·이탈도 (McCullagh §8.3) — \(V \propto \mu^2\) 의 순수 형태