1 도입 — 왜 우도 없이 추론하려 하는가
Ch.2–Ch.8 은 모두 지수족 분포를 출발점으로 삼았다. 이항·포아송·정규·감마·역가우시안은 각기 고유한 로그우도 \(\ell(\beta) = \sum \log f(y_i; \mu_i, \phi)\) 를 갖고, MLE·이탈도·Wald·LRT 가 이 위에서 자연스럽게 정의된다.
그러나 실무에서는 다음과 같은 상황이 빈번하다.
- 자료의 분포에 대한 이론적 근거가 없다: 과거 유사 실험이 없거나, 생성 기제가 알려지지 않은 자료.
- 일부 특징만 안다: 반응이 양수이다, 퍼센트이다, 평균이 커지면 분산도 커진다 — 이런 질적 정보는 있지만 완전한 PDF 는 특정하기 어렵다.
- 이항·포아송 모형을 썼더니 과산포(overdispersion)가 관찰된다: 분산이 \(\mu(1-\mu)/m\) 보다 크거나 \(\mu\) 보다 크다.
이때 지수족 우도를 억지로 가정하면, 점추정은 여전히 합리적일 수 있지만 표준오차와 검정 통계량이 체계적으로 왜곡된다. 이 장은 완전한 우도를 포기하고, 대신 평균-분산 관계만으로 추론 가능한가를 묻는다. 답은 그렇다이며, 이 틀을 준우도(quasi-likelihood)라 부른다(McCullagh & Nelder, 1989, §9.1).
직관. 로그우도의 미분, 즉 스코어 함수 \(U = \partial \ell / \partial \mu\) 가 추정·검정의 거의 모든 1차 점근이론을 떠받친다. 놀랍게도 스코어의 핵심 세 성질은 분포 전체가 아니라 처음 두 모멘트에만 의존한다. 그렇다면 분포를 모른 채 이 세 성질만 갖는 함수를 인위적으로 만들어 스코어처럼 쓰면 된다. 이것이 준우도의 아이디어다.
2 출발점 — 로그우도 미분의 세 성질
관측 \(Y\) 가 평균 \(\mu\) 와 분산 \(\operatorname{var}(Y) = \sigma^2 V(\mu)\) 를 갖는다고 하자. 여기서 \(V(\cdot)\) 는 분산함수(variance function)라 불리는 알려진 함수이고, \(\sigma^2\) 는 (보통 미지의) 산포 모수이다. 분포 \(f\) 는 가정하지 않는다.
함수
\[ U = u(\mu; Y) = \frac{Y - \mu}{\sigma^2 V(\mu)} \]
를 정의한다. 구조적 의미: 분자 \((Y - \mu)\) 는 관측과 평균의 생 잔차, 분모 \(\sigma^2 V(\mu)\) 는 그 잔차의 분산이다. 따라서 \(U\) 는 “잔차를 자신의 분산 역수로 가중한 표준화 신호” — 분산이 큰 구간(정보가 약함) 에서는 가중을 줄이고, 분산이 작은 구간(정보가 강함) 에서는 가중을 키우는 최적 가중 의 일반 형태다. 계산하면 다음 세 성질을 갖는다.
\[ E(U) = 0, \qquad \operatorname{var}(U) = \frac{1}{\sigma^2 V(\mu)}, \qquad - E\!\left(\frac{\partial U}{\partial \mu}\right) = \frac{1}{\sigma^2 V(\mu)}. \]
비교. 만약 로그우도가 존재한다면, 스코어 \(\partial \ell/\partial \mu\) 도 정확히 같은 세 성질을 갖는다(Bartlett 항등식). 즉, 1차 점근이론의 근간이 되는 세 성질을 \(U\) 가 완비한다. 분포 전체를 가정하지 않고도, 평균-분산 관계만으로 스코어에 준하는 함수를 확보한 셈이다.
평균 \(E(Y) = \mu(\beta)\) 와 분산 \(\operatorname{var}(Y) = \sigma^2 V(\mu)\) 만 가정할 때,
\[ U(\beta) = \frac{1}{\sigma^2}\,\mathbf{D}^T \mathbf{V}^{-1}(\mathbf{Y} - \boldsymbol{\mu}) \]
을 준-스코어 함수라 한다. 여기서 \(\mathbf{D}_{ir} = \partial \mu_i / \partial \beta_r\) 이며 \(\mathbf{V} = \operatorname{diag}\{V_1(\mu_1), \ldots, V_n(\mu_n)\}\) (독립 관측의 경우)이다.
왜 이 형태가 나오나. GLM 스코어 \(U_r = \sum (y_i - \mu_i)\,(\partial \mu_i/\partial \beta_r)/\{\sigma^2 V(\mu_i)\}\) 의 정확한 베꼈다. 지수족 GLM 의 스코어 방정식은 오직 \(\mu_i\), \(V(\mu_i)\), \(\partial \mu_i/\partial \beta_r\) 로만 쓰여진다는 사실이 준우도의 존재 이유다 — 분포를 모두 가정하지 않아도 이 세 재료만 있으면 같은 추정방정식을 세울 수 있다.
3 독립 관측의 준우도
3.1 준우도의 정의
독립 관측 하에서 공분산 행렬은 \(\sigma^2 V(\mu) = \sigma^2 \operatorname{diag}\{V_i(\mu_i)\}\) 가 된다. 각 성분 \(Y_i\) 에 대해 준우도는
\[ Q(\mu_i; y_i) = \int_{y_i}^{\mu_i} \frac{y_i - t}{\sigma^2 V_i(t)}\,dt \tag{9.3} \]
로 정의되며, 전체 표본에 대해서는 단순 합
\[ Q(\boldsymbol{\mu}; \mathbf{y}) = \sum_{i=1}^{n} Q_i(\mu_i; y_i) \]
이 된다.
직관. 적분한이다. \(\partial Q / \partial \mu = -(y-\mu)/\{\sigma^2 V(\mu)\} = -U\) 이므로, \(Q\) 를 \(\mu\) 에 대해 미분하면 준-스코어가 그대로 나온다(부호는 로그우도 미분 관행과 맞춘다). 즉, 스코어를 적분하여 우도스러운 스칼라 함수를 복원한 결과가 \(Q\) 이다.
3.2 Table 9.1 — 분산함수와 대응되는 준우도
| 분산함수 \(V(\mu)\) | 준우도 \(Q(\mu; y)\) | 정준모수 \(\theta\) | 대응 분포 |
|---|---|---|---|
| \(1\) | \(-(y-\mu)^2/2\) | \(\mu\) | Normal |
| \(\mu\) | \(y\log\mu - \mu\) | \(\log \mu\) | Poisson |
| \(\mu^2\) | \(-y/\mu - \log\mu\) | \(-1/\mu\) | Gamma |
| \(\mu^3\) | \(-y/(2\mu^2) + 1/\mu\) | \(-1/(2\mu^2)\) | Inverse Gaussian |
| \(\mu^\zeta\ (\zeta\neq 0,1,2)\) | \(\mu^{-\zeta}\!\left(\frac{\mu y}{1-\zeta} - \frac{\mu^2}{2-\zeta}\right)\) | \(\frac{1}{(1-\zeta)\mu^{\zeta-1}}\) | Tweedie 계열 |
| \(\mu(1-\mu)\) | \(y\log\!\left(\frac{\mu}{1-\mu}\right) + \log(1-\mu)\) | \(\log\!\left(\frac{\mu}{1-\mu}\right)\) | Binomial/\(m\) |
| \(\mu^2(1-\mu)^2\) | \((2y-1)\log\!\left(\frac{\mu}{1-\mu}\right) - \frac{y}{\mu} - \frac{1-y}{1-\mu}\) | — | — |
| \(\mu + \mu^2/k\) | \(y\log\!\left(\frac{\mu}{k+\mu}\right) + k\log\!\left(\frac{k}{k+\mu}\right)\) | \(\log\!\left(\frac{\mu}{k+\mu}\right)\) | Negative binomial |
핵심 관찰. 앞 네 행은 모두 실제 분포(지수족 1-모수)의 로그우도와 같다. 다섯 번째 행 \(V(\mu) = \mu^\zeta\) 는 Tweedie 분포족에 대응한다. 여섯 번째 \(V(\mu) = \mu(1-\mu)\) 는 이항 로그우도와 동일하지만, 이 경우에는 \(y\) 가 \(0/1\) 이 아니어도 — 예를 들어 \(y \in [0,1]\) 의 연속 비율이어도 — 그대로 준우도로 쓸 수 있다. 일곱 번째 \(\mu^2(1-\mu)^2\) 는 알려진 분포가 없지만 퍼센트 자료의 분산이 평균이 극단값에 가까울수록 더 작아지는 경우를 표현한다.
자료 유형별 매핑. 각 분산함수는 서로 다른 자료 형태의 “분산이 평균과 함께 어떻게 커지는가” 를 규정한다. \(V = \mu\) 는 카운트 (한 단위 시간·면적당 사건 수) — 평균이 크면 분산도 선형으로 커진다. \(V = \mu(1-\mu)\) 는 비율/이항 확률 — 양 끝에서 분산이 0 에 수렴. \(V = \mu^2\) 는 우측 편포 양수 자료 (금액·반응 시간) — 변동계수가 평균에 무관하게 일정. \(V = \mu^3\) 은 대기 시간·수명 의 역가우스 구조. \(V = \mu + \mu^2/k\) 는 과산포 카운트 의 음이항. 한 줄로 말하면: 자료가 어떤 스케일에서 “평평한 오차” 를 보이는가 를 선택하면 대응하는 \(V\) 가 결정된다.
so what. 분포를 가정하지 않은 채, 분산함수 하나만 지정하면 GLM 추정 엔진이 그대로 돌아간다. 이것이 준우도의 실용적 승리다.
3.3 준-이탈도(Quasi-deviance)
로그우도 비와 유사한 양으로 준-이탈도
\[ D(y; \mu) = -2\sigma^2\, Q(\mu; y) = 2\int_{\mu}^{y} \frac{y-t}{V(t)}\,dt \tag{9.4} \]
가 정의된다. 성질은 다음과 같다.
- \(y = \mu\) 에서 \(0\), 그 외에서 양수이다(분산함수가 양이므로 피적분 함수의 부호가 \((y-t)\) 의 부호와 같다).
- \(\sigma^2\) 에 의존하지 않는 순수 \((y, \mu)\) 함수이다.
- 전체 이탈도 \(D(\mathbf{y}; \boldsymbol{\mu}) = \sum D(y_i; \mu_i)\) 가 모형 비교에 사용된다.
직관. 준-이탈도는 “각 관측이 제 평균에서 얼마나 벗어나 있는가”를 분산함수 \(V\) 로 가중하여 측정한다. \(V\) 가 작은 구간에서는 같은 \(|y-\mu|\) 이라도 이탈도가 크게 잡힌다 — 정밀 관측 구간에서의 벗어남은 더 심각하다.
4 준우도 추정
4.1 추정방정식과 IRLS
회귀계수 \(\beta\) 에 대한 준우도 추정방정식은 \(U(\widehat{\beta}) = 0\) 이다. 즉
\[ \mathbf{D}^T \mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) = \mathbf{0} \tag{9.5} \]
를 푼다. 준-정보행렬은
\[ \mathbf{i}_\beta = \mathbf{D}^T \mathbf{V}^{-1} \mathbf{D}/\sigma^2 \tag{9.6} \]
이며 Fisher 정보와 같은 역할을 한다. 점근 공분산은
\[ \operatorname{cov}(\widehat{\beta}) \simeq \mathbf{i}_\beta^{-1} = \sigma^2 (\mathbf{D}^T \mathbf{V}^{-1} \mathbf{D})^{-1}. \]
반복 절차는 Fisher 스코어링과 동일하다.
\[ \widehat{\beta}_{k+1} = \widehat{\beta}_k + (\widehat{\mathbf{D}}_k^T \widehat{\mathbf{V}}_k^{-1} \widehat{\mathbf{D}}_k)^{-1}\widehat{\mathbf{D}}_k^T \widehat{\mathbf{V}}_k^{-1}(\mathbf{y} - \widehat{\boldsymbol{\mu}}_k). \]
중요 관찰. 이 반복은 \(\sigma^2\) 에 무관하다. 스칼라 \(\sigma^2\) 가 \(\mathbf{V}^{-1}\) 과 \((\mathbf{D}^T \mathbf{V}^{-1} \mathbf{D})^{-1}\) 양쪽에서 상쇄되기 때문이다. 따라서 점추정은 \(\sigma^2\) 를 몰라도 진행된다. IRLS 알고리즘을 그대로 활용할 수 있다 — 실제로 R::glm(family=quasi(...)) 과 statsmodels.GLM(family=...) 는 이 원리를 쓴다.
4.2 1단계 추정량의 선형 표현
참값 \(\beta\) 에서 출발한 1단계 갱신은
\[ \widehat{\beta}_1 = \beta + (\mathbf{D}^T \mathbf{V}^{-1} \mathbf{D})^{-1}\mathbf{D}^T \mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) \tag{9.7} \]
가 된다. 이 식으로부터
- 거의 비편향(\(E(\widehat{\beta}_1) \approx \beta\)),
- 점근 정규성(\(\widehat{\beta}_1\) 이 \((\mathbf{y}-\boldsymbol{\mu})\) 의 선형 결합이며 다수의 독립 성분의 합이므로 CLT 적용),
- 점근 분산 \(\sigma^2 (\mathbf{D}^T \mathbf{V}^{-1} \mathbf{D})^{-1}\)
이 모두 2차 모멘트 가정만으로 따라 나온다. 분포 전체를 쓰지 않고도 MLE 와 같은 점근 결과를 얻는 것이 준우도의 이론적 핵심 성과이다.
4.3 산포 모수의 추정
\(\sigma^2\) 은 우도 기반이 아닌 모멘트 추정으로 구한다.
\[ \widetilde{\sigma}^2 = \frac{1}{n-p}\sum_i \frac{(Y_i - \widehat{\mu}_i)^2}{V_i(\widehat{\mu}_i)} = \frac{X^2}{n-p}, \]
여기서 \(X^2\) 는 일반화 Pearson 통계량이다. 이 추정량은 분포 가정에 덜 민감하여 강건하다(McCullagh & Nelder, 1989, §9.2.3).
왜 이탈도 기반이 아닌가. 일반 MLE 에서는 \(\sigma^2\) 도 \(\partial Q/\partial \sigma^2 = 0\) 으로 구할 수 있지만, 준우도에서 \(Q\) 는 \(\sigma^2\) 를 정규화 상수처럼 포함만 할 뿐 \(\sigma^2\) 에 대한 별도의 정보를 주지 않는다. 따라서 잔차의 모멘트로부터 직접 추정한다.
5 종속 관측으로의 확장
5.1 공분산 구조의 일반화
공분산이 \(\operatorname{cov}(\mathbf{Y}) = \sigma^2 \mathbf{V}(\boldsymbol{\mu})\) 이되 \(\mathbf{V}\) 가 대각이 아닌 경우(예: 종단 자료의 피험자 내 상관, 공간 상관)에도 준-스코어 식 (9.5)는 그대로 사용된다.
\[ U(\beta) = \mathbf{D}^T \mathbf{V}^{-1}(\mathbf{Y} - \boldsymbol{\mu})/\sigma^2 \]
세 성질 \(E(U)=0\), \(\operatorname{cov}(U) = \mathbf{D}^T \mathbf{V}^{-1} \mathbf{D}/\sigma^2\), \(-E(\partial U/\partial \beta) = \mathbf{D}^T \mathbf{V}^{-1} \mathbf{D}/\sigma^2\) 는 모두 유지된다. 점근 분포도 같다.
실무 연결. 이 틀이 바로 Liang & Zeger (1986)의 일반화추정방정식(Generalized Estimating Equations, GEE) 이다. 종단 자료에서 피험자 내 공분산 구조를 “작업 상관”(working correlation)으로 지정하고 준-스코어로 풀면, 상관구조가 틀려도 \(\widehat{\beta}\) 는 여전히 일치 추정량이 된다(샌드위치 표준오차로 보정).
5.2 경로 독립성(Path Independence) 문제
독립 관측과 달리, \(\mathbf{V}\) 가 비대각일 때는 \(U(\beta)\) 가 어떤 스칼라 함수의 그래디언트인가가 자명하지 않다. 구체적으로 \(r \neq s\) 에 대해
\[ \frac{\partial U_r(\beta)}{\partial \beta_s} \neq \frac{\partial U_s(\beta)}{\partial \beta_r} \]
일 수 있다. 이 경우 선적분
\[ Q(\boldsymbol{\mu}; \mathbf{y}, t(s)) = \sigma^{-2}\int_{t(s)=\mathbf{y}}^{t(s)=\boldsymbol{\mu}} (\mathbf{y}-t)^T \{\mathbf{V}(t)\}^{-1} dt(s) \]
는 적분 경로에 의존한다. 즉 스칼라 함수로서의 준우도가 존재하지 않는다. 추정방정식은 쓸 수 있지만, “준우도 함수”라는 스칼라 객체는 없다.
경로 독립 조건(식 9.9). \(\mathbf{W} = \mathbf{V}^{-1}\) 의 편미분이 세 지표에 대해 대칭,
\[ \partial W_{ij}/\partial \mu_k = \partial W_{ik}/\partial \mu_j = \partial W_{jk}/\partial \mu_i, \]
이면 경로 독립이며, 이는 볼록 함수 \(b^*(\mu)\) 가 존재하여 \(\mathbf{V}^{-1} = \partial^2 b^*/\partial \mu^2\) 로 쓸 수 있다는 조건과 동치이다. 이때 정준모수 \(\theta(\mu) = \nabla b^*(\mu)\) 와 누율함수 \(b(\theta)\) 가 정의되어 지수족 구조가 자연스럽게 복원된다.
so what. 경로 독립이 성립하지 않으면 추정방정식은 유효하지만 이탈도·LRT 유사 통계량이 정의되지 않는다. 실무에서는 대부분 독립 관측이거나 단순한 클러스터 상관이므로 이 문제가 드러나지 않지만, 일반 GEE 문헌에서 “로그우도 기반 검정” 대신 Wald 또는 스코어 검정이 선호되는 이론적 배경이 바로 이것이다.
6 예시 — 보리 잎 반점(Leaf-blotch) 데이터
6.1 배경
Wedderburn (1974) 이 제시한 고전 예제이다. 10 품종 보리를 9 사이트에서 재배하고, 잎 반점병(Rhynchosporium secalis) 감염 면적의 비율 \(y \in [0,1]\) 을 측정했다. \(y\) 는 \(0/1\) 이 아닌 연속 비율이다.
6.2 모형 1: 분산함수 \(V(\mu) = \mu(1-\mu)\)
이항 로짓 GLM 을 그대로 쓰되, 반응을 비율로 취급한다.
- 선형예측자: \(\eta_{ij} = \alpha_i + \gamma_j\) (사이트 \(i\), 품종 \(j\)).
- 링크: \(\eta = \log\{\pi/(1-\pi)\}\).
- 분산함수: \(V(\pi) = \pi(1-\pi)\).
이 경우 준우도는 Table 9.1 의 \(y \log\{\mu/(1-\mu)\} + \log(1-\mu)\) 가 되어, 이항 로그우도와 수식 형태가 동일하다. 분포 가정을 안 해도 추정은 이항 GLM 과 같다.
진단. Pearson 잔차 \(r_P = (y-\widehat{\pi})/\sqrt{\widehat{\pi}(1-\widehat{\pi})/m}\) 를 선형예측자에 대해 그리면, 중앙은 0 부근이지만 양 극단(작은 \(\widehat{\pi}\) 또는 큰 \(\widehat{\pi}\))에서 잔차가 계통적으로 크다. 분산함수 \(\mu(1-\mu)\) 가 극단에서 너무 빨리 0 으로 떨어져, 실제 자료의 산포를 과소평가한다는 신호이다.
6.3 모형 2: 분산함수 \(V(\mu) = \mu^2(1-\mu)^2\)
Wedderburn 은 극단에서의 산포 축소를 더 심하게 반영하는
\[ V(\mu) = \mu^2(1-\mu)^2 \]
를 제안한다. 이 분산함수에 대응되는 알려진 분포는 없다 — 준우도의 진정한 위력을 보여주는 대목이다. Table 9.1 일곱 번째 행의 \(Q\) 를 쓰면 모든 GLM 진단이 그대로 적용된다.
- \(\widetilde{\sigma}^2 = 71.2/72 \approx 0.99\) (잔차 자유도 72 기준).
- 진단 그림에서 계통적 잔차 패턴이 크게 감소한다.
- 분산함수가 평균에 대해 더 빠르게 감소하므로, 극단 구간의 정밀도가 상대적으로 높아진다: 낮은 감염률 구간의 품종 대비는 더 정확히, 높은 감염률 구간은 덜 정확히 추정된다.
핵심 교훈. 분포를 바꾸지 않고도 분산함수를 교정하는 것만으로 모형이 크게 개선된다. “이항 유사 자료이지만 이항이 아닌” 영역에서 준우도가 빛을 낸다.
7 코드 예시
7.1 순수 구현 — 준-스코어 반복(IRLS 형태)
분산함수 \(V(\mu) = \mu^2(1-\mu)^2\), 로짓 링크의 준우도 추정을 직접 구현한다.
import numpy as np
def quasi_fit(X, y, V_fn, link_inv, link_deriv, tol=1e-8, max_iter=50):
"""
X: (n, p) 설계행렬
y: (n,) 반응 (비율 in [0,1])
V_fn(mu): 분산함수
link_inv(eta): g^{-1}(eta)
link_deriv(mu): dg/dmu (로짓: 1/(mu*(1-mu)))
"""
n, p = X.shape
beta = np.zeros(p)
eta = X @ beta
mu = link_inv(eta)
for it in range(max_iter):
# 작업 반응과 가중치 (IRLS)
g_prime = link_deriv(mu) # dg/dmu
z = eta + (y - mu) * g_prime # 작업 반응
V_mu = V_fn(mu)
w = 1.0 / (V_mu * g_prime**2) # IRLS 가중치
# 가중 최소제곱
W = np.diag(w)
beta_new = np.linalg.solve(X.T @ W @ X, X.T @ W @ z)
if np.max(np.abs(beta_new - beta)) < tol:
beta = beta_new
break
beta = beta_new
eta = X @ beta
mu = link_inv(eta)
# Pearson 잔차 기반 산포 추정
resid = (y - mu) / np.sqrt(V_fn(mu))
phi = np.sum(resid**2) / (n - p)
cov_beta = phi * np.linalg.inv(X.T @ W @ X)
return beta, phi, cov_beta
# 로짓 링크
link_inv = lambda eta: 1.0 / (1.0 + np.exp(-eta))
link_deriv = lambda mu: 1.0 / (mu * (1.0 - mu))
# 분산함수 mu^2 (1-mu)^2 (Wedderburn 제안)
V_fn = lambda mu: (mu * (1.0 - mu))**2포인트 해설.
- 작업 반응 \(z = \eta + (y-\mu)\,g'(\mu)\) 는 스코어 방정식을 \(\beta\) 에 대해 선형화한 결과이다.
- 가중치 \(w = 1/\{V(\mu)\cdot g'(\mu)^2\}\) 에는 분산함수만 개입한다 — 분포 가정은 전혀 없다.
- \(\phi = X^2/(n-p)\) 로 산포를 추정하고, 공분산에 곱한다.
7.2 프레임워크 — R glm(family=quasi())
# 비율 y, 사이트 site, 품종 variety
# 분산함수 mu(1-mu) 로짓 링크
fit1 <- glm(y ~ site + variety,
family = quasi(link = "logit", variance = "mu(1-mu)"),
data = leaf_blotch)
summary(fit1)
# 분산함수 mu^2 (1-mu)^2 는 내장 지원 없음 → custom family
qvar <- function(mu) (mu * (1 - mu))^2
quasi_custom <- list(
family = "quasi_wedderburn",
link = "logit",
linkfun = function(mu) log(mu/(1-mu)),
linkinv = function(eta) 1/(1+exp(-eta)),
variance = qvar,
dev.resids = function(y, mu, wt) 2*wt*((2*y-1)*log(mu/(1-mu))
- y/mu - (1-y)/(1-mu) + 2),
aic = function(...) NA,
mu.eta = function(eta) {
p <- 1/(1+exp(-eta)); p*(1-p)
},
initialize = expression({mustart <- pmax(pmin(y, 0.999), 0.001)}),
validmu = function(mu) all(mu > 0 & mu < 1)
)
class(quasi_custom) <- "family"
fit2 <- glm(y ~ site + variety, family = quasi_custom, data = leaf_blotch)
summary(fit2)7.3 프레임워크 — Python statsmodels
import statsmodels.api as sm
from statsmodels.genmod.families import Binomial
from statsmodels.genmod.families.varfuncs import VarianceFunction
# 1) 분산함수 mu(1-mu) — 이항과 동일한 준우도
X = sm.add_constant(design) # site + variety 더미
fit1 = sm.GLM(y, X, family=Binomial()).fit(scale='X2') # scale='X2' 로 과산포 추정
print(fit1.summary())
print("phi_hat =", fit1.scale)
# 2) 커스텀 분산함수 mu^2 (1-mu)^2
class VarMuSqOneMinusMuSq(VarianceFunction):
def __call__(self, mu):
return (mu * (1 - mu))**2
def deriv(self, mu):
return 2*mu*(1-mu)*(1 - 2*mu)
fam = Binomial()
fam.variance = VarMuSqOneMinusMuSq()
fit2 = sm.GLM(y, X, family=fam).fit(scale='X2')
print(fit2.summary())결과 해석. scale 속성이 \(\widetilde{\sigma}^2 = X^2/(n-p)\) 이다. 표준오차는 \(\sqrt{\widetilde{\sigma}^2}\) 로 곱해져 이항 모형 대비 확장된다.
8 최적 추정함수와 확장 준우도
8.1 최적 추정함수(Optimal Estimating Function)
준-스코어 \(U(\beta)\) 가 “불편 추정함수” 집합
\[ \mathcal{G} = \{G(\beta; \mathbf{Y}) : E[G] = 0,\ \text{모든 }\beta\text{에서}\} \]
안에서 최적인가를 묻는다. 최적성은 Godambe 정보량
\[ J_G = \{E(\partial G/\partial \beta)\}^T \{\operatorname{cov}(G)\}^{-1} \{E(\partial G/\partial \beta)\} \]
를 최대화하는 \(G\) 로 정의된다. 이는 \(\widehat{\beta}\) 의 점근 분산의 역수에 해당한다.
결과(McCullagh & Nelder, 1989, §9.4). 공분산이 \(\sigma^2 \mathbf{V}(\mu)\) 로 정확히 지정되었다면, 준-스코어 \(U(\beta) = \mathbf{D}^T \mathbf{V}^{-1}(\mathbf{Y}-\boldsymbol{\mu})/\sigma^2\) 는 선형 불편 추정함수 부류 내에서 최적이다. 즉, 같은 1차 모멘트 조건을 만족하는 어떤 선형 함수로도 준-스코어보다 작은 점근 분산을 얻을 수 없다.
직관. 가우스-마코프 정리의 GLM 확장이다. BLUE 가 OLS 가중의 \(\beta\) 와 동일하듯, “best quasi-unbiased estimating function”이 준-스코어와 동일하다.
8.2 확장 준우도(Extended Quasi-likelihood, EQL)
\(\sigma^2\) 를 \(\beta\) 와 함께 우도 방식으로 추정하고 싶을 때 쓰는 조정된 준우도이다.
\[ Q^+(\beta, \phi; \mathbf{y}) = -\tfrac{1}{2}\!\left\{\frac{D(\mathbf{y}; \boldsymbol{\mu})}{\phi} + \log\!\big(2\pi \phi V(y_i)\big)\right\} \]
의 개념으로 (Nelder & Pregibon, 1987) 도입되었다. 이는 \(\mu\) 와 \(\phi\) 를 동시에 최적화할 수 있도록 \(Q\) 를 정규화한 형태이다.
언제 필요한가. 산포 모수 \(\phi\) 자체가 공변량에 따라 변하는 경우 — 예: \(\log \phi_i = z_i^T \gamma\) 를 모형화하는 joint mean-dispersion modelling (Ch.10). 평균 모형과 분산 모형을 번갈아 적합할 때, EQL 이 공통 목적함수 역할을 한다.
주의. \(Q^+\) 는 진짜 로그우도가 아니다. 특히 카운트 자료의 작은 \(\mu\) 영역이나 분산함수가 0 에 접근하는 경계에서 수치적 불안정이 있다. 실무에서는 REML-유사 보정이 함께 쓰인다.
9 이 접근이 없을 때의 문제 — Why It Matters
과산포된 이항·포아송 자료를 이항·포아송 우도로 강행하면 다음이 발생한다.
| 증상 | 원인 | 결과 |
|---|---|---|
| 표준오차 과소 | \(\operatorname{var}(Y) > V(\mu)\) 인데 \(\widetilde{\sigma}^2 = 1\) 강제 | Wald 검정 위신 오해, 거짓 유의 |
| 이탈도 \(\chi^2_{n-p}\) 근사 실패 | 분산이 모형 분산과 다름 | 모형 적합도 검정 왜곡 |
| LRT 기반 모형 선택 불안정 | 우도 자체가 과산포를 반영 못 함 | 불필요한 변수 선택, 과적합 |
준우도 해법.
- 분산함수 \(V(\mu)\) 는 그대로 두되, \(\sigma^2 > 1\) 을 허용한다(준-이항, 준-포아송).
- 점추정은 변하지 않고, 표준오차만 \(\sqrt{\widetilde{\sigma}^2}\) 배 확장된다.
- 모형 비교는 LRT 대신 F-검정 유사 통계량 \((\Delta D/\Delta\text{df})/\widetilde{\sigma}^2\) 을 쓴다.
R::glm(family=quasibinomial) 이나 family=quasipoisson 이 정확히 이 레시피다.
10 응용 분야
| 분야 | 활용 | 구체 예시 |
|---|---|---|
| 보건·역학 | 군집 샘플링 이항 자료 | 가구 단위로 묶인 개인 질병 유무 — 준-이항 로짓 |
| 생태학 | 과산포 카운트 | 덫당 포획 개체수 — 준-포아송 또는 음이항 |
| 마케팅 | 전환율 자료 | 세그먼트별 전환율 (연속 비율) — \(V(\pi) = \pi(1-\pi)\) 준우도 |
| 임상 시험 | 반복 측정 이진 반응 | 준우도 GEE, 교환가능/AR(1) 작업 상관 |
| 제조 품질 | 불량률의 분산-평균 교정 | 극단 \(\mu\) 구간에서 \(\mu^2(1-\mu)^2\) 분산함수 |
| 보험 | 청구건수 과산포 | 준-포아송 또는 음이항 |
공통 구조. 반응의 분포는 복잡·미상이거나 과산포된 지수족 근사이지만, 평균-분산 관계의 질적 형태는 과거 경험으로부터 알 수 있다. 준우도는 이 최소 정보만으로 GLM 기계를 돌린다.