1 들어가며 — 11-0 Overview 의 식들의 깊이
Ch.11 Overview (11-0) 에서 multinomial logit 의 큰 그림 (식 11.1-11.5) 과 ICC (식 11.6-11.7) 의 직관을 제시했다. 본 sub-post 는 같은 식들을 수학적으로 깊이 도출 + 추정 식 (11.8-11.11) 까지.
“§ 11.1 = 식 (11.1-11.2) reference cell 의 정확한 도출 — Gumbel 잡음의 효용 차이가 자연스럽게 logistic cdf 를 만드는 random utility model 의 수학. Bock 의 일반 contrast (식 11.3-11.5) 는 행렬 미적분 (vec, Kronecker, elimination) 으로 모든 contrast 표기 통합. Helmert contrast 가 ordinal 의 또 다른 표현. § 11.1.1 = 잠재 변수의 두 Gumbel 차이가 logistic 분포가 되는 증명 (McCullagh-Nelder 1989) 으로 ICC 분모의 \(\pi^2/3\) 도출. § 11.1.2 = 추정은 §9.6, §10.2.4 의 marginal MLE + Fisher scoring + Gauss-Hermite 그대로 — multinomial 우도 (식 11.8) 만 다름. Score (식 11.10) 의 핵심 항 \(D(y_{ij} - P_{ij})\) 가 ‘contrast × 잔차’ 의 자연스러운 형태.”
2 § 11.1 — Multinomial Logit 의 정확한 도출
2.1 Random Utility Model 로부터의 도출
각 환자가 \(C\) 개 선택지 (범주) 마다 잠재 효용 (latent utility) \(y_{ijc}\) 를 가짐:
\[ y_{ijc} = \mu_{ijc} + \epsilon_{ijc}, \quad c = 1, 2, \ldots, C \]
- \(\mu_{ijc}\): 결정론적 효용 부분 (\(\mu_{ijc} = x_{ij}^\top \beta_c + z_{ij}^\top T_c \theta_i\)).
- \(\epsilon_{ijc}\): 무작위 효용 부분 (관측 못한 개별 차이).
선택 규칙:
\[ Y_{ij} = c \iff y_{ijc} = \max_{c'} y_{ijc'} \]
— “가장 효용 큰 선택지 선택”.
\(\epsilon_{ijc}\) 의 분포 가정에 따라 모형이 결정:
- Type I Extreme Value (Gumbel) 잡음 → multinomial logit (closed form).
- 다변량 정규 잡음 → multinomial probit (적분 필요, 복잡).
- Generalized Extreme Value 잡음 → nested logit (계층적 선택 모형).
Gumbel 의 결정적 장점: \(P(Y = c) = \frac{\exp(\mu_c)}{\sum \exp(\mu_h)}\) 이 closed form.
증명 핵심:
\[ P(Y = c) = P(y_c > y_h \text{ for all } h \neq c) = P(\epsilon_c - \epsilon_h > \mu_h - \mu_c \text{ for all } h) \]
Gumbel 잡음의 cumulative distribution function:
\[ F_\epsilon(z) = \exp(-\exp(-z)) \]
Independence + 분포 형태로 적분:
\[ P(Y = c) = \int_{-\infty}^\infty f_\epsilon(\epsilon_c) \prod_{h \neq c} F_\epsilon(\epsilon_c + \mu_c - \mu_h) \, d\epsilon_c = \frac{\exp(\mu_c)}{\sum_h \exp(\mu_h)} \]
이것이 식 (11.1-11.2) 의 정확한 출처. Random utility + Gumbel = multinomial logit.
경제학·마케팅 표준: 선택 모형 (discrete choice model) 의 토대. 통근 수단 선택, 제품 brand 선택, 직업 선택 등.
2.2 식 (11.1-11.2) — Reference Cell Formulation
\(c = 1\) 을 reference 로 두면 \(\beta_1 = 0, T_1 = 0\) → \(\mu_{ij1} = 0 \Rightarrow \exp(\mu_{ij1}) = 1\):
\[ p_{ijc} = \Pr(Y_{ij} = c \mid \theta) = \frac{\exp(z_{ijc})}{1 + \sum_{h=2}^C \exp(z_{ijh})} \quad (c = 2, \ldots, C) \tag{11.1} \]
여기서 \(z_{ijc} = x_{ij}^\top \beta_c + z_{ij}^\top T_c \theta_i\).
\[ p_{ij1} = \frac{1}{1 + \sum_{h=2}^C \exp(z_{ijh})} \tag{11.2} \]
확률 합 검증:
\[ \sum_{c=1}^C p_{ijc} = \frac{1 + \sum_{h=2}^C \exp(z_{ijh})}{1 + \sum_{h=2}^C \exp(z_{ijh})} = 1 \]
\(\mu_{ij1} = 0\) 으로 reference 를 고정하는 것은 위치 식별 가능성 의 결과:
- 효용은 차이만 의미 — \(y_{ijc} \to y_{ijc} + k\) 는 모든 효용을 같이 이동, 선택 규칙 (최대값) 무관.
- 따라서 한 범주의 효용을 0 으로 고정하면 다른 범주의 효용이 식별 가능.
- Reference 선택은 scale 이 아닌 origin — 어떤 범주를 0 으로 두든 결과는 같음 (모수 회전).
다른 reference 선택:
- 다른 범주 \(c^* \neq 1\) 을 reference 로 두면 모수 변환:
- \(\beta_c^{\text{new}} = \beta_c - \beta_{c^*}\).
- 적합도, 예측 동일.
Reference 선택의 실용 권고:
- 가장 흔한 범주: \(\sum_i \mathbb{1}(Y_i = c)\) 가 최대인 범주. 표본 효율.
- 임상적 기준: “control” 또는 “no event” 범주 (예: no care, undecided, none).
- 해석 단순성: 다른 범주가 reference 와 비교될 때 임상적 의미 명확.
3 Bock (1972) 의 General Contrast Formulation
3.1 식 (11.3-11.4) — 행렬 미적분으로 표현
\(C\) 개 모든 범주에 대해 대칭적 표기:
\[ p_{ijc} = \frac{\exp(z_{ijc})}{\sum_{h=1}^C \exp(z_{ijh})} \quad (c = 1, 2, \ldots, C) \tag{11.3} \]
(식 11.1-11.2 와 동일하지만 reference 0 을 명시적으로 도입하지 않음.)
\[ z_{ijc} = x_{ij}^\top \Gamma d_c + (z_{ij}^\top \otimes \theta_i^\top) J_{r^*}^\top \Lambda d_c \tag{11.4} \]
표기:
- \(D = [d_1, d_2, \ldots, d_C]\): \((C-1) \times C\) contrast matrix. 각 column 이 한 범주에 적용되는 contrast vector.
- \(\Gamma\): \(p \times (C-1)\) 회귀 계수 matrix. 각 column 이 한 contrast 의 회귀 계수.
- \(\Lambda = [v(T_1), v(T_2), \ldots, v(T_{C-1})]\): \(r^* \times (C-1)\) random effects 모수 matrix. 각 column 이 한 contrast 의 Cholesky factor 의 free elements (\(r^* = r(r+1)/2\)).
- \(J_{r^*}\): Magnus (1988) elimination matrix. Random effects 행렬 미분 (§ 9.6.2 의 식 9.54 와 같은 도구).
- \(\otimes\): Kronecker product.
식 (11.4) 가 복잡해 보이지만 핵심 두 단계:
- 회귀 부분 \(x_{ij}^\top \Gamma d_c\):
- \(\Gamma d_c\) 가 contrast \(c\) 의 회귀 계수 (\(p \times 1\)).
- \(x^\top \Gamma d_c\) 가 그 contrast 의 선형 예측자.
- Random effects 부분 \((z_{ij}^\top \otimes \theta_i^\top) J_{r^*}^\top \Lambda d_c\):
- \(\Lambda d_c\) 가 contrast \(c\) 의 Cholesky factor 의 free elements (\(r^* \times 1\)).
- \(J_{r^*}^\top \Lambda d_c\) 가 elimination 으로 free elements 만 추출 후 vec 형태로 변환.
- \((z^\top \otimes \theta^\top)\) 가 \(z'T_c \theta\) 의 전개 — random effects design 과 효과의 곱.
왜 이렇게 복잡한 표기가 필요한가:
- 일반 contrast 표현: reference cell 만 아니라 Helmert, polynomial, 사용자 정의 contrast 모두 처리.
- 추정 알고리즘 (식 11.10) 이 행렬 형태로 단순화.
- 식 (11.4) 가 한 식으로 모든 변형 포함 — 모든 special case 의 공통 framework.
3.2 식 (11.5) — Random Intercept 단순화
\(z_{ij} = 1\) (random intercept), \(\theta_i\) 가 스칼라:
\[ z_{ijc} = x_{ij}^\top \Gamma d_c + \Lambda d_c \theta_i \tag{11.5} \]
여기서 \(\Lambda = [\sigma_1, \sigma_2, \ldots, \sigma_{C-1}]\) (\(1 \times (C-1)\) vector, 각 contrast 의 random intercept SD).
식 (11.5) 의 \(\Lambda d_c \theta_i\) 의 의미:
- \(\Lambda d_c\) 가 contrast \(c\) 의 random intercept SD (\(1 \times 1\) 스칼라).
- \(\theta_i\) 가 표준 정규 random effect.
- → 각 contrast 마다 다른 SD 로 scale 된 random intercept.
Reference cell 사례 (D = diagonal-like matrix in 11-0):
\[ D = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix} \]
- \(d_1 = (0, 0, \ldots, 0)^\top\) → reference 의 효용 = 0.
- \(d_c = (\delta_{c-1, 1}, \delta_{c-1, 2}, \ldots)\) — \(c-1\) 번째 column 이 1, 나머지 0.
→ \(\Gamma d_c\) 가 \(\Gamma\) 의 \((c-1)\) 번째 column = \(\beta_c\). → \(\Lambda d_c\) 가 \(\Lambda\) 의 \((c-1)\) 번째 column = \(\sigma_c\).
식 (11.1-11.2) 의 단순 표기로 reduce.
3.3 Helmert Contrast — Ordinal 의 또 다른 표현
\(C = 4\) 일 때 Helmert contrast matrix:
\[ D = \begin{bmatrix} -3/4 & 1/4 & 1/4 & 1/4 \\ 0 & -2/3 & 1/3 & 1/3 \\ 0 & 0 & -1/2 & 1/2 \end{bmatrix} \]
의미:
- 첫 contrast: 범주 1 vs (2, 3, 4) — “첫 범주 vs 이후 모두”.
- 둘째 contrast: 범주 2 vs (3, 4) — “둘째 범주 vs 이후 모두”.
- 셋째 contrast: 범주 3 vs 4 — “셋째 vs 마지막”.
각 row 의 합 = 0 (zero-sum contrast). Scale 은 contrast coefficient 의 차이로 정규화 (1).
응답이 본질적으로 순서 가 있을 때 Helmert contrast 는 ordinal 모형의 대안:
- 순차적 비교 — “각 범주 vs 그 다음 범주들의 평균”.
- 의학·생존 분야의 자연 순서.
Continuation-Ratio Logit (Ten Have & Uttal 1994) 와의 비교:
- Continuation-ratio: \(\log[P(Y = c) / P(Y > c)]\) — “\(c\) 에 머묾 vs \(c\) 너머 진행” (probability 비율).
- Helmert + multinomial: \(\log[P(Y = c) / P(Y \in \{c+1, \ldots, C\})\) 같은 형태] — 비슷하지만 logit 척도.
차이점:
- Continuation-ratio: probabilities 사이의 비율.
- Helmert (multinomial): logit 자체에 contrast 적용.
- 비슷한 직관, 다른 수학적 형태.
언제 Helmert 사용:
- 응답이 순서 있지만 비례 오즈 가정 위반 (§ 10.2.1 partial PO 의 대안).
- 순차적 진행 (예: “치료 단계의 진행”, “학년의 진급”) 모형화.
- Sparse 데이터 — Helmert 가 reference cell 보다 추정 안정성 좋을 수 있음.
4 § 11.1.1 — ICC 의 Type I Extreme Value 도출
4.1 식 (11.6) — 잠재 변수 표현
\[ y_{ijc} = x_{ij}^\top \beta_c + \sigma_c \theta_i + \epsilon_{ijc} \quad (c = 1, 2, \ldots, C) \tag{11.6} \]
- \(y_{ijc}\): 범주 \(c\) 의 잠재 효용 (extremal concept).
- \(\sigma_c\): 범주 \(c\) 의 random intercept SD.
- \(\epsilon_{ijc} \sim \text{Type I Extreme Value (Gumbel)}\), 독립.
4.2 식 (11.7) — Reference Cell 표현
\(\beta_1 = 0, \sigma_1 = 0\) 으로 두고 reference 대비:
\[ y_{ijc} = x_{ij}^\top \beta_c + \sigma_c \theta_i + (\epsilon_{ijc} - \epsilon_{ij1}) \quad (c = 2, \ldots, C) \tag{11.7} \]
선택 규칙: \(Y = c \iff y_{ijc} > 0\) for all \(c > 1\) + \(y_{ij1} = 0\) 이 최대일 때 reference 선택.
(식 11.7 의 좌변은 reference 대비 차이 — original \(y_{ijc} - y_{ij1}\).)
식 (11.7) 의 잡음 부분 \(\epsilon_{ijc} - \epsilon_{ij1}\) 의 분포 도출.
Gumbel 분포의 정의:
\[ F_{\text{Gumbel}}(z) = \exp(-\exp(-z)) \]
평균 = \(\gamma \approx 0.577\) (Euler-Mascheroni 상수), 분산 = \(\pi^2 / 6\).
두 독립 Gumbel 의 차이의 분포 (McCullagh-Nelder 1989):
\(X, Y \sim \text{i.i.d. Gumbel}\) 이면 \(Z = X - Y\) 의 cdf:
\[ F_Z(z) = P(X - Y \leq z) = \int_{-\infty}^\infty F_X(y + z) f_Y(y) \, dy \]
대입 후 적분:
\[ F_Z(z) = \int_{-\infty}^\infty \exp(-\exp(-(y+z))) \cdot \exp(-y) \exp(-\exp(-y)) \, dy \]
치환 \(u = \exp(-y)\):
\[ F_Z(z) = \int_0^\infty \exp(-u \exp(-z)) \exp(-u) \, du = \int_0^\infty \exp(-u(1 + \exp(-z))) \, du \]
\[ = \frac{1}{1 + \exp(-z)} = \Psi(z) \]
→ 표준 logistic 분포. 분산 \(\pi^2/3\).
이 결과의 의미:
- Gumbel 잡음의 효용 차이가 logistic.
- → multinomial logit 의 link 가 logistic 인 정확한 수학적 근거.
- 두 independent Gumbel 분산 합이 \(2 \cdot \pi^2/6 = \pi^2/3\) — logistic 분산과 일치 (독립 차이의 분산 보존).
4.3 \(C-1\) 개 ICC
식 (11.7) 의 잠재 변수 분산 분해 (range \(c = 2, \ldots, C\)):
\[ V(y_{ijc} - y_{ij1} \mid x) = \sigma_c^2 + V(\epsilon_{ijc} - \epsilon_{ij1}) = \sigma_c^2 + \pi^2/3 \]
ICC:
\[ \widehat{r}_c = \frac{\widehat\sigma_c^2}{\widehat\sigma_c^2 + \pi^2/3} \quad (c = 2, 3, \ldots, C) \]
→ \(C-1\) 개 ICC.
각 범주별 ICC 가 다르면:
- \(r_c\) 큼: 범주 \(c\) vs reference 의 odds 가 환자별 매우 다름. 환자 segmentation 효과적.
- \(r_c\) 작음: 범주 \(c\) 의 응답이 환자 차이로 잘 설명 안 됨 (시점 변동이 더 중요).
Ordinal vs Nominal ICC:
- Ordinal (Ch.10): 단일 ICC — 모든 cumulative logit 에 같은 환자 분산.
- Nominal (Ch.11): \(C-1\) 개 ICC — 각 contrast 별 환자 분산 다름.
예시 — 의료 이용 (no care / outpatient / inpatient):
- \(r_2\) (outpatient vs no care): 0.15 — 외래 결정에 환자 차이 작음.
- \(r_3\) (inpatient vs no care): 0.42 — 입원 결정에 환자 차이 큼.
해석:
- 외래 진료는 일반적인 의료 접근성 (환자 무관).
- 입원은 만성 질환, 보험, 거주지 등 환자 고유 특성 영향.
→ Nominal 모형의 고유 정보 — ordinal 에서 절대 못 얻음.
5 § 11.1.2 — Parameter Estimation 의 깊이
5.1 식 (11.8) — Multinomial Conditional Likelihood
\(Y_i\) 가 환자 \(i\) 의 시점별 응답 vector. Conditional independence 가정:
\[ \ell(Y_i \mid \theta) = \prod_{j=1}^{n_i} \prod_{c=1}^C (p_{ijc})^{y_{ijc}} \tag{11.8} \]
표기:
- \(y_{ijc} = 1\) if \(Y_{ij} = c\), 0 otherwise.
- 각 시점 \(j\) 에서 \(C\) 개 indicator 중 정확히 하나만 1.
- \(p_{ijc}\) = 식 (11.3) 의 multinomial probability.
식 (11.8) 의 형태가 Bernoulli (식 9.39), multinomial ordinal (식 10.16) 와 동일 구조:
- 곱 \(\prod_j \prod_c\).
- Indicator \(y_{ijc}\) 가 한 범주만 살림.
- Conditional likelihood = “관측 범주 확률의 곱”.
§ 9.6 → § 10.2.4 → § 11.1.2 의 연속성:
- Bernoulli: \(C = 2\) 의 multinomial 특수.
- Cumulative ordinal: cell probability 가 인접 cumulative 의 차이.
- Multinomial nominal: cell probability 가 식 (11.3) 의 직접 표현.
세 모형의 conditional likelihood 가 같은 구조 — 같은 추정 알고리즘이 자동 적용.
5.2 식 (11.9) — Marginal Log-Likelihood
\[ \log L = \sum_{i=1}^N \log h(Y_i) = \sum_{i=1}^N \log \int_\theta \ell(Y_i \mid \theta) g(\theta) \, d\theta \tag{11.9} \]
\(g(\theta)\) = 다변량 표준 정규 (\(\mathcal{N}_r(0, I)\)).
식 (11.9) 의 marginal MLE 가 이항·순서형의 직접 일반화. 적분 차원, 분포, 최대화 절차 모두 동일.
차이는:
- \(\ell(Y_i \mid \theta)\) 의 형태 (multinomial 곱).
- Random effects 차원이 더 클 수 있음 (\(C-1\) 개 random intercept).
→ Gauss-Hermite quadrature 의 차원 폭발 (\(Q^r\), § 9.6.3) 이 명목 모형에서 더 심함.
5.3 식 (11.10) — Score 의 일반 형태
\(\Delta\) = \(\Gamma\) 또는 \(\Lambda\) 둘 중 하나 (모수 행렬).
\[ \frac{\partial \log L}{\partial \Delta^\top} = \sum_{i=1}^N h_i^{-1} \int_\theta \left[ \sum_{j=1}^{n_i} D(y_{ij} - P_{ij}) \otimes \partial \Delta \right] \ell_i \, g(\theta) \, d\theta \tag{11.10} \]
표기:
- \(y_{ij}\): \(C \times 1\) indicator vector (\(Y_{ij} = c\) 인 위치만 1, 나머지 0).
- \(P_{ij}\): \(C \times 1\) predicted probability vector — 식 (11.3) 으로 계산.
- \(D\): \((C-1) \times C\) contrast matrix.
- \(\partial \Delta\): 모수별 미분 (식 11.11).
\[ \partial \Gamma = x_{ij}^\top, \quad \partial \Lambda = [J_{r^*}(\theta \otimes z_{ij})]^\top \tag{11.11} \]
\(\Gamma\) 미분 (회귀 계수): \(x_{ij}^\top\) — 표준 회귀 score 의 형태.
\(\Lambda\) 미분 (random effects): 행렬 미분 toolkit (§ 9.6.2 의 식 9.54 와 동일 형태).
식 (11.10) 의 핵심 항 \(D(y_{ij} - P_{ij})\) 분해:
- \(y_{ij} - P_{ij}\): \(C \times 1\) 잔차 vector — “관측 indicator - 예측 확률”.
- \(D\): \((C-1) \times C\) contrast matrix.
- \(D(y_{ij} - P_{ij})\): \((C-1) \times 1\) contrast 별 잔차 — 각 비교의 잔차.
→ Score = 잔차의 contrast × 모수 미분의 가중 평균 (사후 분포로 weighting).
Reference cell 사례 (D = diagonal-like):
- \(D(y_{ij} - P_{ij})\) 의 \(c\) 번째 원소 = \(y_{ijc} - P_{ijc}\) (단, \(c = 2, \ldots, C\)).
- → 단순 잔차 형태.
General contrast 사례 (Helmert 등):
- \(D(y_{ij} - P_{ij})\) 의 \(k\) 번째 원소 = contrast \(k\) 의 잔차.
- → 일반화된 잔차 — contrast 가 어떤 비교에 가중치를 두는지에 따라.
Bernoulli (식 9.44) 와의 연결:
- \(C = 2\) + reference cell → \(D(y_{ij} - P_{ij})\) 가 \((y_{ij1} - P_{ij1})\) 한 원소 = -(잔차).
- → § 9.6 의 식 (9.44) 의 잔차 형태와 일치.
Ordinal cumulative (식 10.19) 와의 차이:
- Cumulative ordinal: 인접 cumulative 차이 형태.
- Multinomial nominal: 직접 잔차 형태.
- 두 모형이 같은 framework 위에서 다른 곱셈/차이 구조.
5.4 Fisher Scoring + Gauss-Hermite
§ 9.6 식 (9.47), § 10.2.4 식 (10.22) 와 같은 Fisher scoring:
\[ \Theta_{l+1} = \Theta_l + I(\Theta_l)^{-1} \frac{\partial \log L}{\partial \Theta_l} \]
정보 행렬:
\[ I(\Theta_l) = \sum_{i=1}^N h_i^{-2} \frac{\partial h_i}{\partial \Theta_l} \left( \frac{\partial h_i}{\partial \Theta_l} \right)^\top \]
추정 알고리즘이 이항 (§ 9.6) → 순서형 (§ 10.2.4) → 명목 (§ 11.1.2) 으로 가면서 모수 차원만 증가, 알고리즘 구조 동일:
| Chapter | 회귀 모수 | Random effects 모수 |
|---|---|---|
| § 9 (Bernoulli) | \(p\) | \(r(r+1)/2\) |
| § 10 (Cumulative ordinal) | \(p + (C-2)\) (절단점) | \(r(r+1)/2\) |
| § 11 (Multinomial) | \(p(C-1)\) | \(r(r+1)/2 \cdot (C-1)\) |
명목 모형의 모수 폭발 — 4 범주 + 5 covariate + 2 random effects = \(5 \cdot 3 = 15\) 회귀 + \(3 \cdot 3 = 9\) random effects = 24 모수.
Gauss-Hermite quadrature 부담:
- 적분 차원 = random effects 의 총 자유 모수 수.
- 명목 모형: \(r \cdot (C-1)\) 차원 잠재 vector.
- \(Q = 5\), \(r = 1\), \(C = 4\) → \(5^3 = 125\) 점.
- \(r = 2\), \(C = 4\) → \(5^6 = 15625\) 점.
→ 5+ 범주 + 다중 random effects 면 Bayesian (Stan) 또는 Laplace approximation 권고.
적합 stability 의 issue:
- 모수 많음 + sparse 데이터 → 일부 모수 수렴 어려움.
- 일부 범주의 표본이 매우 작으면 그 범주의 회귀 계수 추정 불안정.
- 실무에서는 가장 작은 범주 표본 ≥ 30 정도 권고.
6 응용 분야
| 분야 | \(C\) | Reference cell 권고 |
|---|---|---|
| 정치 분석 | 3-5 (republican / democrat / independent / undecided / other) | undecided 또는 가장 작은 |
| 의료 이용 | 3-4 (no care / outpatient / inpatient / emergency) | no care |
| 통근 수단 | 3-5 (car / bus / subway / bike / walk) | 가장 흔한 |
| 정신과 진단 | 3-5 (없음 / 우울 / 불안 / 기분조절 / 정신증) | 없음 |
| 직업 분류 | 4-6 (manager / professional / clerk / labor / unemployed) | 가장 흔한 |
| 마케팅 brand | 3-5 (brand A / B / C / none) | none |
7 코드 예시
7.1 Step 1: Random Utility Model 의 직접 시뮬레이션
import numpy as np
import pandas as pd
def simulate_random_utility(n: int, beta_c: dict, sigma_c: dict,
n_categories: int = 4,
seed: int = 2026) -> pd.DataFrame:
"""McFadden (1974) 의 random utility model 직접 구현.
각 범주의 잠재 효용 = 결정론적 부분 + Gumbel 잡음.
선택 = 최대 효용의 범주.
"""
rng = np.random.default_rng(seed)
# 공변량
x = rng.normal(size=n)
# 각 범주의 효용 계산 (reference c=1: beta=0, sigma=0)
utilities = np.zeros((n, n_categories))
for c in range(1, n_categories): # c = 2, 3, ...
beta = beta_c.get(c, 0.0)
sigma = sigma_c.get(c, 0.0)
upsilon = sigma * rng.normal(size=n)
# Gumbel 잡음 (-log(-log(uniform)))
eps = -np.log(-np.log(rng.uniform(size=n)))
utilities[:, c] = beta * x + upsilon + eps
# Reference 범주
utilities[:, 0] = -np.log(-np.log(rng.uniform(size=n)))
# 최대 효용의 범주 선택
Y = utilities.argmax(axis=1) + 1 # 1, 2, ..., C
return pd.DataFrame({"x": x, "Y": Y})
# 4 범주 시뮬레이션
df = simulate_random_utility(n=5000,
beta_c={2: 0.5, 3: 1.0, 4: -0.5},
sigma_c={2: 0.0, 3: 0.0, 4: 0.0}, # single-level (random effects 없음)
n_categories=4)
print(f"전체: {len(df)}")
print(f"\n범주별 비율:")
print(df['Y'].value_counts().sort_index() / len(df))\(\beta_c\) 효과:
- \(\beta_2 = 0.5\) → 범주 2 의 효용 ↑ → 비율 ↑.
- \(\beta_3 = 1.0\) (가장 큼) → 범주 3 가 가장 큰 비율.
- \(\beta_4 = -0.5\) → 범주 4 의 비율 ↓.
이 패턴이 시뮬레이션 결과에서 확인되면 random utility model 의 정확성 검증.
7.2 Step 2: Multinomial Logit 의 Closed Form 검증
import numpy as np
from scipy.special import softmax
def multinomial_probabilities_closed_form(x: np.ndarray, beta_c: dict,
n_categories: int = 4) -> np.ndarray:
"""식 (11.3) 의 closed form.
Random utility model 의 적분 결과 — Gumbel 잡음 가정으로 도출.
"""
n = len(x)
z = np.zeros((n, n_categories))
for c in range(1, n_categories):
beta = beta_c.get(c, 0.0)
z[:, c] = beta * x
# Reference c=1: z=0
return softmax(z, axis=1)
# 시뮬레이션과 closed form 비교
beta_c = {2: 0.5, 3: 1.0, 4: -0.5}
x_test = np.array([-2, -1, 0, 1, 2])
# Closed form 확률
P_closed = multinomial_probabilities_closed_form(x_test, beta_c, n_categories=4)
print("Closed form 확률 (식 11.3):")
print(P_closed.round(3))
# 시뮬레이션 (큰 표본) 으로 검증
n_sim = 100000
x_sim = np.repeat(x_test, n_sim // len(x_test))
df_sim = simulate_random_utility(n=len(x_sim),
beta_c=beta_c,
sigma_c={2: 0.0, 3: 0.0, 4: 0.0},
n_categories=4)
df_sim['x'] = x_sim
# x 별 빈도
P_sim = df_sim.groupby('x')['Y'].value_counts(normalize=True).unstack(fill_value=0)
print("\n시뮬레이션 비율:")
print(P_sim.round(3))Closed form 확률과 시뮬레이션 비율이 거의 일치 (시뮬레이션 표본 변동 내).
이것이 random utility + Gumbel 잡음 = multinomial logit 의 직접 검증. 식 (11.1-11.2) 의 정확한 도출 시연.
7.3 Step 3: 두 Gumbel 차이가 Logistic 인지 검증
import numpy as np
import matplotlib.pyplot as plt
def sample_gumbel(n: int, seed: int = 2026) -> np.ndarray:
"""Type I Extreme Value (Gumbel) 표본 생성."""
rng = np.random.default_rng(seed)
return -np.log(-np.log(rng.uniform(size=n)))
# 두 독립 Gumbel 의 차이
n = 100000
g1 = sample_gumbel(n, seed=1)
g2 = sample_gumbel(n, seed=2)
diff = g1 - g2
# 분포 통계
print(f"두 Gumbel 의 차이 통계:")
print(f" 평균: {diff.mean():.4f} (예상: 0)")
print(f" 분산: {diff.var():.4f} (예상: pi^2/3 = {np.pi**2/3:.4f})")
print(f" 왜도: {((diff - diff.mean())**3).mean() / diff.std()**3:.4f} (예상: 0, logistic 대칭)")
# 분포 시각화
from scipy.stats import logistic
fig, ax = plt.subplots(figsize=(8, 5))
ax.hist(diff, bins=100, density=True, alpha=0.5, label="Gumbel 1 - Gumbel 2 (시뮬)")
z_range = np.linspace(-10, 10, 200)
ax.plot(z_range, logistic.pdf(z_range), "r-", linewidth=2,
label="표준 logistic pdf (이론)")
ax.set_xlabel("z")
ax.set_ylabel("Density")
ax.set_title("두 독립 Gumbel 의 차이 = 표준 Logistic")
ax.legend()
plt.tight_layout()분포 통계:
- 평균 ≈ 0 (logistic 평균).
- 분산 ≈ \(\pi^2/3\) (logistic 분산).
- 왜도 ≈ 0 (logistic 대칭, Gumbel 비대칭이지만 차이는 대칭).
히스토그램이 logistic pdf 와 일치 → 두 Gumbel 차이 = logistic 의 시각적 확인.
이 결과가 식 (11.7) 의 ICC 분모 \(\pi^2/3\) 의 정확한 출처. 두 Gumbel 차이가 logistic 이라는 사실 위에 nominal mixed-effects 의 ICC 모든 식이 서 있음.
7.4 Step 4: Helmert Contrast 의 직접 적용 (R)
library(nnet)
# 4 범주 ordinal 데이터 (Helmert 적용 시연)
set.seed(2026)
n <- 1000
x <- rnorm(n)
# 시뮬레이션 — 약한 ordinal 패턴 (cumulative log-odds 비슷)
y_latent <- 0.5 * x + rlogis(n)
gamma <- c(-1.5, 0, 1.5)
y <- cut(y_latent, breaks = c(-Inf, gamma, Inf), labels = 1:4, ordered = FALSE)
df <- data.frame(y = y, x = x)
# Helmert contrast 설정
contrasts(df$y) <- contr.helmert(4)
print(contrasts(df$y))
# 출력:
# 첫 contrast: -1, 1, 0, 0 (범주 1 vs 2)
# 둘째 contrast: -1, -1, 2, 0 (범주 1, 2 vs 3)
# 셋째 contrast: -1, -1, -1, 3 (범주 1, 2, 3 vs 4)
# Bock (1975) 의 Helmert 정의와 약간 다름 (R 의 표준 정의)
# Multinomial logit 적합 (Helmert contrast 활용)
fit_helmert <- multinom(y ~ x, data = df, trace = FALSE)
summary(fit_helmert)R 의 contr.helmert 와 Bock (1975, p. 196) 의 Helmert contrast 가 약간 다름:
- R: integer-valued contrasts.
- Bock: scaled contrasts (sum-to-zero, scale to unity).
두 정의 모두 같은 비교 (각 범주 vs 이후 평균) 를 표현하지만 scale 이 다름. 결과 해석 시 contrast scale 고려 필수.
Bock 의 정의가 multinomial 표기에 더 자연 — 식 (11.4) 의 \(\Lambda d_c\) 의 자동 normalization.
R 에서 Bock 정의 직접 사용:
D_bock <- matrix(c(
-3/4, 1/4, 1/4, 1/4,
0, -2/3, 1/3, 1/3,
0, 0, -1/2, 1/2
), nrow=3, byrow=TRUE)
contrasts(df$y) <- t(D_bock)Mixed-effects multinomial with Helmert — mclogit::mblogit 가 contrast 옵션 지원, 또는 직접 contrast matrix 적용 후 적합.
8 관련 주제
선행 지식
- Ch.11 Overview — Nominal GLMM 의 큰 그림 (식 11.1-11.5)
- § 9.6 추정 — Marginal MLE + Gauss-Hermite (Ch.11 에 직접 적용)
- § 9.5.5 c-log-log — Type I extreme value 잠재 잡음
- § 10.2.4 추정 — Cumulative ordinal MLE
- § 9.6.2 Cholesky 미분 — 식 (9.54) 의 nominal 적용
후속 주제 (Ch.11 sub-posts)
- § 11.2 — Health Services Research Example (3 범주 의료 이용 nominal, MHRP 데이터 재방문)
- § 11.3 — Competing Risk Survival Models
- § 11.3.1: Waiting for Organ Transplantation
- § 11.4 — Ch.11 Summary
관련 개념
- McFadden (1974) — Random utility model 원전 (Nobel 2000)
- Bock (1972, 1975) — Multinomial logit + extremal concept + Helmert contrast 정립
- Magnus (1988) — 행렬 미적분 toolkit (식 11.4 의 vec, Kronecker, elimination)
- McCullagh & Nelder (1989) — 두 Gumbel 차이 = logistic 의 정리
- Maddala (1983) — Type I extreme value 분포 표준 reference
- Ten Have & Uttal (1994) — Continuation-ratio logit (Helmert 와 비교)
- Hedeker (2003) — Mixed-effects multinomial general contrast 정립
- Hartzel et al. (2001) — Clustered ordinal/nominal 통합
- § 10.2.1 Partial PO — Helmert 의 ordinal 대안 토대