Independent Observations — 준우도 상세

1 도입 — §9.2 의 문제 설정

Ch.9 개관(08-1)에서 준우도의 핵심 아이디어를 소개했다. 이 포스트는 독립 관측(§9.2) 상황에서의 준우도 이론을 처음 원리부터 상세히 전개한다. 다음을 다룬다.

• 공분산 함수의 형태적 제약 — 왜 대각이어야 하며, 각 \(V_i\) 가 자기 \(\mu_i\) 에만 의존해야 하는가.
• 준-스코어의 세 성질을 직접 유도하고, 로그우도 미분과 일대일 대응을 확인한다.
• Table 9.1 의 모든 분산함수에 대한 준우도 \(Q(\mu; y)\) 를 단계별로 적분한다.
• Newton-Raphson·Fisher 스코어링의 대수를 전개하고, IRLS 의 작업 반응·가중치가 어디서 나오는지 밝힌다.
• 1단계 추정량의 선형화를 통해 비편향성·점근 정규성을 증명한다.
• \(\sigma^2\) 의 Pearson 기반 추정의 논리적 근거.
• Wedderburn (1974)의 보리 잎 반점 자료를 두 분산함수로 비교 적합한 완전한 수치 결과.

전제. 관측 \(\mathbf{Y} = (Y_1, \ldots, Y_n)^T\) 은 서로 독립이며, 평균 벡터 \(\boldsymbol{\mu}(\beta)\) 는 회귀 모수 \(\beta \in \mathbb{R}^p\) 의 매끄러운 함수이다. 공분산은

\[ \operatorname{cov}(\mathbf{Y}) = \sigma^2 \mathbf{V}(\boldsymbol{\mu}) \]

이고, \(\sigma^2\) 는 (미지의) 스칼라 산포 모수이며 \(\beta\) 에 의존하지 않는다.

2 공분산 함수의 형태적 제약

2.1 대각 조건

관측이 독립이므로 공분산 행렬은 대각이다.

\[ \mathbf{V}(\boldsymbol{\mu}) = \operatorname{diag}\{V_1(\boldsymbol{\mu}), \ldots, V_n(\boldsymbol{\mu})\}. \]

이는 자료의 성질에서 자연스럽게 따라 나오는 결과이지 선택이 아니다.

2.2 함수적 독립 조건(식 9.1)

더 강한 조건이 부과된다. 각 \(V_i(\boldsymbol{\mu})\) 는 오직 자기 성분 \(\mu_i\) 에만 의존해야 한다.

\[ \mathbf{V}(\boldsymbol{\mu}) = \operatorname{diag}\{V_1(\mu_1), \ldots, V_n(\mu_n)\}. \tag{9.1} \]

왜 이 조건이 필요한가. 원리적으로는 \(V_1\) 이 \((\mu_1, \mu_2, \ldots)\) 모두에 의존할 수도 있다. 그러나 McCullagh & Nelder 는 두 가지 이유로 이를 배제한다.

1. 물리적 비합리성. 두 관측이 독립인데 한쪽의 분산이 다른 쪽의 평균에 의존하는 기제는 상상하기 어렵다. 독립과 함수 의존 양립은 예외적 구성이다.
2. 기술적(수학적) 필요. §9.2.2 에서 스칼라 함수 \(Q(\mu; y)\) 를 적분으로 정의하려면 피적분 함수가 1차원 경로에서만 의미를 가져야 한다. 함수적 독립이 성립해야 \(Q_i\) 를 각 성분별로 독립적으로 계산하고 합할 수 있다.

반사실 — 조건이 없으면 어떻게 되는가. 만일 \(V_1\) 이 \(\mu_2\) 에도 의존한다면, 선적분 \(\int_{\boldsymbol{\mu}_0}^{\boldsymbol{\mu}} (\mathbf y - \mathbf t)^T \mathbf V(\mathbf t)^{-1} d\mathbf t\) 의 값이 선택한 경로에 따라 달라진다. 그 결과 \(Q(\boldsymbol\mu; \mathbf y)\) 는 다값함수(multi-valued)가 되어 하나의 수치로 정의되지 않는다. 이때는 §9.3 의 종속 관측 틀 — 적분 대신 준-스코어 방정식 자체를 다변수로 세우는 방식 — 으로 넘어가야 한다.

so what. 대부분의 실제 응용에서 \(V_1 = V_2 = \cdots = V\) 라는 동일 분산함수 가정이 성립한다. 이 장의 대수는 그렇지 않아도 유효하지만, 표기의 단순화를 위해 자주 \(V(\mu_i)\) 로 쓴다.

2.3 표기 통일

이 포스트 전체에서 다음 약속을 쓴다.

기호	의미
\(\mu_i = E(Y_i)\)	\(i\)-번째 관측의 평균
\(V(\mu_i)\)	\(i\)-번째 관측의 분산함수(평균의 함수)
\(\operatorname{var}(Y_i) = \sigma^2 V(\mu_i)\)	관측 분산
\(\eta_i = x_i^T \beta\)	선형예측자
\(g(\mu_i) = \eta_i\)	링크 함수
\(\mathbf{D}\)	\(n \times p\) 행렬, \(D_{ir} = \partial \mu_i / \partial \beta_r\)
\(\mathbf{V}\)	\(n \times n\) 대각, \(V_{ii} = V(\mu_i)\)

3 준-스코어의 세 성질 — 직접 유도

3.1 단일 관측의 준-스코어

단일 관측 \(Y\) (평균 \(\mu\), 분산 \(\sigma^2 V(\mu)\)) 에 대해

\[ U = u(\mu; Y) = \frac{Y - \mu}{\sigma^2 V(\mu)}. \]

이 함수의 첫 두 모멘트와 부분도함수의 기댓값을 차례로 계산한다.

3.2 성질 1: \(E(U) = 0\)

\[ E(U) = \frac{1}{\sigma^2 V(\mu)}\, E(Y - \mu) = \frac{1}{\sigma^2 V(\mu)} \cdot 0 = 0. \]

이유. \(E(Y) = \mu\) 자체가 가정이므로 자명하다. 분산함수와 \(\sigma^2\) 은 상수처럼 빠져나온다. 이 성질은 정규 방정식의 불편성을 보장한다 — 추정방정식 \(U(\widehat{\beta}) = 0\) 의 해가 참값 근처에 머문다.

3.3 성질 2: \(\operatorname{var}(U) = 1/\{\sigma^2 V(\mu)\}\)

\[ \operatorname{var}(U) = \frac{1}{\{\sigma^2 V(\mu)\}^2}\, \operatorname{var}(Y - \mu) = \frac{\sigma^2 V(\mu)}{\{\sigma^2 V(\mu)\}^2} = \frac{1}{\sigma^2 V(\mu)}. \]

해석. \(U\) 의 분산은 관측 \(Y\) 의 분산의 역수에 \(\sigma^2\) 을 뺀 형태로 표현된다. 분산이 큰 관측은 \(U\) 가 0 근처에 머물고, 분산이 작은 관측은 \(U\) 가 크게 흔들린다 — 즉, 분산이 작을수록 스코어의 신호가 강하다. IRLS 가중치가 \(1/V(\mu)\) 로 나오는 근본 이유이다.

3.4 성질 3: \(-E(\partial U/\partial \mu) = 1/\{\sigma^2 V(\mu)\}\)

편미분을 계산한다.

\[ \frac{\partial U}{\partial \mu} = \frac{\partial}{\partial \mu}\!\left[\frac{Y - \mu}{\sigma^2 V(\mu)}\right] = \frac{-\sigma^2 V(\mu) - (Y-\mu)\sigma^2 V'(\mu)}{\{\sigma^2 V(\mu)\}^2}. \]

기댓값을 취하면 \(E(Y-\mu) = 0\) 이므로 두 번째 항이 사라지고,

\[ E\!\left(\frac{\partial U}{\partial \mu}\right) = \frac{-1}{\sigma^2 V(\mu)}. \]

따라서

\[ - E\!\left(\frac{\partial U}{\partial \mu}\right) = \frac{1}{\sigma^2 V(\mu)} = \operatorname{var}(U). \]

핵심 관찰. 분산 = \(-\) 기댓값 미분이라는 등식이 성립한다. 이는 로그우도 미분이 만족하는 Bartlett 의 제2 항등식과 동일한 구조다. 준-스코어는 이 관점에서 “2차 Bartlett 항등식을 만족하는 비-우도 함수”라고 볼 수 있다.

정리: 준-스코어의 세 성질

관측 \(Y\) 가 \(E(Y) = \mu\), \(\operatorname{var}(Y) = \sigma^2 V(\mu)\) 를 만족하면 \(U = (Y-\mu)/\{\sigma^2 V(\mu)\}\) 는

\[ 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)} \]

을 만족한다. 세 성질은 로그우도 미분(스코어)이 만족하는 것과 정확히 같다.

3.5 왜 이 성질이 “추론의 거의 전부”인가

1차 점근 이론에서 스코어의 역할은 다음 세 가지다.

• 일치성. \(E(U) = 0\) 이 \(U(\widehat{\theta}) = 0\) 의 해 \(\widehat{\theta}\) 가 참값에 수렴하도록 한다.
• 점근 분산. 점근 분산이 Fisher 정보 \(-E(\partial U/\partial \theta)\) 의 역수로 나온다.
• 표본 간 집계. \(n\) 개 관측의 스코어 합이 CLT 에 의해 정규 근사된다.

이 세 역할은 모두 위의 세 성질만으로 충족된다. 분포 전체의 정보 \(f(y;\theta)\) 는 필요하지 않다. 이것이 준우도가 가능한 이유이다.

4 준우도 함수의 정의

4.1 적분으로서의 \(Q\)

단일 관측의 준우도는

\[ Q(\mu; y) = \int_{y}^{\mu} \frac{y - t}{\sigma^2 V(t)}\,dt. \tag{9.3} \]

적분의 의미. 피적분 함수 \((y-t)/\{\sigma^2 V(t)\}\) 는 “현 위치 \(t\) 에서의 준-스코어”이다. 하한 \(y\) 에서 시작해 상한 \(\mu\) 까지 적분하면, \(\mu\) 에서의 누적 준우도가 된다. 미분하면

\[ \frac{\partial Q}{\partial \mu} = \frac{y - \mu}{\sigma^2 V(\mu)} = U, \]

즉 \(Q\) 를 \(\mu\) 에 대해 미분하면 준-스코어가 복원된다. 이는 로그우도와 스코어의 관계와 같다.

하한을 왜 \(y\) 로 잡나. 하한의 선택은 \(\mu\) 에 무관한 상수만 더해준다. \(y\) 를 택하면 \(Q(\mu = y; y) = 0\) 이 되어 “데이터가 정확히 관측값을 맞추면 준우도가 0”이라는 자연스러운 기준점이 생긴다.

4.2 전체 표본의 준우도

관측이 독립이므로

\[ Q(\boldsymbol{\mu}; \mathbf{y}) = \sum_{i=1}^n Q_i(\mu_i; y_i). \]

\(\mathbf{y}\) 에 대한 적분은 필요 없다. 독립성이 이 단순한 합산을 가능하게 한다.

5 Table 9.1 — 모든 항목 유도

Table 9.1 의 각 분산함수에 대해 \(Q(\mu; y)\) 를 직접 계산한다. \(\sigma^2\) 는 상수이므로 \(\sigma^{-2}\) 을 흡수하여 다음 적분을 계산한다.

\[ Q^*(\mu; y) = \int_y^\mu \frac{y-t}{V(t)}\,dt. \]

(실제 \(Q = Q^*/\sigma^2\) 이다. 표에서는 이 정규화가 관례적으로 생략된다.)

5.1 \(V(\mu) = 1\) — Normal

\[ Q^* = \int_y^\mu (y - t)\,dt = \left[yt - \tfrac{t^2}{2}\right]_y^\mu = y\mu - \tfrac{\mu^2}{2} - \tfrac{y^2}{2} = -\tfrac{(y-\mu)^2}{2}. \]

확인. 정규분포 로그우도 \(-\tfrac{1}{2}\log(2\pi\sigma^2) - \tfrac{(y-\mu)^2}{2\sigma^2}\) 의 \(\mu\)-의존 부분과 일치한다. 정준모수 \(\theta = \mu\).

5.2 \(V(\mu) = \mu\) — Poisson

\[ Q^* = \int_y^\mu \frac{y-t}{t}\,dt = \int_y^\mu\!\left(\frac{y}{t} - 1\right)dt = y\log(\mu/y) - (\mu - y). \]

상수 항(데이터만의 함수)을 \(y\) 관련 상수로 묶어내면 \(\mu\)-의존 부분은

\[ Q^*(\mu; y) = y\log\mu - \mu. \]

확인. 포아송 로그우도 \(y\log\mu - \mu - \log(y!)\) 의 \(\mu\)-의존 부분과 일치. 정준모수 \(\theta = \log\mu\).

5.3 \(V(\mu) = \mu^2\) — Gamma

\[ Q^* = \int_y^\mu \frac{y-t}{t^2}\,dt = \int_y^\mu\!\left(\frac{y}{t^2} - \frac{1}{t}\right)dt = -\frac{y}{\mu} - \log\mu + \frac{y}{y} + \log y. \]

\(\mu\)-의존 부분은

\[ Q^*(\mu; y) = -\frac{y}{\mu} - \log\mu. \]

확인. 감마 로그우도의 \(\mu\)-의존 부분(\(\nu\) 로 재모수화 시). 정준모수 \(\theta = -1/\mu\).

5.4 \(V(\mu) = \mu^3\) — Inverse Gaussian

\[ Q^* = \int_y^\mu \frac{y-t}{t^3}\,dt = \int_y^\mu\!\left(\frac{y}{t^3} - \frac{1}{t^2}\right)dt = -\frac{y}{2\mu^2} + \frac{1}{\mu} + \frac{y}{2y^2} - \frac{1}{y}. \]

\(\mu\)-의존 부분:

\[ Q^*(\mu; y) = -\frac{y}{2\mu^2} + \frac{1}{\mu}. \]

정준모수 \(\theta = -1/(2\mu^2)\).

5.5 \(V(\mu) = \mu^\zeta\), \(\zeta \neq 0, 1, 2\) — Tweedie 계열

\[ Q^* = \int_y^\mu\!\left(\frac{y}{t^\zeta} - t^{1-\zeta}\right)dt = y\cdot\frac{t^{1-\zeta}}{1-\zeta}\bigg|_y^\mu - \frac{t^{2-\zeta}}{2-\zeta}\bigg|_y^\mu. \]

\(\mu\)-의존 부분은

\[ Q^*(\mu; y) = \frac{y\mu^{1-\zeta}}{1-\zeta} - \frac{\mu^{2-\zeta}}{2-\zeta} = \mu^{-\zeta}\!\left(\frac{\mu y}{1-\zeta} - \frac{\mu^2}{2-\zeta}\right). \]

의미. Tweedie 분포족의 분산함수가 \(V(\mu) = \phi\mu^p\) 이며 \(p = 2\) 에서 감마, \(p = 1\) 에서 포아송, \(1 < p < 2\) 에서 복합 포아송-감마(claims 모델)가 된다.

5.6 \(V(\mu) = \mu(1-\mu)\) — Binomial/\(m\)

\[ Q^* = \int_y^\mu \frac{y-t}{t(1-t)}\,dt. \]

부분분수 분해 \(\frac{1}{t(1-t)} = \frac{1}{t} + \frac{1}{1-t}\) 를 쓰면 피적분 함수는

\[ \frac{y-t}{t(1-t)} = y\!\left(\frac{1}{t} + \frac{1}{1-t}\right) - \frac{t}{t(1-t)}\cdot(t) \quad\text{(재정리)}. \]

간단히 정리하면

\[ \frac{y-t}{t(1-t)} = \frac{y}{t} - \frac{1-y}{1-t} \]

(좌우 \(t(1-t)\) 곱해 검산: \(y(1-t) - (1-y)t = y - yt - t + yt = y - t\). 맞음.)

적분하면

\[ \int \frac{y}{t} dt - \int \frac{1-y}{1-t}dt = y\log t + (1-y)\log(1-t) + C. \]

\(\mu\)-의존 부분:

\[ Q^*(\mu; y) = y\log\mu + (1-y)\log(1-\mu) = y\log\!\left(\frac{\mu}{1-\mu}\right) + \log(1-\mu). \]

확인. 이항 로그우도의 형태와 동일. 정준모수는 로짓 \(\theta = \log\{\mu/(1-\mu)\}\). 중요한 점은 \(y\) 가 반드시 \(0/1\) 일 필요가 없고, \([0, 1]\) 의 연속 비율이어도 준우도로 유효하다는 것이다.

5.7 \(V(\mu) = \mu^2(1-\mu)^2\) — 분포 없음

이 분산함수는 알려진 분포에 대응하지 않는다. 부분분수 분해를 먼저 한다.

\[ \frac{1}{t^2(1-t)^2} = \frac{A}{t} + \frac{B}{t^2} + \frac{C}{1-t} + \frac{D}{(1-t)^2}. \]

계수 비교로 \(A = 2, B = 1, C = 2, D = 1\) 를 얻는다. 피적분 함수는

\[ \frac{y-t}{t^2(1-t)^2} = (y-t)\!\left[\frac{2}{t} + \frac{1}{t^2} + \frac{2}{1-t} + \frac{1}{(1-t)^2}\right]. \]

적분 후 \(\mu\)-의존 부분을 정리하면 (대수가 길어 결과만 제시)

\[ Q^*(\mu; y) = (2y - 1)\log\!\left(\frac{\mu}{1-\mu}\right) - \frac{y}{\mu} - \frac{1-y}{1-\mu}. \]

주의. 이 함수는 \(\mu \in (0, 1)\) 에서만 정의되며 경계 \(0, 1\) 에서 발산한다. 관측된 비율 \(y_i = 0\) 이 있으면 이탈도를 보통의 방식으로 정의하기 어렵다(Wedderburn의 관찰).

직관. \(\mu^2(1-\mu)^2 = [\mu(1-\mu)]^2\) 는 이항 분산의 제곱이다. 즉 극단 \(\mu\) 에서 분산이 더 빨리 작아진다. 보리 잎 반점처럼 극단 비율에서 변동이 과도하게 감소하는 자료에 적합하다.

5.8 \(V(\mu) = \mu + \mu^2/k\) — Negative Binomial

\(k\) 가 고정된 경우에 한해 이는 음이항 분산함수이다. 부분분수

\[ \frac{1}{t(1 + t/k)} = \frac{1}{t} - \frac{1}{k + t} \]

를 이용해

\[ Q^*(\mu; y) = y\log\!\left(\frac{\mu}{k+\mu}\right) + k\log\!\left(\frac{k}{k+\mu}\right) \]

를 얻는다. 정준모수 \(\theta = \log\{\mu/(k+\mu)\}\).

실무 주의. \(k\) 가 알려지지 않으면 이는 완전한 준우도가 아니다. 별도 추정(예: Pearson 모멘트)이 필요하다.

6 준-이탈도

\(\sigma^2\) 의존성을 제거하기 위해 이탈도를 다음과 같이 정의한다.

\[ D(y; \mu) = -2\sigma^2 Q(\mu; y) = 2\int_\mu^y \frac{y - t}{V(t)}\,dt. \tag{9.4} \]

6.1 핵심 성질

성질	이유
\(D(y; y) = 0\)	적분 구간이 한 점이므로
\(D(y; \mu) \geq 0\)	피적분 함수의 부호가 \((y-t)\) 와 일치
\(D\) 는 \(\sigma^2\) 에 의존하지 않음	정의에서 \(\sigma^2\) 가 상쇄
\(D\) 는 \(y\) 와 \(\mu\) 의 순수 함수	관측과 적합값만으로 계산 가능

6.2 예제: 포아송 이탈도

\(V(t) = t\) 이므로

\[ D(y;\mu) = 2\int_\mu^y \frac{y-t}{t}\,dt = 2[y\log t - t]_\mu^y = 2\!\left[y\log(y/\mu) - (y-\mu)\right]. \]

GLM 의 포아송 이탈도와 정확히 일치한다. 준우도의 이탈도가 지수족 이탈도와 동일한 대수적 형태로 나옴을 확인할 수 있다.

6.3 예제: 이항 이탈도

\(V(t) = t(1-t)\) 이므로

\[ D(y;\mu) = 2\!\left[y\log(y/\mu) + (1-y)\log\!\left(\frac{1-y}{1-\mu}\right)\right]. \]

6.4 전체 이탈도

\(D(\mathbf{y}; \boldsymbol{\mu}) = \sum_i D(y_i; \mu_i)\). 모형 비교에서

\[ F = \frac{\{D(\mathbf{y}; \widehat{\boldsymbol{\mu}}_0) - D(\mathbf{y}; \widehat{\boldsymbol{\mu}}_1)\}/(p_1 - p_0)}{\widetilde{\sigma}^2} \]

는 근사 \(F_{p_1-p_0,\, n-p_1}\) 분포를 따른다. LRT 대신 이 \(F\)-검정 통계량이 과산포 상황에서 일관된 모형 비교를 제공한다.

7 파라미터 추정 — Newton-Raphson·Fisher 스코어링

7.1 준-스코어 방정식

회귀 모수 \(\beta\) 에 대한 총 준-스코어는 연쇄법칙에 의해

\[ U_r(\beta) = \sum_{i=1}^n \frac{Y_i - \mu_i}{\sigma^2 V(\mu_i)}\cdot\frac{\partial \mu_i}{\partial \beta_r} = \sum_i \frac{(Y_i - \mu_i)D_{ir}}{\sigma^2 V(\mu_i)}. \]

행렬 표현으로

\[ U(\beta) = \frac{1}{\sigma^2}\mathbf{D}^T \mathbf{V}^{-1}(\mathbf{Y} - \boldsymbol{\mu}). \tag{9.5} \]

추정방정식은 \(U(\widehat{\beta}) = \mathbf{0}\), 즉 \(\mathbf{D}^T \mathbf{V}^{-1}(\mathbf{Y} - \boldsymbol{\mu}) = \mathbf{0}\). \(\sigma^2\) 은 양변에서 상쇄되어 점추정에는 영향을 주지 않는다.

7.2 준-정보 행렬

세 성질 중 성질 2 와 3 에 의해 \(\operatorname{cov}(U) = -E(\partial U/\partial \beta)\) 이다. 계산하면

\[ \mathbf{i}_\beta = \operatorname{cov}(U) = \frac{1}{\sigma^2}\mathbf{D}^T \mathbf{V}^{-1}\mathbf{D}. \tag{9.6} \]

Fisher 정보의 대응. 일반 MLE 에서 Fisher 정보는 스코어의 공분산이다. 준우도에서는 준-스코어의 공분산이 그 자리를 대신한다.

7.3 Fisher 스코어링 반복

Newton-Raphson 은 해시안 \(\partial U/\partial \beta\) 를 쓰지만, 독립 관측 하에서는 기댓값을 쓴 Fisher 스코어링이 더 안정적이다.

\[ \widehat{\beta}_{k+1} = \widehat{\beta}_k + \mathbf{i}_\beta^{-1}(\widehat{\beta}_k)\cdot U(\widehat{\beta}_k) = \widehat{\beta}_k + (\mathbf{D}_k^T \mathbf{V}_k^{-1}\mathbf{D}_k)^{-1}\mathbf{D}_k^T \mathbf{V}_k^{-1}(\mathbf{y} - \boldsymbol{\mu}_k). \]

\(\mathbf{D}_k\) 와 \(\mathbf{V}_k\) 는 \(\widehat{\beta}_k\) 에서 평가된 값이다. \(\sigma^2\) 은 \(\mathbf{V}^{-1}\) 와 \(\mathbf{i}_\beta^{-1}\) 에서 상쇄된다.

7.4 IRLS 형태 — 작업 반응과 가중치의 유도

스코어 방정식 \(\mathbf{D}^T \mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) = 0\) 을 \(\beta\) 에서 선형화한다. \(\mu_i = g^{-1}(\eta_i)\), \(\eta_i = x_i^T\beta\) 이므로

\[ \frac{\partial \mu_i}{\partial \beta_r} = \frac{d\mu_i}{d\eta_i}\cdot x_{ir} = \frac{1}{g'(\mu_i)}\cdot x_{ir}. \]

즉 \(\mathbf{D} = \mathbf{W}_0\mathbf{X}\), 여기서 \(\mathbf{W}_0 = \operatorname{diag}\{1/g'(\mu_i)\}\). 대입하면

\[ \mathbf{X}^T \mathbf{W}_0 \mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) = 0. \]

선형화: \(g(\mathbf{y}) \approx g(\boldsymbol{\mu}) + g'(\boldsymbol{\mu})\cdot(\mathbf{y} - \boldsymbol{\mu})\) 이므로 작업 반응

\[ z_i = \eta_i + (y_i - \mu_i)g'(\mu_i) \]

을 정의하면 \((y_i - \mu_i) = (z_i - \eta_i)/g'(\mu_i)\). 대입하여

\[ \mathbf{X}^T \mathbf{W}_0 \mathbf{V}^{-1}\!\left[\frac{1}{g'(\boldsymbol{\mu})}(\mathbf{z} - \boldsymbol{\eta})\right] = 0 \;\Longrightarrow\; \mathbf{X}^T \mathbf{W}(\mathbf{z} - \mathbf{X}\beta) = 0. \]

IRLS 가중치는

\[ W_{ii} = \frac{1}{V(\mu_i)[g'(\mu_i)]^2}. \]

이 방정식은 작업 반응 \(\mathbf{z}\) 에 대한 가중 최소제곱 문제 \(\mathbf{X}^T \mathbf{W}\mathbf{X}\beta = \mathbf{X}^T \mathbf{W}\mathbf{z}\) 와 동일하다. 매 반복마다 \(\mu_i\), \(z_i\), \(W_{ii}\) 를 갱신하고 가중 최소제곱을 푼다.

핵심 관찰. 가중치와 작업 반응이 오직 \(V(\mu)\) 와 \(g(\mu)\) 만으로 결정된다. 분포 가정은 어디에도 등장하지 않는다. 이것이 준우도 IRLS 가 GLM IRLS 와 완전히 동일한 엔진을 쓰는 이유다.

7.5 정준 링크 하에서의 단순화

\(g = (b^*)'\) 인 정준 링크에서는 \(g'(\mu) = 1/V(\mu)\) 가 성립한다(지수족의 기본 항등식). 대입하면

\[ W_{ii} = \frac{1}{V(\mu_i)\cdot V(\mu_i)^{-2}} = V(\mu_i), \]

\[ z_i = \eta_i + (y_i - \mu_i)/V(\mu_i). \]

가중치와 분산함수가 정확히 같은 형태로 나타나면서, 식이 극적으로 간소화된다. 이것이 “정준 링크에서 GLM 이 아름답다”는 말의 실체다.

8 1단계 추정량의 점근 성질

참값 \(\beta\) 에서 시작한 Fisher 스코어링 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} \]

\(\widehat{\beta}_1\) 은 계산 가능한 통계량은 아니지만(참값을 모르므로) 이론적 분석에 유용하다.

8.1 비편향성

\(E(\mathbf{y} - \boldsymbol{\mu}) = \mathbf{0}\) 이므로

\[ E(\widehat{\beta}_1) = \beta + (\mathbf{D}^T \mathbf{V}^{-1}\mathbf{D})^{-1}\mathbf{D}^T \mathbf{V}^{-1}\cdot\mathbf{0} = \beta. \]

단, 이는 1단계에 한한 엄밀한 비편향이다. 수렴한 \(\widehat{\beta}\) 는 일반적으로 \(O(n^{-1})\) 편향을 갖는다(Cox & Snell 유형 보정 필요).

8.2 점근 분산

\[ \operatorname{var}(\widehat{\beta}_1) = (\mathbf{D}^T \mathbf{V}^{-1}\mathbf{D})^{-1}\mathbf{D}^T \mathbf{V}^{-1}\cdot \operatorname{var}(\mathbf{y})\cdot \mathbf{V}^{-1}\mathbf{D}(\mathbf{D}^T \mathbf{V}^{-1}\mathbf{D})^{-1}. \]

\(\operatorname{var}(\mathbf{y}) = \sigma^2 \mathbf{V}\) 를 대입하면 중간 항이 \(\sigma^2\mathbf{V}^{-1}\mathbf{V}\mathbf{V}^{-1} = \sigma^2\mathbf{V}^{-1}\) 로 축약되어

\[ \operatorname{var}(\widehat{\beta}_1) = \sigma^2(\mathbf{D}^T \mathbf{V}^{-1}\mathbf{D})^{-1} = \mathbf{i}_\beta^{-1}. \]

8.3 점근 정규성

\(\widehat{\beta}_1\) 은 \((\mathbf{y} - \boldsymbol{\mu})\) 의 선형 결합이고, 각 \((y_i - \mu_i)\) 는 평균 0 분산 \(\sigma^2 V(\mu_i)\) 의 독립 확률변수이다. 적절한 Lindeberg 조건 하에 CLT 가 적용되어

\[ \widehat{\beta}_1 \stackrel{d}{\to} N(\beta,\; \sigma^2(\mathbf{D}^T \mathbf{V}^{-1}\mathbf{D})^{-1}). \]

so what. 오직 처음 두 모멘트만 가정하고도 분포 수준의 점근 결과를 얻었다. 완전한 우도 함수 없이 추론·검정이 가능한 정확한 이유다.

9 \(\sigma^2\) 의 추정 — Pearson 통계량

9.1 추정량

일반화 Pearson 통계량

\[ X^2 = \sum_i \frac{(Y_i - \widehat{\mu}_i)^2}{V_i(\widehat{\mu}_i)} \]

을 이용하여

\[ \widetilde{\sigma}^2 = \frac{X^2}{n - p}. \]

9.2 왜 이탈도 기반이 아닌가

일반 GLM 에서는 이탈도 \(D\) 가 근사적으로 \(\sigma^2 \chi^2_{n-p}\) 분포를 따르므로 \(D/(n-p)\) 도 \(\sigma^2\) 의 추정량이다. 준우도에서는 이 근사가 덜 안정적이다.

• \(D\) 의 카이제곱 근사는 지수족 분포를 전제로 유도되는데, 준우도는 분포를 가정하지 않는다.
• Pearson \(X^2\) 은 2차 모멘트만 쓰므로 분포 가정에 덜 민감하다.
• 실증적으로 \(X^2/(n-p)\) 가 \(D/(n-p)\) 보다 강건하게 \(\sigma^2\) 을 추정한다(McCullagh & Nelder, 1989, §9.2.3).

9.3 분자유도에 대한 유의점

\(n - p\) 는 잔차 자유도이다. 이는 추정된 \(p\) 개 회귀 모수가 잔차 제곱의 유효 자유도를 \(p\) 만큼 감소시킨다는 근사에 기반한다. 비선형 모형이나 제약이 있는 경우에는 별도 조정이 필요하다.

10 예제 — Leaf-blotch 완전 분석

10.1 자료

Wedderburn (1974) 이 보고한 보리 잎 반점(Rhynchosporium secalis) 감염 자료. 9 개 사이트 × 10 품종 = 90 관측. 반응 \(y\) 는 감염된 잎 면적의 비율(\([0, 1]\) 연속). 모든 관측이 독립(서로 다른 시험구)이라 가정한다.

10.2 모형

\[ \operatorname{logit}(\mu_{ij}) = \alpha_i + \gamma_j,\qquad \operatorname{var}(Y_{ij}) = \sigma^2 V(\mu_{ij}). \]

주효과만(사이트 9개, 품종 10개), 교호작용 없음. 잔차 자유도 \(n - p = 90 - (1 + 8 + 9) = 72\).

10.3 적합 1 — \(V(\mu) = \mu(1-\mu)\)

이항 분산함수로 family=quasi(link="logit", variance="mu(1-mu)").

통계량	값
잔차 이탈도	\(D = 6.13\)
Pearson 통계량	\(X^2 = 6.39\)
잔차 자유도	\(n - p = 72\)
\(\widetilde{\sigma}^2 = X^2/(n-p)\)	\(0.089\)

해석. \(\widetilde{\sigma}^2 = 0.089 \ll 1\) 이다. 이는 자료가 이항 기댓값보다 훨씬 작은 변동을 보인다는 뜻이다. 카운트 자료가 아니므로 \(\sigma^2 \approx 1\) 을 기대할 이유는 없다. 문제는 \(\sigma^2\) 의 크기 자체가 아니라, 분산-평균 관계가 \(\mu(1-\mu)\) 로 잘 포착되는가이다.

품종 효과(변종 1 기준 contrast, logit 스케일).

품종	1	2	3	4	5	6	7	8	9	10
\(\widehat{\gamma}_j\)	0.00	0.15	0.69	1.05	1.62	2.37	2.57	3.34	3.50	4.25
SE	0.00	0.72	0.67	0.65	0.63	0.61	0.61	0.60	0.60	0.60

품종 1–3 이 저항성이 가장 높고, 품종 8–10 이 가장 취약하다.

진단 그림 요약(Fig. 9.1 of McCullagh & Nelder).

• Pearson 잔차 \(r_P = (y - \widehat{\mu})/\sqrt{\widehat{\mu}(1-\widehat{\mu})}\) 를 \(\widehat{\eta}\) 에 대해 플롯.
• 양 극단(\(\widehat{\eta}\) 이 매우 크거나 작을 때)에서 잔차가 계통적으로 크다.
• \(\mu(1-\mu)\) 가 \(\mu \to 0, 1\) 에서 너무 빨리 0 으로 떨어져, 자료의 실제 산포를 과소평가한다.

10.4 적합 2 — \(V(\mu) = \mu^2(1-\mu)^2\)

Wedderburn 의 제안: 극단에서 분산이 더 빨리 감소하는 함수를 쓴다. \(Q(\mu; y) = (2y-1)\log\{\mu/(1-\mu)\} - y/\mu - (1-y)/(1-\mu)\).

통계량	값
Pearson 가중 잔차 제곱합	\(71.2\)
잔차 자유도	\(72\)
\(\widetilde{\sigma}^2\)	\(0.99\)

해석. 이제 \(\widetilde{\sigma}^2 \approx 1\) 으로, 분산 구조가 자료에 잘 부합한다. 이탈도는 \(y = 0\) 관측들로 인해 표준 방식으로 정의되지 않으므로 Pearson 쪽에 의존한다.

품종 효과 (재적합).

품종	1	2	3	4	5	6	7	8	9	10
\(\widehat{\gamma}_j\)	0.00	0.47	0.08	0.95	1.35	1.33	2.34	3.26	3.14	3.89
SE	0.00	0.47	0.47	0.47	0.47	0.47	0.47	0.47	0.47	0.47

관찰. 모든 품종 대비의 표준오차가 정확히 같다. 이는 IRLS 가중치가 상수 1 이 되어 분석이 직교화되기 때문이다(정준 링크+ 이 분산함수 조합의 특이성).

• 저 감염률 품종의 상대 정밀도 증가: 이전 모델에서 SE 0.72 → 이제 0.47.
• 고 감염률 품종의 상대 정밀도 감소: 이전 0.60 → 이제 0.47.
• 전체 표준오차의 평균 수준은 유지되지만, 분포가 달라진다.

품종 순위 변화. 첫 모델에서 품종 3 이 중간(0.69)이었는데 새 모델에서는 최하위에 가까움(0.08). 자료의 실제 분산 구조를 정확히 반영하면 품종 간 서열이 바뀔 수 있다는 실무적 경고.

진단 요약(Fig. 9.2).

• \(\widehat{\eta}\) 에 대한 Pearson 잔차 패턴이 현저히 개선됨.
• 잔여 편의 3 개 이상치: (품종 4, 사이트 3) = 3.01, (품종 5, 사이트 7) = 2.51, (품종 6, 사이트 8) = 2.24.
• 전반적 구조적 이탈 없음.

10.5 핵심 메시지

1. 분포를 바꾸지 않고 분산함수만 교체해도 결과가 크게 달라질 수 있다.
2. 잔차 진단은 분산함수 설정의 필수 도구이다.
3. 분산함수가 자료와 맞지 않으면, 점추정은 그럭저럭이어도 표준오차·검정 통계량이 체계적으로 왜곡된다.
4. 준우도는 이 모든 진단·교정을 분포 가정 없이 허용한다.

11 코드 예시

11.1 Python 저수준 — IRLS 직접 구현

import numpy as np

def quasi_irls(X, y, V_fn, link_inv, link_deriv, tol=1e-8, max_iter=100):
    """
    준우도 IRLS.
    X: (n, p) 설계행렬
    y: (n,) 반응
    V_fn(mu): 분산함수
    link_inv(eta): g^{-1}
    link_deriv(mu): g'(mu) = d eta / d mu
    """
    n, p = X.shape
    # 초기화: 링크의 정의역에 맞는 보수적 시작값
    mu = np.clip(y, 0.001, 0.999) if y.max() <= 1 else np.maximum(y, 0.1)
    eta = np.log(mu / (1 - mu))  # 로짓 가정 (일반화 시 수정)

    for it in range(max_iter):
        g_prime = link_deriv(mu)
        z = eta + (y - mu) * g_prime      # 작업 반응
        V = V_fn(mu)
        w = 1.0 / (V * g_prime**2)        # IRLS 가중치
        W = np.diag(w)
        beta = np.linalg.solve(X.T @ W @ X, X.T @ W @ z)
        eta_new = X @ beta
        mu_new = link_inv(eta_new)
        if np.max(np.abs(mu_new - mu)) < tol:
            mu, eta = mu_new, eta_new
            break
        mu, eta = mu_new, eta_new

    # 산포 추정
    pearson = np.sum((y - mu)**2 / V_fn(mu))
    phi = pearson / (n - p)
    cov_beta = phi * np.linalg.inv(X.T @ W @ X)
    se_beta = np.sqrt(np.diag(cov_beta))
    return {"beta": beta, "mu": mu, "phi": phi, "se": se_beta, "pearson": pearson}

# 로짓 링크 + V = mu(1-mu) (준-이항)
link_inv = lambda eta: 1 / (1 + np.exp(-eta))
link_deriv = lambda mu: 1 / (mu * (1 - mu))
V_binom = lambda mu: mu * (1 - mu)

# Wedderburn 제안 V = mu^2(1-mu)^2
V_wedd = lambda mu: (mu * (1 - mu))**2

11.2 Python — statsmodels 실무 코드

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.genmod.families import Binomial
from statsmodels.genmod.families.varfuncs import VarianceFunction

# 예: leaf_blotch DataFrame 가정
# 열: y (비율), site (범주), variety (범주)
X = pd.get_dummies(leaf_blotch[["site", "variety"]], drop_first=True).astype(float)
X = sm.add_constant(X)
y = leaf_blotch["y"].clip(0.001, 0.999)  # 경계 회피

# 1) V = mu(1-mu) — 준-이항
fit1 = sm.GLM(y, X, family=Binomial()).fit(scale="X2")
print(fit1.summary())
print(f"phi_hat = {fit1.scale:.4f}")  # ~ 0.089

# 2) V = mu^2(1-mu)^2 — Wedderburn 분산
class VarWedderburn(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 = VarWedderburn()
fit2 = sm.GLM(y, X, family=fam).fit(scale="X2")
print(fit2.summary())
print(f"phi_hat = {fit2.scale:.4f}")  # ~ 0.99

11.3 R — glm(family=quasi) 와 커스텀 family

library(stats)

# 1) 내장 quasi: V = mu(1-mu)
fit1 <- glm(y ~ site + variety,
            family = quasi(link = "logit", variance = "mu(1-mu)"),
            data = leaf_blotch)
summary(fit1)
# Dispersion parameter ~ 0.089

# 2) 커스텀 family: V = mu^2(1-mu)^2
quasi_wedd <- function() {
  linkfun  <- function(mu) log(mu/(1-mu))
  linkinv  <- function(eta) 1/(1 + exp(-eta))
  mu.eta   <- function(eta) {p <- 1/(1+exp(-eta)); p*(1-p)}
  variance <- function(mu) (mu*(1-mu))^2
  dev.resids <- function(y, mu, wt) {
    # Wedderburn 준-이탈도 근사 (y=0,1 에서 유한화 필요)
    eps <- 1e-8
    mu <- pmax(pmin(mu, 1-eps), eps)
    2 * wt * ((2*y-1)*log(mu/(1-mu)) - y/mu - (1-y)/(1-mu) + 2)
  }
  initialize <- expression({
    mustart <- pmax(pmin(y, 1 - 1e-3), 1e-3)
  })
  aic <- function(...) NA
  structure(list(family = "quasi_wedderburn", link = "logit",
                 linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta,
                 variance = variance, dev.resids = dev.resids,
                 aic = aic, initialize = initialize,
                 validmu = function(mu) all(mu > 0 & mu < 1)),
            class = "family")
}

fit2 <- glm(y ~ site + variety, family = quasi_wedd(), data = leaf_blotch)
summary(fit2)
# Dispersion parameter ~ 0.99

11.4 진단 코드

import matplotlib.pyplot as plt

# Pearson 잔차 vs 선형예측자
def diagnostic_plots(fit, V_fn, title):
    mu_hat = fit.fittedvalues
    eta_hat = fit.model.family.link(mu_hat)
    r_pearson = (fit.model.endog - mu_hat) / np.sqrt(V_fn(mu_hat))
    # 산포 보정
    r_std = r_pearson / np.sqrt(fit.scale)

    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    axes[0].scatter(eta_hat, r_std, alpha=0.6)
    axes[0].axhline(0, color='gray', linestyle='--')
    axes[0].set(xlabel=r'$\hat\eta$', ylabel='Std. Pearson residual', title=title)

    axes[1].scatter(np.log(V_fn(mu_hat)), r_std, alpha=0.6)
    axes[1].axhline(0, color='gray', linestyle='--')
    axes[1].set(xlabel=r'$\log V(\hat\mu)$', ylabel='Std. Pearson residual')
    plt.tight_layout()
    plt.show()

첫 번째 플롯(잔차 vs \(\widehat{\eta}\))에서 패턴이 사라지고, 두 번째(잔차 vs \(\log V\))에서 기울기가 없으면 분산함수가 적절하다.

12 관련 주제

선행 지식

• Quasi-likelihood Functions 개관 — 분포 없이 평균-분산만으로 추론
• GLM 이론 기초 — 지수족·정준연결·이탈도·IRLS
• GLM 적합 알고리즘 — IRLS 의 완전한 유도
• GLM 잔차의 세 유형 — Pearson·Anscombe·Deviance

같은 장(章) 후속

• Dependent Observations — 경로 독립성과 GEE (McCullagh §9.3) — placeholder
• Optimal Estimating Functions — Godambe 정보와 최적성 (McCullagh §9.4) — placeholder
• Extended Quasi-likelihood — 평균·산포 결합 모형 (McCullagh §9.5) — placeholder

관련 개념

• 이항 자료의 과산포 — Over-dispersion
• 다항 자료의 과산포 — Dirichlet-multinomial·\(\sigma^2\) 추정
• 감마 GLM — 분산함수 \(V(\mu) = \mu^2\)
• 로그선형 모형과 포아송 과산포

참고 문헌

• McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd Ed.), §9.2 (pp. 323–332).
• Wedderburn, R. W. M. (1974). Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method. Biometrika, 61(3), 439–447.
• Nelder, J. A., & Pregibon, D. (1987). An extended quasi-likelihood function. Biometrika, 74(2), 221–232.