1 들어가며 — § 9.6 의 직접 확장

§ 9.6 (이항 GLMM 추정) 의 framework — marginal MLE + Gauss-Hermite quadrature + Fisher scoring + Empirical Bayes — 가 ordinal mixed-effects 에 거의 그대로 적용된다. § 10.2.4 는 그 일반화의 세부.

§ 9.6 와 § 10.2.4 의 비교

항목	§ 9.6 (이항)	§ 10.2.4 (순서형)
Conditional likelihood	Bernoulli 곱 (식 9.39)	Multinomial 곱 (식 10.16)
Censoring 처리	없음	식 (10.17) — survival 응용
자유 모수	$\beta, \sigma_v$	$\beta, T, \gamma_c, \alpha_c, \tau$
$c$-가변 모수	없음	$\gamma_c$ ($C-1$ 개), $\alpha_c$ ($h \times (C-1)$ 개)
Score 형태	잔차 × 공변량 (식 9.44)	인접 cumulative 차이 (식 10.19)
적분	Gauss-Hermite (식 9.55)	동일 (식 10.18)
Update 알고리즘	Fisher scoring (식 9.47)	동일 (식 10.22)

→ 새로 배울 것: 절단점 미분, multinomial 우도, censoring 처리. 나머지는 § 9.6 그대로.

한 줄 요약

“§ 10.2.4 = ordinal 의 cell 확률을 인접 cumulative 의 차이로 정의 (식 10.14) → multinomial 우도 (식 10.16) → 우도를 random effects 로 적분 → marginal MLE. Score 는 두 종류 미분: c-무관 모수 ($\beta$, $T$) 는 §9.6 과 같은 형태로 인접 pdf 차이가 가중치, c-가변 모수 ($\gamma_c, \alpha_c$) 는 Kronecker delta $\delta_{cc'}$ 로 해당 cumulative 만 영향. Scaling $\tau$ 는 chain rule 의 추가 한 단계. Fisher scoring (식 10.22) 과 Gauss-Hermite quadrature 가 §9.6 의 toolkit 그대로 적용. 새로운 알고리즘이 아니라 인접 cumulative 의 차이라는 ordinal 특화 한 단계만 추가.”

2 § 10.2.4 — Cumulative Ordinal 의 Marginal MLE

2.1 식 (10.14) — Cell 확률

Conditional cell probability

랜덤 효과 $\theta$ 가 주어진 상태에서 응답이 범주 $c$ 에 떨어질 확률:

\[ \Pr(Y_{ij} = c \mid \theta) = P_{ijc} - P_{ij, c-1} \tag{10.14} \]

$P_{ijc}$ 는 cumulative 확률 (식 10.1):

\[ P_{ijc} = \frac{1}{1 + \exp(-z_{ijc})} = \Psi(z_{ijc}) \quad \text{(logit)} \]

경계: $P_{ij0} = 0$ ($\gamma_0 = -\infty$), $P_{ijC} = 1$ ($\gamma_C = \infty$).

따라서:

$P(Y_{ij} = 1) = P_{ij1} - 0 = P_{ij1}$.
$P(Y_{ij} = C) = 1 - P_{ij, C-1}$.

직관 — “범주 확률은 인접 cumulative 의 차이”

Cumulative 확률 $P_{ijc} = P(Y \leq c)$ 가 정의된 후, 범주 $c$ 자체의 확률은 인접 cumulative 의 차이 — 잠재 변수 framework 에서 자연스러움 (Figure 10.1 의 두 절단점 사이 면적).

3 범주 사례:

$P(Y = 1) = P_{ij1} - 0 = \Psi(\gamma_1 - x^\top \beta)$.
$P(Y = 2) = P_{ij2} - P_{ij1} = \Psi(\gamma_2 - x^\top \beta) - \Psi(\gamma_1 - x^\top \beta)$.
$P(Y = 3) = 1 - P_{ij2} = 1 - \Psi(\gamma_2 - x^\top \beta)$.

이 단순한 차이가 ordinal 모형의 모든 우도·미분의 토대.

왜 cumulative 표기를 쓰는가:

범주 확률을 직접 모형화하려면 각 범주에 별도의 회귀가 필요 — 명목 모형 (Ch.11) 의 길.
Cumulative 표기는 단조 누적의 자연스러움을 활용 — 절단점이 위치 정보, 회귀 계수가 경사 정보.
→ 절약적 (proportional odds) + 해석 직관적.

2.2 식 (10.15) — General $z_{ijc}$

통합된 선형 예측자

본 sub-post 는 가장 일반적 모형 (location-scale, partial PO 포함, § 10.2.2 의 식 10.8) 으로 진행:

\[ \log \left[ \frac{P_{ijc}}{1 - P_{ijc}} \right] = z_{ijc} = \frac{\gamma_c - (x_{ij}^\top \beta + u_{ij}^\top \alpha_c + z_{ij}^\top T \theta_i)}{\exp(w_{ij}^\top \tau)} \quad (c = 1, \ldots, C-1) \tag{10.15} \]

특수 경우:

Proportional odds: $\alpha_c = 0, \tau = 0$.
Partial proportional odds: $\alpha_c \neq 0, \tau = 0$.
Location-scale: $\alpha_c \neq 0, \tau \neq 0$.

직관 — 한 식에 모든 모형 family

식 (10.15) 의 가치: 하나의 식이 ordinal mixed-effects 의 모든 변형을 포함.

추정 알고리즘을 일반 모형으로 짜두면:

사용자가 $\alpha_c = 0$ 으로 제약 → proportional odds 모형 적합.
사용자가 $\tau = 0$ 으로 제약 → partial PO 또는 standard 모형 적합.
사용자가 자유 추정 → location-scale 모형 적합.

→ 하나의 코드로 다양한 모형 적합. Hedeker MIXOR 같은 소프트웨어의 설계 철학.

본 § 10.2.4 의 score 와 Fisher scoring 도 식 (10.15) 에 대해 도출 — 특수 경우는 자동 처리.

2.3 식 (10.16) — Conditional Likelihood (Multinomial)

식 (10.16) — 환자별 조건부 우도

같은 환자 $i$ 의 $n_i$ 개 시점 응답을 결합:

\[ \ell(Y_i \mid \theta) = \prod_{j=1}^{n_i} \prod_{c=1}^{C} (P_{ijc} - P_{ij, c-1})^{y_{ijc}} \tag{10.16} \]

표기:

$Y_{ij}$: 환자 $i$ 의 시점 $j$ 에서의 응답 ($\in \{1, 2, \ldots, C\}$).
$y_{ijc}$: indicator. $Y_{ij} = c$ 면 1, 아니면 0.
즉 각 시점에서 $C$ 개 indicator 중 정확히 하나만 1.

직관 — Multinomial 의 곱 형태

$y_{ijc}$ 의 정의 때문에 각 시점에서 $C$ 개 항 중 하나만 살아남음 (지수가 1):

$Y_{ij} = c^*$ 라면: 우도 contribution = $P_{ijc^*} - P_{ij, c^*-1}$ (cell 확률).
다른 $c \neq c^*$ 항: $(\cdot)^0 = 1$ (생략).

따라서 식 (10.16) 은 사실상:

\[ \ell(Y_i \mid \theta) = \prod_{j=1}^{n_i} P(Y_{ij} = Y_{ij}^{\text{obs}} \mid \theta) \]

— “각 시점의 관측 범주 확률의 곱”.

§ 9.6 와의 비교:

§ 9.6 식 (9.39) 의 Bernoulli 곱:

\[ \ell(Y_i \mid \theta) = \prod_{j=1}^{n_i} \Psi(z_{ij})^{Y_{ij}} [1 - \Psi(z_{ij})]^{1-Y_{ij}} \]

이것이 식 (10.16) 의 $C = 2$ 특수 경우. Bernoulli 의 두 항 ($p$, $1-p$) 이 ordinal 의 $C$ 항으로 일반화.

Conditional independence 가정 그대로: 같은 환자의 시점이 $\theta$ 조건부로 독립. § 9.6 의 가정과 동일.

2.4 식 (10.17) — Censoring 확장 (Survival)

식 (10.17) — Right-censoring 처리

§ 10.2.3 의 survival 응용에서, 일부 응답이 censored:

\[ \ell(Y_i \mid \theta) = \prod_{j=1}^{n_i} \prod_{c=1}^{C} \left[ (P_{ijc} - P_{ij, c-1})^{d_{ij}} (1 - P_{ijc})^{1 - d_{ij}} \right]^{y_{ijc}} \tag{10.17} \]

표기:

$d_{ij} = 1$: 시점 $Y_{ij}$ 에서 사건 발생 (event).
$d_{ij} = 0$: 시점 $Y_{ij}$ 에서 censoring.

추가 절단점: $\gamma_{C+1} = \infty$, $P_{ij, C+1} = 1$. 이로써 $c = 1, \ldots, C$ 자유 추정 가능 (censoring 영역 포함).

직관 — Event 와 censoring 의 우도 contribution

$Y_{ij} = c^*$ 인 환자의 contribution (식 10.17 의 단일 항으로 간단화):

$d_{ij} = 1$ (event at $c^*$): contribution = $P_{ij c^*} - P_{ij, c^*-1}$ — cell 확률 (식 10.16 와 같음).
$d_{ij} = 0$ (censored at $c^*$): contribution = $1 - P_{ij c^*}$ — survival 확률 (“$c^*$ 너머에 사건이 일어날 확률” 의 표현).

해석:

Event 환자: 정확히 그 시점에 사건 발생 → cell 확률.
Censored 환자: 그 시점까지 관측 후 종료 → 그 시점 너머 (어디든) 사건 발생할 가능성 모두 합산 = survival 확률.

식 (10.17) 은 식 (10.16) 의 일반화:

모든 환자가 $d_{ij} = 1$ (사건 발생) 이면 식 (10.16) 으로 reduce.
Censoring 이 있는 종단 분석 (대부분의 임상 RCT) 에서 식 (10.17) 사용.

이것이 discrete-time PH mixed-effects 모형의 정확한 우도. § 10.2.3 의 framework 가 추정 단계에서 자연스럽게 통합됨.

2.5 식 (10.18) — Marginal Likelihood

Marginal 우도 — 적분으로 random effects 제거

§ 9.6 의 식 (9.40) 와 동일한 형태:

\[ h(Y_i) = \int_\theta \ell(Y_i \mid \theta) g(\theta) \, d\theta \tag{10.18} \]

$g(\theta)$ = 다변량 표준 정규 ($\mathcal{N}_r(0, I)$).

표본 전체의 marginal log-likelihood:

\[ \log L = \sum_{i=1}^N \log h(Y_i) \]

→ Marginal MLE 의 표준 절차. § 9.6 와 동일.

직관 — § 9.6 와 정확히 같은 적분

§ 9.6 의 식 (9.40) 와 식 (10.18) 의 차이는 conditional likelihood $\ell$ 의 형태뿐:

§ 9.6: Bernoulli 곱 (식 9.39).
§ 10.2.4: Multinomial 곱 (식 10.16) 또는 censoring 포함 (식 10.17).

적분 자체는 같은 차원 ($\theta$ 의 random effects 차원), 같은 분포 ($g(\theta)$ = standard normal). → 같은 Gauss-Hermite quadrature 적용.

식 (10.18) 의 적분이 닫힌 형태로 안 풀리는 이유 — § 9.6 와 동일. Logistic cdf 의 비선형성 + 곱의 적분.

핵심 메시지: ordinal mixed-effects 의 적분 추정이 이항보다 본질적으로 어렵지 않다. 단지 conditional likelihood 의 항이 더 많을 뿐.

3 Score 의 일반 형태와 미분

3.1 식 (10.19) — Score 의 일반 형태

모수 vector 와 score

모수 vector:

\[ \eta^\top = [\gamma_c : \alpha_c : \beta : \tau : T] \]

5 개 그룹: 절단점, partial PO 계수, 비례 계수, scaling 계수, Cholesky factor.

Score:

\[ \frac{\partial \log L}{\partial \eta} = \sum_{i=1}^N [h(Y_i)]^{-1} \int_\theta \frac{\partial \ell_i}{\partial \eta} \, g(\theta) \, d\theta \tag{10.19} \]

여기서 $\ell_i = \ell(Y_i \mid \theta) = \prod_j \prod_c (p_{ijc})^{y_{ijc}}$ (식 10.20), $p_{ijc} = P_{ijc} - P_{ij,c-1}$ (식 10.21).

직관 — 식 (10.19) 의 구조

식 (10.19) 의 형태는 § 9.6 식 (9.43) 과 정확히 같음:

환자별 합 + 사후 분포 weighting.
$h^{-1}(Y_i)$ 로 정규화 ($h$ 가 marginal likelihood).
$\int (\partial \ell_i / \partial \eta) g(\theta) d\theta$ 가 사후 기댓값의 적분.

새로 계산할 것은 $\partial \ell_i / \partial \eta$ — conditional likelihood 의 미분. Chain rule 로 단계별 도출.

Chain rule:

\[ \frac{\partial \ell_i}{\partial \eta} = \sum_{j, c} y_{ijc} \frac{\partial \log p_{ijc}}{\partial \eta} \cdot \ell_i \]

→ 핵심은 $\partial \log p_{ijc} / \partial \eta$. 이것이 모수 종류에 따라 다름.

3.2 미분의 두 종류 — c-무관 vs c-가변

$p_{ijc}$ 의 분해

$p_{ijc} = P_{ijc} - P_{ij, c-1} = \Psi(z_{ijc}) - \Psi(z_{ij, c-1})$.

미분:

\[ \frac{\partial p_{ijc}}{\partial \eta} = \frac{\partial \Psi(z_{ijc})}{\partial \eta} - \frac{\partial \Psi(z_{ij, c-1})}{\partial \eta} \]

Chain rule:

\[ \frac{\partial \Psi(z_{ijc})}{\partial \eta} = \psi(z_{ijc}) \cdot \frac{\partial z_{ijc}}{\partial \eta} \]

여기서 $\psi$ 는 logistic pdf.

c-무관 모수의 미분 ($\beta$, $T$)

$\beta$ 와 $T$ 는 모든 cumulative logit 에 같은 값으로 들어감 — $z_{ijc}$ 의 미분이 $c$ 에 따라 변하지 않음.

\[ \frac{\partial \log p_{ijc}}{\partial \eta} = \frac{\psi(z_{ijc}) - \psi(z_{ij, c-1})}{\Psi(z_{ijc}) - \Psi(z_{ij, c-1})} \cdot \frac{\partial z_{ijc}}{\partial \eta} \]

$z_{ijc}$ 의 모수별 미분:

\[ \frac{\partial z_{ijc}}{\partial \beta} = \frac{-x_{ij}}{\exp(w_{ij}^\top \tau)} \]

\[ \frac{\partial z_{ijc}}{\partial v(T)} = \frac{-J_r (\theta \otimes z_{ij})}{\exp(w_{ij}^\top \tau)} \]

여기서 $v(T)$ 는 Cholesky factor 의 자유 모수 vector ($r(r+1)/2$ 개), $J_r$ 는 elimination matrix (§ 9.6.2 의 식 9.54 와 동일).

가중치의 직관적 해석: $\frac{\psi(z_{ijc}) - \psi(z_{ij, c-1})}{\Psi(z_{ijc}) - \Psi(z_{ij, c-1})}$ 는 “두 절단점 사이 cell 의 평균 pdf 기울기 / cell 확률” — 잠재 변수 척도에서 $z$ 변화에 대한 cell 확률 변화의 비율.

이항 ($C = 2$) 에서는 $\Psi(z_{ij,1}) = \Psi(z_{ij}), \Psi(z_{ij,0}) = 0$ 이라 가중치가 단순화 → § 9.6 식 (9.44) 의 잔차 형태로 reduce. § 9.6 의 자연 일반화.

c-가변 모수의 미분 ($\gamma_c$, $\alpha_c$)

$\gamma_c$ 는 cumulative logit $c$ 에만 등장 — 다른 logit 의 미분은 0. Kronecker delta 로 표현:

\[ \frac{\partial \log p_{ijc}}{\partial \eta_{c'}} = \frac{\delta_{c c'} \psi(z_{ijc}) \cdot \frac{\partial z_{ijc}}{\partial \eta_{c'}} - \delta_{c-1, c'} \psi(z_{ij, c-1}) \cdot \frac{\partial z_{ij, c-1}}{\partial \eta_{c'}}}{\Psi(z_{ijc}) - \Psi(z_{ij, c-1})} \]

여기서 $\delta_{cc'} = 1$ if $c = c'$, else 0.

$\eta_{c'}$ 의 미분:

\[ \frac{\partial z_{ijc}}{\partial \gamma_{c'}} = \frac{1}{\exp(w_{ij}^\top \tau)} \quad (\text{c-가변}, \delta \text{ 만 살림}) \]

\[ \frac{\partial z_{ijc}}{\partial \alpha_{c'}} = \frac{-u_{ij}}{\exp(w_{ij}^\top \tau)} \]

Kronecker delta 의 의미: 각 cell $c$ 의 확률은 두 절단점 $\gamma_{c-1}, \gamma_c$ 만 관여. 다른 절단점 $\gamma_{c'}$ ($c' \neq c-1, c$) 에는 의존 안 함. delta 가 이 구조를 자동 표현.

도식:

범주 1: 절단점 γ_1 만 관여 (γ_0 = -∞)
범주 2: 절단점 γ_1, γ_2
범주 3: 절단점 γ_2, γ_3
...
범주 C: 절단점 γ_{C-1} (γ_C = ∞)

이 구조 때문에 score 행렬이 띠 형태 (band) — 인접 절단점에만 비대각 원소.

3.3 Scaling $\tau$ 의 미분

$\tau$ 미분의 chain rule

$z_{ijc} = \frac{\gamma_c - (\cdots)}{\exp(w^\top \tau)}$ 의 분모 $\exp(w^\top \tau)$ 가 $\tau$ 에 의존.

$\partial \exp(w^\top \tau) / \partial \tau = \exp(w^\top \tau) \cdot w$.

Quotient rule:

\[ \frac{\partial z_{ijc}}{\partial \tau} = \frac{-[\exp(w^\top \tau) w_{ij}] [\gamma_c - (\cdots)]}{[\exp(w^\top \tau)]^2} = -w_{ij}^\top z_{ijc} \]

(분자의 $[\gamma_c - (\cdots)] / \exp(w^\top \tau) = z_{ijc}$ 사용.)

결과:

\[ \frac{\partial \log p_{ijc}}{\partial \tau} = \frac{\psi(z_{ij, c-1}) z_{ij, c-1} - \psi(z_{ijc}) z_{ijc}}{\Psi(z_{ijc}) - \Psi(z_{ij, c-1})} \cdot w_{ij} \]

직관 — Scaling 미분의 형태

$\beta, T$ 의 미분과 비교하면:

$\beta$ 미분 가중치: $[\psi(z_{ijc}) - \psi(z_{ij, c-1})] / \Delta P$.
$\tau$ 미분 가중치: $[\psi(z_{ij, c-1}) z_{ij, c-1} - \psi(z_{ijc}) z_{ijc}] / \Delta P$.

차이: $\tau$ 미분에는 $z_{ij,c-1}, z_{ijc}$ 의 절댓값이 곱해짐 — scaling 이 잠재 변수의 polynomial 효과를 만들기 때문.

이 추가 한 단계 ($\times z$) 가 $\tau$ 추정의 수치적 어려움을 약간 늘림 — outlier 환자에 더 민감.

실무: scaling 모수는 $\beta, T$ 보다 추정 안정성이 약간 낮음. 표본 크기가 부족하면 수렴 어려움 가능. Bayesian 우회가 자주 권고.

4 Fisher Scoring — 식 (10.22)

식 (10.22) — Update 알고리즘

§ 9.6 의 식 (9.47) 와 동일한 형태:

\[ \Theta_{l+1} = \Theta_l + I(\Theta_l)^{-1} \frac{\partial \log L}{\partial \Theta_l} \tag{10.22} \]

정보 행렬:

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

→ Score 의 outer product 합. 양정치 보장.

직관 — § 9.6 와 같은 알고리즘

Ordinal 추정 알고리즘이 이항과 거의 동일:

초기값 설정: 단순 ordinal logit (single-level) 으로 $\beta, \gamma$ 추정. Random effects 분산 = 1 또는 small value.
각 iteration:
- 현재 모수에서 cell 확률 $p_{ijc}$ 계산.
- Conditional likelihood $\ell_i$ 계산 (식 10.16 또는 10.17).
- Marginal likelihood $h_i$ 계산 (식 10.18, Gauss-Hermite quadrature).
- Score $\partial \log L / \partial \eta$ 계산 (식 10.19, 위의 미분 공식).
- Information matrix $I(\Theta)$ 계산.
- 식 (10.22) 로 update.
수렴 검사: $|_{l+1} - l| < $ 임계값 또는 $|L{l+1} - L_l| < $ 임계값.

수렴 후:

표준오차: $\text{SE}(\widehat\Theta) = \sqrt{\text{diag}(I^{-1})}$.
Empirical Bayes random effects: 식 (9.48) 의 ordinal 버전 (cell 확률로 변경).
LR test: $-2 \log L$ 차이의 $\chi^2$.

→ § 9.6 의 toolkit 그대로 적용. 새로 짤 알고리즘 없음.

4.1 Gauss-Hermite Quadrature 의 적용

적분의 합 근사

§ 9.6.3 의 식 (9.55-9.59) 와 동일하게 Gauss-Hermite quadrature 로 marginal likelihood 적분 근사:

\[ h(Y_i) \approx \sum_{q=1}^{Q^r} \ell(Y_i \mid B_q) A(B_q) \]

여기서 $B_q$ 는 quadrature 노드, $A(B_q)$ 는 가중치, $r$ 은 random effects 차원.

각 노드 $B_q$ 에서:

\[ z_{ijcq} = \frac{\gamma_c - (x_{ij}^\top \beta + u_{ij}^\top \alpha_c + z_{ij}^\top T B_q)}{\exp(w_{ij}^\top \tau)} \]

$\theta$ 자리에 $B_q$ 대입.

Conditional cell probability:

\[ P_{ijcq} = \Psi(z_{ijcq}), \quad p_{ijcq} = P_{ijcq} - P_{ij, c-1, q} \]

Conditional likelihood:

\[ \ell(Y_i \mid B_q) = \prod_j \prod_c (p_{ijcq})^{y_{ijc}} \]

차원 폭발과 Adaptive Quadrature

§ 9.6.3 의 caveat — $Q^r$ 폭발 — 이 ordinal 에도 적용:

$r$ (random effects)	$Q = 5$	$Q = 10$	$Q = 20$
1	5	10	20
2	25	100	400
3	125	1000	8000

차원 ≥ 3 이면 adaptive quadrature 권고 (§ 9.6.3 의 식 9.60-9.61).

Ordinal 의 추가 부담:

각 노드에서 $C-1$ 개 cumulative 와 $C$ 개 cell 확률 계산.
미분도 $C-1$ 개 절단점에 대해 별도.
→ 같은 random effects 차원에서 이항보다 약 $C$ 배 더 많은 계산.

실무에서는 6-7 범주 이상 ordinal 모형의 mixed-effects 적합이 매우 느림 — Bayesian (Stan, brms) 가 더 효율적일 수 있음.

5 응용 분야

분야	추정 권고
임상 ordinal RCT (3-5 범주)	R `ordinal::clmm` (Laplace) 또는 SAS `PROC NLMIXED` (quadrature)
Mixed PO + 시간 가변 공변량	brms (Bayesian)
Partial PO mixed-effects	Hedeker MIXOR 또는 brms `cs()`
Location-scale mixed-effects	SAS `PROC NLMIXED` (직접 코딩) 또는 brms (sigma 모형화)
6+ 범주 + multiple random effects	brms / Stan 또는 INLA
Discrete-time PH	R `glmer cloglog` (pooled dichotomous)

6 코드 예시

6.1 Step 1: Cell 확률과 Conditional Likelihood (numpy)

import numpy as np
from scipy.special import expit


def cell_probs(z_cuts: np.ndarray) -> np.ndarray:
    """식 (10.14): cumulative 차이로 cell 확률.

    z_cuts: shape (n, C-1) 의 z_ijc 행렬.
    Returns: shape (n, C) 의 cell 확률.
    """
    n, C_minus_1 = z_cuts.shape
    P_cum = expit(z_cuts)  # cumulative 확률
    # P_{ij,0} = 0, P_{ij,C} = 1
    P_full = np.column_stack([np.zeros(n), P_cum, np.ones(n)])
    # cell 확률 = 인접 차이
    cell = P_full[:, 1:] - P_full[:, :-1]
    return cell


# 단순 사례: 4 범주 + single 환자
gamma = np.array([-1.5, 0, 1.5])
beta = 0.5
x = np.array([0, 1])  # 두 환자
z = gamma - x[:, None] * beta  # shape (2, 3)

cell = cell_probs(z)
print("Cell 확률:")
print(f"  x=0: {cell[0].round(3)}")
print(f"  x=1: {cell[1].round(3)}")
print(f"  합 검증: {cell.sum(axis=1)}")

검증 포인트

각 환자의 cell 확률 합이 1.
$x = 1$ 그룹은 양의 $\beta$ 효과로 더 높은 범주에 확률 이동.

이 함수가 식 (10.14) 의 직접 구현. 모든 ordinal 우도 계산의 토대.

6.2 Step 2: Conditional Likelihood with Censoring (식 10.17)

def conditional_likelihood_with_censoring(z_cuts: np.ndarray,
                                            Y: np.ndarray, d: np.ndarray) -> np.ndarray:
    """식 (10.17): censoring 포함 conditional likelihood (single observation).

    z_cuts: shape (n, C-1).
    Y: shape (n,) 응답 (1, ..., C).
    d: shape (n,) censoring indicator (1=event, 0=censored).
    """
    n, C_minus_1 = z_cuts.shape
    C = C_minus_1 + 1
    P_cum = expit(z_cuts)
    P_full = np.column_stack([np.zeros(n), P_cum, np.ones(n)])

    likelihood = np.zeros(n)
    for i in range(n):
        c = int(Y[i]) - 1  # 0-indexed
        if d[i] == 1:
            # Event: cell probability
            likelihood[i] = P_full[i, c+1] - P_full[i, c]
        else:
            # Censored: 1 - P_{ij, c}
            likelihood[i] = 1 - P_full[i, c+1]
    return likelihood


# 시뮬레이션 — 4 시점 (= 4 범주), 일부 censored
np.random.seed(2026)
n = 100
z_cuts_demo = np.random.normal(0, 1, size=(n, 3))
Y_demo = np.random.choice([1, 2, 3, 4], size=n)
d_demo = np.random.choice([0, 1], size=n, p=[0.3, 0.7])

L = conditional_likelihood_with_censoring(z_cuts_demo, Y_demo, d_demo)
print(f"평균 conditional likelihood: {L.mean():.4f}")
print(f"Censored 환자 (d=0) 예시:")
mask = (d_demo == 0)
print(f"  Y={Y_demo[mask][:5]}, L={L[mask][:5].round(4)}")

Censoring 처리의 검증

Event 환자: cell probability = $P_{Y} - P_{Y-1}$.
Censored 환자: survival probability = $1 - P_{Y}$.
Censored 의 $L$ 가 보통 더 큼 (한 cell 이 아니라 모든 미래 가능성을 합).

이 함수가 식 (10.17) 의 직접 구현. discrete-time survival mixed-effects 의 토대.

6.3 Step 3: Marginal Likelihood + Gauss-Hermite (numpy)

import numpy as np
from numpy.polynomial.hermite_e import hermegauss
from scipy.special import expit


def marginal_likelihood_ordinal(Y_i: np.ndarray, x_i: np.ndarray,
                                 beta: np.ndarray, gamma: np.ndarray,
                                 sigma_v: float, Q: int = 20) -> float:
    """단일 환자의 marginal likelihood (random intercept ordinal mixed).

    식 (10.18) 를 Gauss-Hermite quadrature 로 근사.
    Y_i: 시점별 응답 (n_i,) — 1, ..., C.
    x_i: 시점별 공변량 (n_i, p).
    """
    nodes, weights = hermegauss(Q)
    weights = weights / np.sqrt(2 * np.pi)

    h = 0.0
    for q in range(Q):
        # 각 노드에서 z_ijc 계산
        # x β + sigma_v θ — random intercept만, x_i β 가 모든 시점 공통 절편
        eta_i = x_i @ beta + sigma_v * nodes[q]  # shape (n_i,)
        # cumulative cuts: z_ijc = gamma_c - eta_ij
        z_cuts = gamma[None, :] - eta_i[:, None]  # shape (n_i, C-1)

        # cell 확률
        n_i, C_minus_1 = z_cuts.shape
        P_cum = expit(z_cuts)
        P_full = np.column_stack([np.zeros(n_i), P_cum, np.ones(n_i)])
        cell = P_full[:, 1:] - P_full[:, :-1]  # (n_i, C)

        # 관측 cell 의 확률
        log_lik = 0.0
        for j in range(n_i):
            c = int(Y_i[j]) - 1
            log_lik += np.log(cell[j, c] + 1e-300)

        h += weights[q] * np.exp(log_lik)
    return h


# 단일 환자 사례
Y_demo = np.array([2, 3, 3, 4, 4])  # 5 시점, 4 범주
x_demo = np.column_stack([np.ones(5), np.arange(5) * 0.3])  # 절편 + 시간
beta_demo = np.array([0.0, 0.5])
gamma_demo = np.array([-1.5, 0, 1.5])
sigma_v_demo = 0.8

h = marginal_likelihood_ordinal(Y_demo, x_demo, beta_demo, gamma_demo, sigma_v_demo)
print(f"Marginal likelihood h(Y_i) = {h:.6e}")
print(f"Log marginal = {np.log(h):.4f}")

적분 검증

이 코드가 § 9.6 식 (9.58) 의 ordinal 버전. 같은 Gauss-Hermite 노드와 가중치 사용.

수렴 검증: $Q$ 를 5, 10, 20, 40 으로 늘려가며 결과 안정성 확인. $Q \geq 10$ 이면 보통 충분.

이 함수를 표본 전체 $\sum_i \log h(Y_i)$ 로 모아 marginal log-likelihood. 그 결과를 모수 ($\beta, \gamma, \sigma_v$) 의 함수로 보고 최대화.

6.4 Step 4: 실무 적합 (R `ordinal::clmm`)

library(ordinal)


# 종단 ordinal 데이터 시뮬레이션 (random intercept)
set.seed(2026)
n_subjects <- 200
n_times <- 5
df <- expand.grid(subject = 1:n_subjects, time = 0:(n_times-1))
df$drug <- rep(rbinom(n_subjects, 1, 0.5), each = n_times)
upsilon <- rep(rnorm(n_subjects, 0, 1.0), each = n_times)

# 잠재 변수 + 4 범주
y_latent <- 0.5 - 0.4 * df$time + 0.3 * df$drug + upsilon + rlogis(nrow(df))
gamma <- c(-2, -0.5, 1.5)
df$y <- cut(y_latent, breaks = c(-Inf, gamma, Inf),
            labels = 1:4, ordered = TRUE)

# Random intercept proportional odds
fit <- clmm(y ~ time + drug + (1 | subject), data = df, link = "logit")
summary(fit)

# 추정 모수
cat("\n추정값:\n")
cat("  Threshold (gamma): ", round(fit$alpha, 3), "\n")  # 절단점
cat("  Coefficients (beta): ", round(fit$beta, 3), "\n")  # 회귀
cat("  Random intercept SD: ", round(sqrt(VarCorr(fit)$subject[1,1]), 3), "\n")

# ICC
sigma_v_sq <- VarCorr(fit)$subject[1, 1]
icc <- sigma_v_sq / (sigma_v_sq + pi^2 / 3)
cat("  ICC: ", round(icc, 3), "\n")

R clmm 의 추정 알고리즘

ordinal::clmm 의 기본 추정:

Laplace approximation (Naïve Gauss-Hermite 보다 빠름).
nAGQ > 1 옵션으로 adaptive Gauss-Hermite 가능 (단 random intercept 만).
다중 random effects 는 Laplace 만.

대안 패키지:

glmmTMB: AD + Laplace, ordinal 가능.
brms: Bayesian (Stan), 가장 유연.
SAS PROC NLMIXED: adaptive quadrature, 식 (10.15) 의 모든 변형 지원.

일반 권고:

빠른 prototype: clmm.
정확한 결과: SAS NLMIXED (frequentist) 또는 brms (Bayesian).
Partial PO 또는 location-scale: SAS 또는 brms.

7 관련 주제

선행 지식

Ch.10 Overview — Cumulative ordinal framework
§ 10.2 ~ 10.2.1 — Proportional + partial proportional odds
§ 10.2.2 — Location-scale 모형 (식 10.15 의 통합)
§ 10.2.3 — Discrete-time survival (censoring 식 10.17 의 응용)
§ 9.6 추정 — 이항 GLMM 의 marginal MLE (toolkit 의 토대)

후속 주제 (Ch.10 sub-posts)

§ 10.3 — NIMH 4 범주 ordinal 분석 (Ch.9 이항 분석과 비교)
§ 10.4 — 노숙자 보건서비스 (section 8 certificate, partial PO 응용)
§ 10.5 — Ch.10 Summary

관련 개념

Hedeker, Siddiqui & Hu (2000) — Censored ordinal mixed-effects MLE 도출
Bock & Aitkin (1981) — Marginal MLE + Gauss-Hermite quadrature 원전
Magnus (1988) — 행렬 미분 toolkit (elimination matrix)
Pinheiro & Bates (1995) — Adaptive quadrature
§ 9.6.2 Cholesky 미분 — 식 (9.54) 의 ordinal 적용
Ch.11 GLMM 명목 — 더 큰 모수 차원의 추정 (범주별 별도 회귀)

§ 10.2.4 — Cumulative Ordinal Mixed-Effects 추정: Marginal MLE 와 Fisher Scoring

1 들어가며 — § 9.6 의 직접 확장

2 § 10.2.4 — Cumulative Ordinal 의 Marginal MLE

2.1 식 (10.14) — Cell 확률

2.2 식 (10.15) — General \(z_{ijc}\)

2.3 식 (10.16) — Conditional Likelihood (Multinomial)

2.4 식 (10.17) — Censoring 확장 (Survival)

2.5 식 (10.18) — Marginal Likelihood

3 Score 의 일반 형태와 미분

3.1 식 (10.19) — Score 의 일반 형태

3.2 미분의 두 종류 — c-무관 vs c-가변

3.3 Scaling \(\tau\) 의 미분

4 Fisher Scoring — 식 (10.22)

4.1 Gauss-Hermite Quadrature 의 적용

5 응용 분야

6 코드 예시

6.1 Step 1: Cell 확률과 Conditional Likelihood (numpy)

6.2 Step 2: Conditional Likelihood with Censoring (식 10.17)

6.3 Step 3: Marginal Likelihood + Gauss-Hermite (numpy)

6.4 Step 4: 실무 적합 (R `ordinal::clmm`)

7 관련 주제

\(r\) (random effects)	\(Q = 5\)	\(Q = 10\)	\(Q = 20\)
1	5	10	20
2	25	100	400
3	125	1000	8000

1 들어가며 — § 9.6 의 직접 확장

2 § 10.2.4 — Cumulative Ordinal 의 Marginal MLE

2.1 식 (10.14) — Cell 확률

2.2 식 (10.15) — General \(z_{ijc}\)

2.3 식 (10.16) — Conditional Likelihood (Multinomial)

2.4 식 (10.17) — Censoring 확장 (Survival)

2.5 식 (10.18) — Marginal Likelihood

3 Score 의 일반 형태와 미분

3.1 식 (10.19) — Score 의 일반 형태

3.2 미분의 두 종류 — c-무관 vs c-가변

3.3 Scaling \(\tau\) 의 미분

4 Fisher Scoring — 식 (10.22)

4.1 Gauss-Hermite Quadrature 의 적용

5 응용 분야

6 코드 예시

6.1 Step 1: Cell 확률과 Conditional Likelihood (numpy)

6.2 Step 2: Conditional Likelihood with Censoring (식 10.17)

6.3 Step 3: Marginal Likelihood + Gauss-Hermite (numpy)

6.4 Step 4: 실무 적합 (R ordinal::clmm)

7 관련 주제

6.4 Step 4: 실무 적합 (R `ordinal::clmm`)