Some Applications Involving Polytomous Data

매칭 쌍 명목 반응 · 순서형 반응 · 치즈 맛 실험 (McCullagh §7.5)

조건부 우도 이론을 다범주 반응(polytomous data)에 적용한다. 명목형 매칭 쌍에서 quasi-symmetry/Bradley-Terry 모형을 유도하고, 순서형 반응에서 비례 오즈 모형의 조건부 추정을 전개한 뒤, 치즈 맛 실험으로 반복 계산 과정을 시연한다.

Statistics
GLM
저자

Kwangmin Kim

공개

2026년 04월 18일

1 이 장의 위치

§7.4에서 이항 자료(binary data)에 조건부 우도를 적용했다면, §7.5는 반응이 3개 이상의 범주를 가지는 다범주 자료(polytomous data)로 무대를 확장한다.

핵심 구조는 동일하다: 매칭 쌍(matched pairs) 설계에서 쌍별 기저 확률 \(\lambda_i\) 가 장해 모수이고, 충분통계량으로 조건부를 취해 관심 모수 \(\Delta\) 만 남긴다. 차이는 반응 척도에 따라 조건부 분포와 추정 전략이 달라진다는 점이다.

절 반응 척도 조건부 모형 핵심 결과
§7.5.1 명목(nominal) \(k(k-1)/2\) 개 독립 이항 Quasi-symmetry = Bradley-Terry
§7.5.2 순서(ordinal) \(k-1\) 개 비중심 초기하 가중 추정 방정식(quasi-likelihood)
§7.5.3 순서(예제) 치즈 맛 실험 조건부 vs 무조건부 MLE 비교

2 매칭 쌍: 명목형 반응 (§7.5.1)

2.1 설정

피험자를 쌍(pair)으로 매칭하고, 각 피험자에게서 \(k\) 개 범주 중 하나의 반응을 관측한다. 쌍 \(i\) 에서 대조군의 로그 반응 확률은

\[ \lambda_i = (\lambda_{i1}, \ldots, \lambda_{ik}), \]

처치군의 로그 반응 확률은

\[ \lambda_i + \Delta = (\lambda_{i1} + \Delta_1, \ldots, \lambda_{ik} + \Delta_k). \]

직관: \(\lambda_i\) 는 쌍마다 완전히 다를 수 있는 “기저 프로필”이고, \(\Delta = (\Delta_1, \ldots, \Delta_k)\) 는 처치가 각 범주의 로그 확률을 일정하게 밀어 올리거나 내리는 공통 처치 효과이다. 이항 자료(§7.4)에서 \(\lambda_i\) 가 스칼라였다면, 여기서는 \(k\) 차원 벡터로 확장된 것이다.

범주 \(j\) 의 반응 확률은 소프트맥스 형태이다:

  • 대조군: \(\Pr(R = j) = \exp(\lambda_{ij}) / \sum_r \exp(\lambda_{ir})\)
  • 처치군: \(\Pr(R = j) = \exp(\lambda_{ij} + \Delta_j) / \sum_r \exp(\lambda_{ir} + \Delta_r)\)

2.2 충분통계량과 조건부 분포

한 쌍의 관측 반응을 지시벡터(indicator vector) \(Z_1, Z_2\) 로 표현한다:

\[ Z_j = \begin{cases} 1 & R = j \\ 0 & \text{otherwise} \end{cases} \]

벡터합 \(Z_\bullet = Z_1 + Z_2\) 가 \(\lambda\) 의 충분통계량이다.

두 가지 경우:

  • \(Z_\bullet = (0, \ldots, 2, \ldots, 0)\): 두 피험자가 같은 범주에 응답. \(Z_\bullet\) 이 \((R_1, R_2)\) 를 완전히 결정하므로 조건부 분포가 퇴화(degenerate)한다. 이 쌍은 조건부 우도에 기여하지 않는다 — 이항 매칭 쌍에서 일치하는 쌍을 버리는 것과 같은 원리이다.

  • \(Z_\bullet = (0, \ldots, 1_i, \ldots, 1_j, \ldots, 0)\): 두 피험자가 범주 \(i, j\) (\(i \neq j\))에 각각 응답. 가능한 배치는 \((R_1, R_2) = (i, j)\) 또는 \((j, i)\) 뿐이다.

\(i \neq j\) 인 경우의 조건부 확률은

\[ \Pr(R_1 = i \mid Z_\bullet) = \frac{e^{\lambda_i} e^{\lambda_j + \Delta_j}}{e^{\lambda_i} e^{\lambda_j + \Delta_j} + e^{\lambda_j} e^{\lambda_i + \Delta_i}} = \frac{e^{\Delta_j}}{e^{\Delta_i} + e^{\Delta_j}}. \]

직관 — 왜 lambda가 사라지는가

분자와 분모 모두 \(e^{\lambda_i + \lambda_j}\) 를 인수로 가지므로 약분된다. 이것은 이항 매칭 쌍에서 \(Y_1 | Y_\bullet\) 이 \(\lambda\) 에 무관한 것과 정확히 같은 메커니즘이다. 다범주로 확장했을 뿐, 충분통계량 조건부 소거의 원리는 동일하다.

2.3 Quasi-Symmetry 모형

모든 비일치 쌍 \((i, j)\) 에 대해 위 조건부 확률을 모으면, \(Y_{ij}\) (반응이 \((i, j)\) 인 쌍의 수)에 대해

\[ Y_{ij} \sim B(m_{ij},\, \pi_{ij}), \qquad i < j, \]

\[ \operatorname{logit}(\pi_{ij}) = \Delta_j - \Delta_i. \]

여기서 \(m_{ij} = Y_{ij} + Y_{ji}\) 는 범주 \(i, j\) 에 응답한 쌍의 대칭 합계이다.

핵심 관찰 — 차이만 남는다. 조건부 분포의 로짓에 \(\Delta_i, \Delta_j\) 의 절대값이 아니라 차이 \(\Delta_j - \Delta_i\) 만 등장한다. 각 쌍의 기저 선호 수준 \(\lambda_{ij}\) (장해 모수) 가 조건화로 소거되고, 범주 간 상대적 선호 척도 만 남은 것이다. 그 결과 \(\{\Delta_j\}\) 에 공통 상수를 더해도 모든 \(\pi_{ij}\) 가 불변 — 이 구조적 비식별성 때문에 보통 \(\Delta_1 = 0\) 같은 정규화 제약을 덧붙인다.

직관: 조건부 우도는 \(k(k-1)/2\) 개의 독립 이항 인수의 곱이다. 각 인수는 “범주 \(i\) 와 \(j\) 에 응답한 쌍들 중에서, 대조군이 \(i\) 이고 처치군이 \(j\) 일 확률”을 로짓 모형으로 표현한다.

이 모형을 quasi-symmetry 모형이라 부르며, Caussinus (1965)가 처음 제안했다. 인구 이동(migration) 연구에서는 같은 모형이 중력 모형(gravity model) 으로 불린다.

2.4 Bradley-Terry 모형과의 동일성

모형 \(\operatorname{logit}(\pi_{ij}) = \Delta_j - \Delta_i\) 는 Bradley-Terry (1952) 모형과 형식적으로 동일하다.

  • \(k\) 명의 선수가 쌍별 대결을 한다
  • \(\pi_{ij}\) = 선수 \(i\) 가 선수 \(j\) 를 이길 확률
  • \(\Delta_j\) = 선수 \(j\) 의 “능력” (strength parameter)

직관: 매칭 쌍 실험에서 “처치군이 범주 \(j\) 에 응답”하는 것과 “선수 \(j\) 가 이기는 것”은 수학적으로 동일한 구조이다. 체스 Elo 레이팅, 스포츠 순위 시스템, 검색 엔진의 쌍별 비교 등 광범위한 응용이 이 모형에 기반한다.

2.5 모형 행렬의 특이한 구조

\(\Delta_j - \Delta_i\) 공식에 대응하는 모형 행렬 \(X\) 는 절편(intercept)을 포함하지 않으며, 상수 벡터가 \(X\) 의 열공간에 속하지 않는다.

\(k = 3\) 인 경우를 명시적으로 쓰면

\[ \begin{pmatrix} \operatorname{logit} \pi_{12} \\ \operatorname{logit} \pi_{13} \\ \operatorname{logit} \pi_{23} \end{pmatrix} = \begin{pmatrix} 1 & -1 & 0 \\ 1 & 0 & -1 \\ 0 & 1 & -1 \end{pmatrix} \begin{pmatrix} \Delta_1 \\ \Delta_2 \\ \Delta_3 \end{pmatrix}. \]

이 \(3 \times 3\) 행렬은 랭크 2이고, 세 열의 합이 \(\mathbf{0}\) 이다.

직관: 각 행은 “차이”를 나타내므로 수준 상수(level constant)가 사라진다. 이는 ANOVA에서 처치 효과의 합이 0이 되도록 제약하는 것과 같은 원리이다. 따라서 \(\Delta\) 에는 항상 하나의 제약이 필요하다: \(\Delta_1 = 0\) (기준 범주) 또는 \(\sum \Delta_j = 0\) (합 제약).

2.6 모형 적합도 검정

quasi-symmetry 모형의 잔차 이탈도 또는 Pearson \(X^2\) 통계량의 자유도는

\[ \frac{(k-1)(k-2)}{2}. \]

\(k = 3\) 이면 자유도 1, \(k = 4\) 이면 자유도 3이다. 자유도가 비교적 작으므로 모형 검정의 검정력이 높지 않을 수 있다.

3 매칭 쌍: 순서형 반응 (§7.5.2)

§7.5.1 의 명목형 반응은 범주 간 순서가 없어 각 범주의 선호 강도를 독립된 모수 \(\delta_j\) 로 따로 추정했다. 그러나 리커트 척도나 등급 평가처럼 범주가 자연 순서 (매우 싫음 < 싫음 < 좋음 < 매우 좋음) 를 가지면 이 설계는 정보를 낭비한다. \(k-1\) 개의 독립 모수 대신, 순서 구조를 반영한 단일 처치 효과 \(\Delta\) 로 압축하면 추정량의 자유도가 극적으로 증가한다. 이것이 Ch.5 의 비례 오즈 모형(proportional odds model)을 가져오는 동기이다.

수학적으로 핵심 변화는 이렇다. 명목형은 범주 지표 \(Y_{1j}\) 를 주변 합 \(s_j\) 로 따로 조건부 하면 각 절단점이 독립 이항 조건부 분포를 준다. 반면 순서형은 “\(\le j\)” 형태의 누적 확률 을 모형화하므로 누적 합 \(Z_{1j}\) 와 \(Z_{1,j+1}\) 이 구조적으로 종속이고, 이 상호 의존이 조건부 분포에도 남는다. 따라서 §7.4·§7.5.1 처럼 조건부 로그우도를 단순한 덧셈 형태로 분해할 수 없다 — 이 난점이 이후 가중 추정 방정식(quasi-likelihood) 으로 이어진다.

3.1 비례 오즈 모형의 조건부 버전

반응 범주가 순서를 가질 때 (예: 매우 싫음 < 싫음 < 좋음 < 매우 좋음), Ch.5의 비례 오즈 모형(proportional odds model)을 적용한다.

두 독립 다항 벡터 \(Y_1 \sim B(m_1, \pi_1)\), \(Y_2 \sim B(m_2, \pi_2)\) 에서 누적 확률 \(\gamma_{ij} = \Pr(\text{범주} \leq j)\) 가

\[ \operatorname{logit} \gamma_{1j} = \theta_j, \qquad \operatorname{logit} \gamma_{2j} = \theta_j - \Delta, \qquad j = 1, \ldots, k-1 \]

을 만족한다.

직관: \(\theta_1, \ldots, \theta_{k-1}\) 은 기저 누적 로그 오즈이고, \(\Delta\) 는 모든 절단점(cut-point)에서 동일하게 작용하는 공통 처치 효과이다. “비례 오즈”라는 이름은 \(\Delta\) 가 어떤 절단점에서든 같은 오즈비 \(\psi = e^\Delta\) 를 산출하기 때문이다.

관심 모수는 \(\Delta\) 하나이고, \(k-1\) 개의 기저 모수 \(\theta_j\) 가 장해 모수이다.

3.2 누적 합과 조건부 분포

범주별 합계 \(s_j = y_{1j} + y_{2j}\) 에서 누적 합을 구성한다:

\[ Z_{1j} = Y_{11} + \cdots + Y_{1j}, \quad S_j = s_1 + \cdots + s_j = Z_{\cdot j}. \]

\(Z_{ij} \sim B(m_i, \gamma_{ij})\) 이므로 \(S_j\) 로 조건부를 취하면 각 \(j\) 에서 비중심 초기하분포가 된다:

\[ Z_{1j} \mid S_j \sim H(m, S_j;\, \psi), \qquad \psi = e^\Delta. \]

조건부 독립의 부재 — 핵심 난점

이항 자료(§7.4.2)에서는 \(n\) 개 층의 조건부 분포가 독립이어서 조건부 로그우도가 단순한 합이었다. 그러나 순서형 반응에서는 \(Z_{1j} \mid S_j\) 들이 결합 조건부 분포로서 \(\theta\) 에도 의존한다. 즉, \(\{Z_{1j}\}\) 전체를 \(S = (S_1, \ldots, S_k)\) 로 조건부하면 \(\theta\) 가 완전히 소거되지 않는다.

이것이 명목형(§7.5.1)과의 본질적 차이이다. 명목형에서는 조건부 분포가 정확히 \(\Delta\) 만의 함수가 되었지만, 순서형에서는 그렇지 않다.

3.3 가중 추정 방정식 (Quasi-Likelihood 접근)

정확한 조건부 우도를 구성할 수 없으므로, 각 \(j\) 에서의 주변(marginal) 조건부 분포 \(Z_{1j} \mid S_j\) 를 개별적으로 활용하는 추정 방정식을 구성한다.

왜 정확한 조건부 우도는 불가능한가. \(k-1\) 개의 비중심 초기하 분포 \(Z_{1j} \mid S_j\) (\(j = 1, \ldots, k-1\)) 들이 서로 다른 표본 공간 위에 정의되고 — 각 \(S_j\) 는 그 자신 이전 \(\{S_1, \ldots, S_{j-1}\}\) 에 조건부로 의존 — 결합 분포를 닫힌 형태(closed form)로 쓰려면 \(\{\theta_j\}\) 들이 소거되지 않고 재등장한다. 즉 “모든 \(S_j\) 로 한 번에 조건부” 를 취해도 장해 모수가 완전히 사라지지 않는다. 이것이 §7.4.2 의 독립 층(stratum) 구조와 본질적으로 다른 점이며, 주변 조건부 분포만 개별적으로 사용하는 준우도 접근이 실무적으로 유일한 타협이 되는 이유이다.

\(\chi_{1j}(\psi) = E(Z_{1j} \mid S_j;\, \psi)\) 를 비중심 초기하의 조건부 기댓값이라 하면, 각 \(j\) 에서의 잔차 \(Z_{1j} - \chi_{1j}(\psi)\) 는 \(S_j\) 에 조건부로 (따라서 무조건부로도) 기댓값 0이다.

GLM의 우도 방정식이 데이터의 선형함수인 것에 착안하여, 가중 합

\[ U(\psi; Z) = \sum_{j=1}^{k-1} w_j^* \{Z_{1j} - \chi_{1j}(\psi)\} \]

를 구성한다. 가중치 \(w_j^*\) 는 \(-\partial U / \partial \Delta\) 의 기댓값이 \(U\) 의 분산과 같도록 선택한다.

직관: 이것은 Ch.9의 준우도(quasi-likelihood)와 동일한 철학이다. “정확한 우도를 모르지만, 평균과 분산 관계는 안다”는 상황에서 최적의 추정 방정식을 구성하는 것이다.

3.4 최적 가중치

가중치의 벡터 형태는

\[ w^* = V^{-1} d, \]

여기서

  • \(V\): \(Z_{11}, \ldots, Z_{1,k-1}\) 의 공분산 행렬
  • \(d_j = \partial \chi_{1j}(\psi) / \partial \Delta = \operatorname{var}(Z_{1j} \mid S_j;\, \psi)\)

\(V\) 는 “서로 다른 표본공간에서 정의된 확률변수들의 공분산”이라는 개념적 어려움이 있다. 실용적 해법으로 §7.3의 근사 초기하 공분산 행렬 \(\widetilde{\Sigma}\) 를 사용한다:

\[ \widetilde{V} = L \widetilde{\Sigma} L^T, \]

여기서 \(L\) 은 누적합을 형성하는 하삼각 행렬이다.

\(\widetilde{\Sigma}\) 의 가장 간단한 일반화역행렬을 선택하면

\[ \widetilde{\Sigma}^{-} = \frac{m_\cdot - 1}{m_\cdot} \operatorname{diag}\left\{\mu_{1j}^{-1} + \mu_{2j}^{-1}\right\}, \]

여기서 \(\mu_{1j} = \chi_{1j} - \chi_{1,j-1}\) 은 처치군의 범주별 기대 도수이다.

3.5 스코어 통계량과 점근 분산

최종 추정 방정식은

\[ U(\psi; Z) = \frac{m_\cdot - 1}{m_\cdot} \sum_{j=1}^{k} (d_j - d_{j-1}) \left(\frac{1}{\mu_{1j}} + \frac{1}{\mu_{2j}}\right)(y_{1j} - \mu_{1j}), \]

또는 근사적으로

\[ U(\psi; Z) \approx \sum_{j=1}^{k-1} \frac{\zeta_j + \zeta_{j+1}}{\zeta_\cdot} (Z_{1j} - \chi_{1j}), \]

여기서 \(\zeta_j = (\mu_{1j}^{-1} + \mu_{2j}^{-1})^{-1}\) 이고 \(\zeta_\cdot = \sum_j \zeta_j\) 이다.

\(U(\hat{\psi}_c; Z) = 0\) 의 해 \(\hat{\Delta}_c\) 의 점근 분산은

\[ \operatorname{var}(\hat{\Delta}_c) \approx \zeta_\cdot \left\{\sum_{j=1}^{k-1} (\zeta_j + \zeta_{j+1}) d_j\right\}^{-1} \approx \frac{3(m_\cdot - 1)}{m_\cdot \zeta_\cdot} \left\{1 - \sum \left(\zeta_j / \zeta_\cdot\right)^3\right\}^{-1}, \]

이 분산은 무조건부 MLE의 분산과 본질적으로 같다 (Exercise 5.3).

직관: 가중치를 잘 선택하면 조건부 추정의 효율 손실이 거의 없다. 이는 순서형 자료에서도 조건부 접근이 “비싸지 않다”는 것을 의미한다 — 편향 교정의 이득은 챙기면서 효율은 거의 보존한다.

4 예제: 치즈 맛 실험 (§7.5.3)

4.1 자료 배경

Ch.5 Table 5.1의 치즈 맛 비교 실험에서 처음 두 행 (첨가물 A와 B의 비교)을 사용한다. 반응은 9개 순서 범주 (1 = 매우 나쁨 ~ 9 = 매우 좋음)이고, 각 그룹 표본 크기는 \(m_1 = m_2 = 52\) 이다.

무조건부 MLE: \(\hat{\Delta} = -3.028\) (s.e. \(\approx 0.455\)).

4.2 반복 계산 과정

Table 7.3은 \(\Delta_0 = -3.028\) 에서 시작하는 한 사이클의 계산을 보여준다.

Table 7.3 — 조건부 추정 한 사이클
\(S_j\) \(Z_{1j}\) \(d_j\) \(\chi_{1j}\) \(\mu_{1j}\) \(\zeta_j\) \(w_j^*\)
6 0 0.2842 0.3021 0.3021 0.2869 0.0615
15 0 0.8220 0.8997 0.5976 0.5579 0.1307
28 1 1.8891 2.2860 1.3863 1.2385 0.3288
46 8 3.6281 6.5994 4.3134 3.2798 0.5105
61 16 3.4019 14.5594 7.9600 3.7359 0.4485
75 24 1.9867 25.4325 10.8731 2.4285 0.3050
95 43 0.4461 43.4791 18.0466 1.7626 0.1580
103 51 0.0440 51.0462 7.5671 0.4095 0.0330
104 52 — 52.0 0.9538 0.0441 —
52.0 13.7437

각 열의 의미:

  • \(S_j\): 누적 범주 합계 \(s_1 + \cdots + s_j\)
  • \(Z_{1j}\): 처치군의 누적 관측 \(Y_{11} + \cdots + Y_{1j}\)
  • \(d_j = \operatorname{var}(Z_{1j} \mid S_j)\): 오즈비 \(\psi_0 = e^{-3.028}\) 에서의 초기하 분산
  • \(\chi_{1j} = E(Z_{1j} \mid S_j)\): 조건부 기댓값
  • \(\mu_{1j} = \chi_{1j} - \chi_{1,j-1}\): 범주별 기대 도수
  • \(\zeta_j = (\mu_{1j}^{-1} + \mu_{2j}^{-1})^{-1}\): 가중치 계산용 조화 기대 도수
  • \(w_j^* = (\zeta_j + \zeta_{j+1}) / \zeta_\cdot\): 최적 가중치

4.3 갱신 계산

스코어 통계량:

\[ U = \sum w_j^* (Z_{1j} - \chi_{1j}) = 0.2881. \]

정보량:

\[ \sum w_j^* d_j = 4.8016. \]

Newton-Raphson 갱신:

\[ \hat{\Delta}_c = -3.028 + \frac{0.2881}{4.8016} = -2.9680. \]

한 사이클 더 반복하면

\[ \hat{\Delta}_c = -2.9743, \qquad \operatorname{var}(\hat{\Delta}_c) = (4.9047)^{-1} = 0.2039, \qquad \operatorname{s.e.} = 0.4515. \]

4.4 조건부 vs 무조건부 비교

항목 무조건부 조건부
\(\hat{\Delta}\) \(-3.028\) \(-2.974\)
s.e. 0.455 0.452
차이 — 0.054 (s.e.의 12%)

차이는 표준오차의 12%에 불과하며 실질적으로 결론에 영향을 주지 않는다. 관례대로 조건부 추정값이 원점 쪽으로 수축한다.

직관: 이 실험은 각 그룹 52명으로 비교적 큰 표본이므로 조건부-무조건부 괴리가 작다. 표본이 매우 희소한 경우에만 두 방법의 차이가 실질적이 된다 — 이는 §7.4.1의 Table 7.1 (7명짜리 표)에서 이미 확인한 바이다.

5 Python 구현

5.1 순수 구현: 명목형 매칭 쌍 (Quasi-Symmetry / Bradley-Terry)

코드 
import numpy as np
from scipy.optimize import minimize
from scipy.special import expit  # logistic function

# ── §7.5.1: 명목형 매칭 쌍 ────────────────────────
# 예: k=3 범주, 비일치 쌍 관측 (Y_{ij}: 대조군=i, 처치군=j인 쌍 수)
Y = np.array([
    [0, 15, 10],   # 대조군=1: 처치군=1,2,3
    [8, 0, 12],    # 대조군=2
    [5, 7, 0],     # 대조군=3
])

k = Y.shape[0]

# m_{ij} = Y_{ij} + Y_{ji} (대칭 합계), 상삼각만
pairs = []
for i in range(k):
    for j in range(i + 1, k):
        m_ij = Y[i, j] + Y[j, i]
        y_ij = Y[i, j]  # 대조군=i인 쌍 수
        pairs.append((i, j, y_ij, m_ij))

# 조건부 로그우도: logit(pi_{ij}) = Delta_j - Delta_i
def neg_cond_loglik(delta):
    """음의 조건부 로그우도 (Delta_1=0 제약)"""
    D = np.zeros(k)
    D[1:] = delta  # Delta_1 = 0 기준
    ll = 0.0
    for i, j, y_ij, m_ij in pairs:
        if m_ij == 0:
            continue
        eta = D[j] - D[i]
        p = expit(eta)
        ll += y_ij * np.log(p + 1e-15) + (m_ij - y_ij) * np.log(1 - p + 1e-15)
    return -ll

# 최적화
res = minimize(neg_cond_loglik, x0=np.zeros(k - 1), method='Nelder-Mead')
delta_hat = np.zeros(k)
delta_hat[1:] = res.x

print("=== Quasi-Symmetry (Bradley-Terry) MLE ===")
for j in range(k):
    print(f"  Delta_{j+1} = {delta_hat[j]:.4f}")

# 적합도: Pearson X^2
X2 = 0.0
df = (k - 1) * (k - 2) // 2
for i, j, y_ij, m_ij in pairs:
    if m_ij == 0:
        continue
    eta = delta_hat[j] - delta_hat[i]
    p = expit(eta)
    mu = m_ij * p
    X2 += (y_ij - mu)**2 / (mu * (1 - p) + 1e-15)

print(f"  Pearson X^2 = {X2:.3f} (df = {df})")

5.2 순수 구현: 순서형 매칭 쌍 (가중 추정 방정식)

코드 
import numpy as np
from math import comb, exp, log

def hyper_mean_var(m1, m2, s, psi):
    """비중심 초기하 H(m, s; psi)의 평균과 분산"""
    lo = max(0, s - m2)
    hi = min(m1, s)
    # log-sum-exp로 정규화
    log_terms = []
    for t in range(lo, hi + 1):
        lt = np.log(comb(m1, t)) + np.log(comb(m2, s - t)) + t * np.log(psi)
        log_terms.append(lt)
    log_terms = np.array(log_terms)
    max_lt = np.max(log_terms)
    log_P0 = max_lt + np.log(np.sum(np.exp(log_terms - max_lt)))

    vals = np.arange(lo, hi + 1, dtype=float)
    log_probs = log_terms - log_P0
    probs = np.exp(log_probs)
    mu = np.sum(vals * probs)
    var = np.sum(vals**2 * probs) - mu**2
    return mu, var

# ── §7.5.3: 치즈 맛 실험 (A vs B) ─────────────────
# 범주 1~9, 각 그룹 52명
# 범주별 도수: group A (y1) and group B (y2) from Table 5.1 rows 1-2
y1 = np.array([0, 0, 1, 7, 8, 8, 19, 8, 1])  # additive A
y2 = np.array([6, 9, 12, 11, 7, 6, 1, 0, 0])  # additive B
m1, m2 = 52, 52
k = len(y1)

# 누적합
S = np.cumsum(y1 + y2)   # S_j
Z1 = np.cumsum(y1)        # Z_{1j}

# Newton-Raphson iteration
delta = -3.028  # 무조건부 MLE에서 시작

for iteration in range(5):
    psi = exp(delta)
    chi_list = []
    d_list = []
    mu1_list = []

    # j = 0..k-2 (k-1개 절단점)
    for j in range(k - 1):
        chi_j, var_j = hyper_mean_var(m1, m2, S[j], psi)
        chi_list.append(chi_j)
        d_list.append(var_j)

    # mu_{1j} = chi_{1j} - chi_{1,j-1}
    mu1_list = [chi_list[0]]
    for j in range(1, k - 1):
        mu1_list.append(chi_list[j] - chi_list[j - 1])
    mu1_list.append(m1 - chi_list[-1])  # 마지막 범주

    # zeta_j = (1/mu1_j + 1/mu2_j)^{-1}
    zeta_list = []
    for j in range(k):
        mu1j = max(mu1_list[j], 1e-10)
        s_j = (y1 + y2)[j]
        mu2j = max(s_j - mu1_list[j] if j < k - 1 else m2 - (m1 - chi_list[-1]) + mu1_list[-1], 1e-10)
        # mu2j from cumulative: mu2_j = s_j - mu1_j
        mu2j = max(s_j - mu1j, 1e-10)
        zeta_list.append(1.0 / (1.0/mu1j + 1.0/mu2j))
    zeta_dot = sum(zeta_list)

    # w_j = (zeta_j + zeta_{j+1}) / zeta_dot, j = 0..k-2
    w_list = []
    for j in range(k - 1):
        w_list.append((zeta_list[j] + zeta_list[j + 1]) / zeta_dot)

    # Score: U = sum w_j * (Z1_j - chi_j)
    U = sum(w_list[j] * (Z1[j] - chi_list[j]) for j in range(k - 1))

    # Information: sum w_j * d_j
    info = sum(w_list[j] * d_list[j] for j in range(k - 1))

    delta_new = delta + U / info
    print(f"Iter {iteration}: delta = {delta:.4f} -> {delta_new:.4f}, U = {U:.4f}")

    if abs(delta_new - delta) < 1e-6:
        delta = delta_new
        break
    delta = delta_new

se = 1.0 / np.sqrt(info)
print(f"\n조건부 MLE: Delta_c = {delta:.4f}, s.e. = {se:.4f}")
print(f"무조건부 MLE: Delta   = -3.028,  s.e. = 0.455")
print(f"차이: {abs(delta - (-3.028)):.4f} ({abs(delta - (-3.028))/se:.1f}% of s.e.)")

5.3 statsmodels/scipy 활용

코드 
import numpy as np
import statsmodels.api as sm
from statsmodels.miscmodels.ordinal_model import OrderedModel

# ── §7.5.1: Bradley-Terry via GLM ─────────────────
# 비일치 쌍 데이터를 이항 GLM로 적합
Y = np.array([
    [0, 15, 10],
    [8, 0, 12],
    [5, 7, 0],
])
k = Y.shape[0]

# 이항 반응과 모형 행렬 구성
y_vec = []
n_vec = []
X_rows = []

for i in range(k):
    for j in range(i + 1, k):
        m_ij = Y[i, j] + Y[j, i]
        if m_ij == 0:
            continue
        y_vec.append(Y[i, j])
        n_vec.append(m_ij)
        # logit(pi_{ij}) = Delta_j - Delta_i
        # 기준: Delta_1 = 0, 자유 모수: Delta_2, Delta_3
        row = np.zeros(k - 1)
        if j > 0:
            row[j - 1] = -1  # -Delta_i term ... 수정 필요
        if i > 0:
            row[i - 1] = 1   # +Delta_i term
        # 정확히: logit = Delta_j - Delta_i
        # Delta_1=0 제약하에 i,j >= 1
        row2 = np.zeros(k - 1)
        for idx in range(k - 1):
            col = idx + 1  # Delta_{col+1}
            if col == j:
                row2[idx] = 1
            if col == i:
                row2[idx] = -1
        X_rows.append(row2)

y_arr = np.array(y_vec, dtype=float)
n_arr = np.array(n_vec, dtype=float)
X_mat = np.array(X_rows)

# GLM (이항, 로짓 링크, 절편 없음)
model = sm.GLM(
    y_arr / n_arr, X_mat,
    family=sm.families.Binomial(),
    var_weights=n_arr
)
result = model.fit()
print("=== Bradley-Terry via GLM ===")
print(result.summary())

6 R 구현

코드 
# ── §7.5.1: Bradley-Terry ─────────────────────────
# BradleyTerry2 패키지 사용
library(BradleyTerry2)

# 예: 3개 범주(선수), 쌍별 승패
wins <- matrix(c(0, 15, 10,
                 8, 0, 12,
                 5, 7, 0), nrow = 3, byrow = TRUE)
rownames(wins) <- colnames(wins) <- c("A", "B", "C")

bt <- BTm(outcome = 1, player1 = rep(1:3, each = 3),
          player2 = rep(1:3, times = 3),
          formula = ~ player, id = "player",
          data = data.frame(
            player1 = factor(rep(c("A","B","C"), each = 3)),
            player2 = factor(rep(c("A","B","C"), times = 3))
          ))

# 수동 이항 GLM 구현
pairs <- data.frame(
  y = c(wins[1,2], wins[1,3], wins[2,3]),
  m = c(wins[1,2]+wins[2,1], wins[1,3]+wins[3,1], wins[2,3]+wins[3,2]),
  d2 = c(1, 0, -1),   # Delta_2 coefficient
  d3 = c(0, 1, 1)     # Delta_3 coefficient (logit = D_j - D_i)
)
# logit(pi_12) = D_2 - D_1 = D_2 (D_1=0)
# logit(pi_13) = D_3 - D_1 = D_3
# logit(pi_23) = D_3 - D_2
pairs$d2 <- c(1, 0, -1)
pairs$d3 <- c(0, 1, 1)

fit_bt <- glm(cbind(y, m - y) ~ d2 + d3 - 1, data = pairs,
              family = binomial(link = "logit"))
summary(fit_bt)

# ── §7.5.2: 순서형 --- 비례 오즈 조건부 추정 ─────
# VGAM 패키지로 무조건부 비례 오즈
library(VGAM)

# 치즈 맛 데이터 (A vs B, 9 범주)
y1 <- c(0, 0, 1, 7, 8, 8, 19, 8, 1)  # additive A
y2 <- c(6, 9, 12, 11, 7, 6, 1, 0, 0)  # additive B

dat <- data.frame(
  resp = factor(rep(1:9, 2), ordered = TRUE),
  group = rep(c("A", "B"), each = 9),
  count = c(y1, y2)
)

# 비례 오즈 (무조건부)
fit_po <- vglm(resp ~ group, family = cumulative(parallel = TRUE),
               weights = count, data = dat)
summary(fit_po)

7 핵심 요약

조건부 우도의 다범주 자료 3대 응용

  • 명목형 매칭 쌍: 벡터합 \(Z_\bullet\) 으로 조건부 취하면 \(k(k-1)/2\) 개의 독립 이항이 되고, \(\operatorname{logit}(\pi_{ij}) = \Delta_j - \Delta_i\) (quasi-symmetry = Bradley-Terry). 일치 쌍은 우도에 기여하지 않는다.
  • 순서형 매칭 쌍: 누적 합계 \(S_j\) 로 조건부 취하면 각 절단점에서 비중심 초기하가 되지만, 결합 분포는 \(\theta\) 에도 의존한다. 가중 추정 방정식(quasi-likelihood)으로 \(\hat{\Delta}_c\) 를 구하며, 효율 손실은 거의 없다.
  • 치즈 맛 실험: 조건부 MLE \(\hat{\Delta}_c = -2.974\) (s.e. 0.452)은 무조건부 MLE \(\hat{\Delta} = -3.028\) (s.e. 0.455)과 표준오차의 12% 이내로 차이가 있으며, 결론에 실질적 영향이 없다.

Subscribe

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