Kwangmin Kim - Likelihood Functions for Polytomous GLMs

1 왜 “우도” 를 따로 다루는가

§5.2 에서 네 종류의 모형 구조를 봤고, §5.3 에서 다항분포의 성질을 정리했다. 이제 이 둘을 엮는 기계가 필요하다 — 우도함수(likelihood function). 다항 반응의 모든 추정·검정은 결국 다음 하나의 로그우도에서 나온다.

\[ \ell(\boldsymbol{\pi}; \mathbf{y}) = \sum_{i, j} y_{ij} \log \pi_{ij} \]

이 포스트는 McCullagh & Nelder (1989, §5.4) 의 세 소절 — 로그우도(§5.4.1), 모수 추정 (§5.4.2), 이탈도(§5.4.3) — 을 따라가며 “왜 이 식이 이 형태인가”를 직관적으로 풀어낸다. §5.5 의 과산포 처리와 §5.6 의 치즈 실험 예제까지 함께 다룬다.

2 §5.4.1 로그우도 — 제약 하의 미분

2.1 설정

\(n\) 개의 독립 다항 관측 \(\mathbf{y}_1, \ldots, \mathbf{y}_n\) 이 있고, 각 \(\mathbf{y}_i = (y_{i1}, \ldots, y_{ik})\) 는 \(\sum_j y_{ij} = m_i\) 로 합이 고정. 확률벡터 \(\boldsymbol{\pi}_i\) 는 \(\sum_j \pi_{ij} = 1\) 제약.

일차 단계에서는 모형을 특정하지 않는다. \(\pi_{ij}\) 를 자유 모수처럼 다뤄 로그우도와 그 미분을 일반적으로 구한 뒤, 나중에 비례 오즈나 로그선형 같은 구체적 모형 (예: \(\text{logit}\,\gamma_{ij} = \mathbf{x}_i^\top \boldsymbol{\beta}\)) 을 대입한다. 이 방식은 모형과 우도를 독립적으로 설계할 수 있게 해 준다.

2.2 로그우도 — 식 (5.14)

한 관측의 기여는 (다항 계수 \(\binom{m_i}{\mathbf{y}_i}\) 는 모수 무관 상수이므로 제거)

\[ \ell(\boldsymbol{\pi}_i; \mathbf{y}_i) = \sum_j y_{ij} \log \pi_{ij}. \]

전체 로그우도는 독립 가정으로 단순 합:

\[ \ell(\boldsymbol{\pi}; \mathbf{y}) = \sum_{i j} y_{ij} \log \pi_{ij}. \tag{5.14} \]

직관 — 이 식의 형태: “관측된 범주의 로그 확률을 관측 횟수만큼 더한다” 는 단순한 구조. \(y_{ij}\) 가 클수록 그 범주의 확률을 높여 주는 방향으로 \(\pi_{ij}\) 를 밀어붙이는 힘이 커진다. 이것이 MLE 가 경험 비율 \(\tilde{\pi}_{ij} = y_{ij}/m_i\) 로 수렴하는 이유이다 (나중에 §5.4.3 에서 확인).

2.3 제약 하의 미분 — 스코어 식 (5.15)

\(\sum_j \pi_{ij} = 1\) 제약 하에 \(\pi_{ij}\) 로 미분하면

\[ \frac{\partial \ell}{\partial \pi_{ij}} = \frac{y_{ij} - m_i \pi_{ij}}{\pi_{ij}}. \]

어떻게 나오는가: \(\ell_i = \sum_j y_{ij} \log \pi_{ij}\) 를 직접 미분하면 \(y_{ij}/\pi_{ij}\) 가 나오지만, 제약을 다루기 위해 한 범주(예: \(\pi_{ik} = 1 - \sum_{j < k} \pi_{ij}\))를 나머지로 표현하고 연쇄법칙을 적용하면 위 형태가 된다. \(y_{ij}/\pi_{ij} - m_i\) 를 공통분모로 묶은 것으로도 해석할 수 있다.

행렬 형태:

\[ \frac{\partial \ell}{\partial \boldsymbol{\pi}_i} = m_i \boldsymbol{\Sigma}_i^{-} (\mathbf{y}_i - m_i \boldsymbol{\pi}_i) = m_i \boldsymbol{\Sigma}_i^{-} (\mathbf{y}_i - \boldsymbol{\mu}_i). \tag{5.15} \]

여기서 \(\boldsymbol{\Sigma}_i = m_i\{\mathrm{diag}(\boldsymbol{\pi}_i) - \boldsymbol{\pi}_i \boldsymbol{\pi}_i^\top\}\) 의 일반화역행렬로 보통 대각 선택 \(\boldsymbol{\Sigma}_i^- = \mathrm{diag}\{1/(m_i \pi_{ij})\}\) 를 쓴다.

일반화역행렬 선택이 결과에 영향을 안 주는 이유: 잔차 \(\mathbf{y}_i - \boldsymbol{\mu}_i\) 가 \(\mathbf{1}^\top(\mathbf{y}_i - \boldsymbol{\mu}_i) = 0\) 을 만족하므로, 상수 벡터 방향으로의 차이(§5.3.3에서 본 \(\boldsymbol{\Sigma}^- - c\mathbf{1}\mathbf{1}^\top\) 자유도)가 이차형식에서 사라진다. 이것은 일반 GLM 에서도 반복되는 원리이다.

왜 스코어가 \(\mathbf{y} - \boldsymbol{\mu}\) 의 가중합인가: 모든 지수족 GLM 에서 성립하는 보편적 사실이다. 다항도 지수족의 특수 사례이므로 같은 구조가 나타난다. “관측-예측 잔차에 분산으로 가중” 한 것이 스코어이다.

2.4 누적 표현 — 식 (5.16)

누적 확률 \(\gamma_{ij} = \sum_{\ell \le j} \pi_{i\ell}\) 로 다시 쓰면

\[ \frac{\partial \ell}{\partial \gamma_{ij}} = \frac{y_{ij} - m_i \pi_{ij}}{\pi_{ij}} - \frac{y_{i,j-1} - m_i \pi_{i,j-1}}{\pi_{i,j-1}}, \quad 1 < j < k, \]

\[ \frac{\partial \ell}{\partial \boldsymbol{\gamma}_i} = m_i \boldsymbol{\Gamma}_i^{-}(\mathbf{z}_i - m_i \boldsymbol{\gamma}_i). \tag{5.16} \]

체인 룰: \(\partial \ell / \partial \gamma_{ij} = \partial \ell / \partial \pi_{ij} - \partial \ell / \partial \pi_{i,j-1}\).

왜 두 표현을 모두 아는 것이 중요한가: 모형이 \(\boldsymbol{\pi}\) 로 직접 표현되면(명목형 로그선형) 식 (5.15) 를, \(\boldsymbol{\gamma}\) 로 표현되면(순서형 비례 오즈) 식 (5.16) 을 쓴다. 수학적으로 동치이지만 계산 복잡도가 다르다 — 모형과 어울리는 파라미터화를 택하면 공식이 훨씬 단순해진다.

두 표현의 관계 요약

표현	파라미터	잔차	g-inverse	어울리는 모형
셀 확률	\(\pi_{ij}\)	\(y_{ij} - m_i \pi_{ij}\)	\(\mathrm{diag}\{1/(m_i \pi_{ij})\}\)	명목(기준범주), 로그선형
누적 확률	\(\gamma_{ij}\)	\(z_{ij} - m_i \gamma_{ij}\)	삼중대각 Jacobi	순서(비례오즈), 비례위험

두 표현이 주는 로그우도 값과 MLE 는 동일하다. 다만 계산 식이 다르다.

3 §5.4.2 모수 추정 — 두 가지 모형군

모수 \(\boldsymbol{\beta}\) 에 대한 스코어 방정식은 원리상 단순하다.

\[ \frac{\partial \ell}{\partial \beta_r} = \sum_{ij} \frac{\partial \ell}{\partial \pi_{ij}} \frac{\partial \pi_{ij}}{\partial \beta_r} \]

또는 동등하게

\[ \frac{\partial \ell}{\partial \beta_r} = \sum_{ij} \frac{\partial \ell}{\partial \gamma_{ij}} \frac{\partial \gamma_{ij}}{\partial \beta_r}. \]

구체적 형태는 모형 이 결정한다. 두 주요 사례를 본다.

3.1 사례 1 — 비례 오즈 모형의 스코어

모형 (5.1): \(\text{logit}\,\gamma_{ij} = \theta_j - \boldsymbol{\beta}^\top \mathbf{x}_i\).

\(\boldsymbol{\beta}^* = (\theta_1, \ldots, \theta_{k-1}, \beta_1, \ldots, \beta_p)\) 으로 모두 모은 \(p^* = p + k - 1\) 차원 모수이다. 디자인 행렬 \(\mathbf{X}^*\) 는 \(nk \times p^*\) 로, \((i, j)\) 행이 \((0, \ldots, \underbrace{1}_{\text{위치 } j}, \ldots, 0, \mathbf{x}_i^\top)\) 이다. 각 관측 \(i\) 에 대한 \(k-1\) 행 블록이

\[ [\mathbf{I}_{k-1} : \mathbf{1}\,\mathbf{x}_i^\top] \]

형태를 갖는다. 스코어는

\[ \frac{\partial \ell}{\partial \beta_r^*} = \sum_{ij} x_{ijr}^* \gamma_{ij}(1 - \gamma_{ij}) \frac{\partial \ell}{\partial \gamma_{ij}}. \]

\(\gamma_{ij}(1 - \gamma_{ij})\) 의 역할: 로짓 역함수 \(\gamma = \text{expit}(\eta)\) 의 도함수가 \(\gamma(1 - \gamma)\) 이므로 체인 룰의 기여가 나타난 것. 이것은 이항 로지스틱에서의 스코어 가중치와 정확히 같은 형태 — 비례 오즈가 “여러 이항을 묶은 것”으로 환원되는 이유이다.

직관 — 경계 근처의 정보 집중: \(\gamma_{ij}(1 - \gamma_{ij})\) 는 \(\gamma_{ij} = 1/2\) 에서 최대, \(0\) 이나 \(1\) 근처에서 \(0\). 즉 경계 확률이 중간값일 때 그 경계가 모수 추정에 가장 많이 기여한다. 극단 범주(거의 비어있는 범주)는 정보를 거의 주지 않는다 — 치즈 실험에서 극단 셀의 적합도가 낮게 보이는 이유이다.

3.2 사례 2 — 로그선형 모형의 스코어 (적률법 동치)

로그선형 모형 \(\log \pi_{ij} = \sum_r x_{ijr}^* \beta_r^*\) (식 5.5~5.9 포함) 의 경우 스코어 방정식은 특히 단순하다.

\[ \sum_{ij} x_{ijr}^* (y_{ij} - \hat{\mu}_{ij}) = 0, \quad r = 1, \ldots, p^*. \]

이것이 최대우도 = 적률법(method of moments) 의 등식이다. “특정 선형조합 \(\sum x_{ijr}^* y_{ij}\) 를 그 기대값과 일치시켜라” 는 지시.

직관 — 어떤 통계량이 일치되는가:

모형	적률법이 일치시키는 통계량
식 (5.8) (점수 × 공변량)	행 합, 열 합, 교호작용 \(\sum_{ij} x_{ir} s_j y_{ij}\)
식 (5.9) (기준범주 로짓)	각 범주별 공변량 합
식 (5.6) (공변량 효과 없음)	행 합, 열 합만 일치

이 구조 덕분에 로그선형 다항 모형은 충분통계량(sufficient statistic) 이 관측된 선형조합 이다. MLE 계산은 이들 선형조합이 일치할 때까지 \(\hat{\boldsymbol{\pi}}\) 를 조정하는 것이다. 반복 비례 적합(iterative proportional fitting) 이 이 아이디어의 고전 알고리즘이다.

3.3 실무 알고리즘 — IRLS 의 다변량 확장

두 사례 모두 반복 재가중 최소제곱(IRLS) 으로 풀린다. 이항 GLM 의 \(W_i = \mu_i(1 - \mu_i)/m_i\) 단일 스칼라가, 다항에서는 \(k \times k\) 행렬 \(\boldsymbol{W}_i \propto \mathrm{diag}(\boldsymbol{\mu}_i) - \boldsymbol{\mu}_i \boldsymbol{\mu}_i^\top / m_i\) 로 바뀔 뿐 구조는 동일하다. 반복식:

\[ \boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + (\mathbf{X}^{*\top} \boldsymbol{W}^{(t)} \mathbf{X}^*)^{-1} \mathbf{X}^{*\top} (\mathbf{y} - \boldsymbol{\mu}^{(t)}). \]

비례 오즈와 같은 순서형 모형에서는 디자인 행렬의 블록 구조 때문에 계산이 다소 복잡하지만 근본은 같은 뉴턴-라프슨 반복이다. 로그선형의 경우 정준 연결(canonical link)이므로 관측 정보 = 기대 정보가 되어 피셔 스코어 = 뉴턴-라프슨 이다.

4 §5.4.3 이탈도 (Deviance) — 적합도 측정

4.1 유도 — 포화 모형과의 차이

다항 로그우도가 최대가 되는 점은 포화 모형(saturated model) 에서이다. 제약 \(\sum_j \pi_{ij} = 1\) 하에 \(\sum_j y_{ij} \log \pi_{ij}\) 를 최대화하면

\[ \tilde{\pi}_{ij} = \frac{y_{ij}}{m_i} \]

(경험 비율)가 나온다. 이것이 도달 가능한 최대 로그우도.

이탈도는 포화 로그우도와 현재 모형의 로그우도 차이의 2배:

\[ D(\mathbf{y}; \hat{\boldsymbol{\pi}}) = 2\ell(\tilde{\boldsymbol{\pi}}; \mathbf{y}) - 2\ell(\hat{\boldsymbol{\pi}}; \mathbf{y}) = 2 \sum_{ij} y_{ij} \log \tilde{\pi}_{ij} - 2 \sum_{ij} y_{ij} \log \hat{\pi}_{ij}. \]

공통 정리:

\[ D(\mathbf{y}; \hat{\boldsymbol{\pi}}) = 2 \sum_{ij} y_{ij} \log\!\left(\frac{y_{ij}}{\hat{\mu}_{ij}}\right). \]

4.2 수식의 직관적 읽기

이 식은 Ch.6 에서 다시 등장할 Kullback–Leibler 발산(KL divergence) 의 \(2m\) 배이다.

\[ D = 2 \sum_{ij} y_{ij} \log\!\left(\frac{y_{ij}/m_i}{\hat{\mu}_{ij}/m_i}\right) = 2 \sum_i m_i \, \mathrm{KL}(\tilde{\boldsymbol{\pi}}_i \,\|\, \hat{\boldsymbol{\pi}}_i). \]

경험 분포 \(\tilde{\boldsymbol{\pi}}_i\) 와 적합 분포 \(\hat{\boldsymbol{\pi}}_i\) 의 정보량 차이
완벽 적합이면 \(\tilde{\boldsymbol{\pi}}_i = \hat{\boldsymbol{\pi}}_i\) 이므로 \(D = 0\)
두 분포가 멀수록 \(D\) 가 크다

직관 — 왜 포화 모형이 기준인가: 포화 모형은 데이터를 완벽히 기억하는 모형이다. “데이터가 말하는 만큼 다 믿으면 로그우도가 얼마가 되는가?” 가 기준선이고, 우리가 제안한 구조적 모형이 그 기준선에서 얼마나 떨어졌는지가 이탈도이다. 관심은 이 거리가 “통계적 우연의 범위 내”인지 여부.

4.3 분포 근사

조건(기대도수 \(\hat{\mu}_{ij}\) 가 충분히 크고 과산포 없음) 하에 \(D\) 는 대략 \(\chi^2_\nu\) (자유도 \(\nu = n(k-1) - p\)) 분포를 따른다.

\(n(k-1) - p\) 의 의미: 전체 관측 자유도 \(n(k-1)\) (각 \(\mathbf{y}_i\) 가 \(k-1\) 자유도) 에서 적합된 모수 수 \(p\) 를 뺀 것.

한계 — 언제 \(\chi^2\) 근사가 깨지는가: 기대도수가 작은 셀(예: \(\hat{\mu}_{ij} < 5\)) 이 많거나, 다항 관측이 “희소(sparse)” 한 경우 (거의 \(m_i = 1\)). 이 경우 Pearson \(X^2\) 과 \(D\) 모두 \(\chi^2\) 근사가 나쁘며, 이탈도의 차이(두 중첩 모형 간 LR 통계량) 는 여전히 잘 작동하지만 절대값을 적합도 검정으로 쓰는 것은 주의해야 한다.

4.4 Pearson \(X^2\) 와의 비교

\(X^2 = \sum_{ij}(y_{ij} - \hat{\mu}_{ij})^2 / \hat{\mu}_{ij}\) 이 되는데, \(y_{ij} \approx \hat{\mu}_{ij}\) 근처에서 Taylor 전개하면 \(D \approx X^2\). 둘 다 같은 \(\chi^2\) 극한을 갖는다. 차이: \(D\) 는 로그우도 기반이어서 모형 비교(LRT)에 자연스럽고, \(X^2\) 는 분산 기반 잔차의 합이어서 개별 잔차 분석에 유용하다.

5 §5.5 과산포 (Over-dispersion)

5.1 문제 설정

실제 데이터의 분산이 다항 가정보다 크다면(군집 샘플링, 패널리스트 상관 등), 다항 공분산 \(\boldsymbol{\Sigma}\) 를 산포 인자 \(\sigma^2\) 로 부풀린다.

\[ \mathrm{E}(\mathbf{Y}) = m\boldsymbol{\pi}, \quad \mathrm{Cov}(\mathbf{Y}) = \sigma^2 \boldsymbol{\Sigma}. \tag{5.17} \]

\(\sigma^2 > 1\) 이면 과산포(over-dispersion), \(\sigma^2 < 1\) 이면 과소산포(under-dispersion). 군집 샘플링에서 cluster 내 양의 상관이 있으면 \(\sigma^2 > 1\) 이 전형적이다.

5.2 추정 — 점추정은 그대로, 표준오차만 보정

\(\hat{\boldsymbol{\beta}}\) 자체는 영향을 받지 않는다 — 스코어 방정식이 일차 모멘트만 사용
\(\mathrm{Cov}(\hat{\boldsymbol{\beta}})\) 에 \(\sigma^2\) 을 곱해 부풀린다
\(\sigma^2\) 은 Pearson \(X^2\) 로부터 추정:

\[ \hat{\sigma}^2 = \frac{X^2}{n(k-1) - p} = \frac{X^2}{\text{residual d.f.}}. \tag{5.18} \]

직관: “모형의 점추정은 여전히 데이터의 중심을 맞추지만, 그 중심의 불확실성을 과소평가하고 있었다. 잔차의 총 이차합이 기대치보다 크면 그 배수로 표준오차를 키운다.” 이 접근법은 quasi-likelihood 의 철학(Ch.9) 과 정확히 같은 논리이다.

왜 이탈도 대신 \(X^2\) 를 쓰는가: Pearson \(X^2\) 은 잔차의 제곱합 구조를 가져 \(\sigma^2\) 추정에 편향이 적고, 희소 데이터에서도 \(\hat{\beta}\) 와 거의 독립이다. 이탈도는 로그 비율 구조라 \(\sigma^2\) 추정량으로는 덜 안정적.

6 §5.6.1 예제 — 치즈 맛 실험

6.1 데이터

Newell 박사의 실험: 네 가지 첨가물 A, B, C, D 를 52명 패널리스트에 대해 9점 척도 (1 = strong dislike ~ 9 = excellent taste) 로 평가. 각 첨가물 \(n_i = 52\), 총 \(N = 208\).

관측 분포(Table 5.1 요약):

첨가물	I	II	III	IV	V	VI	VII	VIII	IX	합
A	0	0	1	7	8	8	19	8	1	52
B	6	9	12	11	7	6	1	0	0	52
C	1	1	6	8	23	7	5	1	0	52
D	0	0	0	1	3	7	14	16	11	52

시각적으로도 질적 순위 (D, A, C, B) 가 보인다. 모형은 이 순위를 정량화하는 도구.

6.2 비례 오즈 모형 적합

모형 (5.1):

\[ \text{logit}\,\gamma_{ij} = \theta_j - \beta_i, \quad j = 1, \ldots, 8, \; i \in \{A, B, C, D\} \]

식별성으로 \(\hat{\beta}_D = 0\). MLE:

첨가물	\(\hat{\beta}\)	SE
A	\(-1.613\)	0.378
B	\(-4.965\)	0.474
C	\(-3.323\)	0.425
D	\(0.0\) (기준)	—

해석:

\(\beta\) 가 양수 → D 기준 대비 높은 점수로 이동
모든 \(\hat{\beta}_i < 0\) → D 가 제일 좋음
순위 \(D > A > C > B\) (대수적으로도 확인)
\(\hat{\beta}_B - \hat{\beta}_A = -3.35\) → B 가 A 보다 9점 척도의 어느 경계를 잘라도 동일하게 약 \(e^{3.35} \approx 28.5\) 배 더 낮은 점수 쪽 오즈. 이것이 비례 오즈 해석의 힘.

6.3 이탈도와 적합도

모형	이탈도	자유도
귀무(\(\beta = 0\))	168.8	24
비례 오즈	20.31	21

\(\Delta D = 148.5\) on 3 df — 첨가물 효과가 극도로 유의. \(\chi^2_3\) 분포의 0.999 분위수보다 훨씬 크다.

잔차 분석의 주의: 극단 셀의 \(y_{ij}\) 가 작아 \(\chi^2\) 근사가 나쁘다. 표준화 잔차 \((y_{ij} - \hat{y}_{ij}) / [m_i \hat{\pi}_{ij}(1 - \hat{\pi}_{ij})]^{1/2}\) 기준 2.0 초과가 두 개 (셀 (1,4) = 2.23, (2,6) = 2.30). 누적 표현 \(z_{ij}\) 에 기반한 잔차를 쓰면 이 괴리가 사라지며, McCullagh 는 이것이 개념적으로 더 적절하다고 주장한다 — \(k-1\) 개 자유도에 정확히 대응하기 때문.

6.4 비례 오즈 가정 점검 — 식 (5.4) 확장 적합

\(\text{logit}\,\gamma_{ij} = (\theta_j - \beta_i)/\exp(\tau_i)\) 로 확장하면 \(\Delta D = 3.3\) on 3 df → 유의하지 않음. 분산 구조가 첨가물 간 다르다는 증거 없음 → 비례 오즈 가정 유지.

6.5 과산포 확인

\(\hat{\sigma}^2 = 20.9 / 21 \approx 1.00\). 즉 과산포 없음. \(\sigma^2 = 1\) 가정이 정당화된다.

6.6 실무적 주의 — 관측 독립성 위반

같은 52명 패널리스트가 네 첨가물을 모두 평가 → 반응 간 양의 상관 \(\rho\) 존재 가능성. 이 경우 대조(contrast) 의 분산은 독립 가정 대비 \((1 - \rho^2)\) 배로 축소되어, 본 계산의 SE 는 상한으로 해석된다. 결론은 보수적(conservative) — 효과 크기는 더 확신할 만하다.

6.7 범주 병합 실험

\(\{1, 2, 3, 4\}\) 병합 + \(\{7, 8, 9\}\) 병합 → 4범주로 축소. 새 추정:

\(\hat{\boldsymbol{\beta}} = (-1.34, -4.57, -3.07, 0)\)

평균 0.7 SE 정도 감소, 분산 약 19% 증가, 상관은 거의 불변.

메시지: 범주 통합 불변성은 완전 불변 을 의미하지 않는다. 모수가 추정하는 대상은 같지만 추정 효율 은 범주 수에 따라 달라진다. 범주가 많을수록 정보가 많지만, 극단 셀이 비어 \(\chi^2\) 근사가 나빠지는 트레이드오프가 있다.

7 코드 예시 — 치즈 데이터 재현

7.1 Step 1: 순수 Python — 비례 오즈 로그우도와 뉴턴 최적화

import numpy as np
from scipy.optimize import minimize
from scipy.special import expit

# 치즈 실험 데이터 — (첨가물, 범주) 도수
counts = np.array([
    [0, 0, 1,  7,  8,  8, 19,  8,  1],  # A
    [6, 9,12, 11,  7,  6,  1,  0,  0],  # B
    [1, 1, 6,  8, 23,  7,  5,  1,  0],  # C
    [0, 0, 0,  1,  3,  7, 14, 16, 11],  # D
], dtype=float)
n_treat, k = counts.shape       # 4, 9
m = counts.sum(axis=1)          # 각 첨가물 52명


def neg_loglik(params):
    # params: theta raw (k-1) + beta (n_treat-1)  — D 를 기준 0
    raw = params[:k-1]
    # 단조증가 강제: theta_1, theta_2 = theta_1 + exp(.), ...
    theta = np.concatenate([[raw[0]], raw[0] + np.cumsum(np.exp(raw[1:]))])
    beta = np.concatenate([params[k-1:], [0.0]])   # D = 0 기준

    ll = 0.0
    for i in range(n_treat):
        gamma = expit(theta - beta[i])                # k-1
        gamma = np.concatenate([[0.0], gamma, [1.0]])  # 양끝 경계
        pi = np.diff(gamma)                            # 범주 확률 k
        ll += np.sum(counts[i] * np.log(np.maximum(pi, 1e-12)))
    return -ll


init = np.concatenate([np.linspace(-3, 3, k-1)[0:1],
                       np.log(np.ones(k-2) * 0.5),
                       np.zeros(n_treat - 1)])
res = minimize(neg_loglik, init, method="BFGS")

raw = res.x[:k-1]
theta_hat = np.concatenate([[raw[0]], raw[0] + np.cumsum(np.exp(raw[1:]))])
beta_hat = np.concatenate([res.x[k-1:], [0.0]])
names = ["A", "B", "C", "D"]
for name, b in zip(names, beta_hat):
    print(f"beta_{name} = {b:+.3f}")

# 이탈도 계산
ll_fit = -res.fun
tilde_pi = counts / m[:, None]
mask = counts > 0
ll_sat = np.sum(counts[mask] * np.log(tilde_pi[mask]))
D = 2 * (ll_sat - ll_fit)
df = n_treat * (k - 1) - len(res.x)
print(f"Deviance = {D:.2f}  on  {df} df")

실행하면 교재 추정치에 근접한 \(\hat{\beta} \approx (-1.61, -4.96, -3.32, 0)\) 와 이탈도 약 20 에 도달한다 — 작은 수치 차이는 수치 최적화의 수렴 허용 오차.

7.2 Step 2: statsmodels / R 로 실무 적합

import pandas as pd
import numpy as np
from statsmodels.miscmodels.ordinal_model import OrderedModel

# long-format 재구성
rows = []
names = ["A", "B", "C", "D"]
for i, name in enumerate(names):
    for j in range(k):
        rows.extend([{"cheese": name, "rating": j + 1}] * int(counts[i, j]))
df = pd.DataFrame(rows)

X = pd.get_dummies(df["cheese"], drop_first=True).astype(float)  # B, C, D 대비
mod = OrderedModel(df["rating"], X, distr="logit")
res = mod.fit(method="bfgs", disp=False)
print(res.summary())

library(MASS)
data <- # ... expand from Table 5.1
fit <- polr(factor(rating) ~ cheese, data = data, method = "logistic")
summary(fit)
anova(polr(factor(rating) ~ 1), fit)   # LRT for additive effect

8 자주 걸리는 함정

함정	증상	처방
제약 무시 미분	스코어가 \(\sum_j = 0\) 안 됨	Lagrange 또는 자유 모수 \(k-1\) 개
로그선형에 정준 아닌 링크 오해	Fisher ≠ 관측 정보 혼동	로그선형은 정준 링크 → 둘 동일
극단 셀 잔차에 과도한 의미 부여	\(\chi^2\) 근사 나쁜 영역	누적 잔차 \(z_{ij}\) 기반 분석 선호
이탈도를 희소 데이터의 적합도로 사용	\(\chi^2\) 근사 붕괴	이탈도 차이로 중첩 모형 비교만
\(\sigma^2 = 1\) 맹신	군집·반복측정 상황 오염	\(X^2/df\) 로 \(\hat{\sigma}^2\) 추정
패널리스트 반복측정 무시	SE 과대 또는 과소평가	혼합효과 모형 또는 GEE
범주 병합이 결과 바꾼다고 불안	효율만 바뀜, 모수는 동일	의도적 민감도 분석으로 보고

9 관련 주제

선행 지식

후속 주제 (placeholder)

관련 개념

로그선형 모형 (Ch.6) — 충분통계량 = 주변합
이항 자료의 우도함수 — \(k=2\) 특수 사례
KL Divergence 와 이탈도 — 정보기하적 해석

10 참고문헌

McCullagh, P. & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.), §5.4–§5.6. Chapman & Hall.
McCullagh, P. (1980). Regression models for ordinal data. JRSS B, 42(2), 109–142.
Wedderburn, R. W. M. (1974). Quasi-likelihood functions, generalized linear models and the Gauss–Newton method. Biometrika, 61, 439–447.
Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley.