1 들어가며 — Baseline 의 한계를 mixed-effects 로 깨기

§ 9.7~9.7.1 에서 NIMH 정신분열증 데이터에 fixed-effects logistic 을 적용한 결과, Drug × Time 효과가 통계적으로 비유의 (p = .11) — 그러나 표 9.2 의 OR 패턴 (1.54 → 0.49 → 0.25 → 0.29) 은 명확한 약효를 시사. 이 모순의 원인은 종단 데이터에서 같은 환자의 시점이 독립이라는 가정이 깨졌기 때문.

§ 9.7.2 는 환자 절편을 랜덤화 (random intercept) 로, § 9.7.3 는 환자 절편과 기울기를 모두 랜덤화 (random intercept + trend) 로 확장한다. 같은 데이터, 같은 변수, 다른 모형 — 결론이 어떻게 달라지는지 직접 비교.

한 줄 요약

“§ 9.7.2 = random intercept 추가만으로 Drug × Time 가 -0.418 → -1.015 로 더 큰 효과 + p = .11 → p < .001 로 유의화. 환자 이질성 (\(\widehat\sigma_v = 2.12\), ICC = 0.58) 이 매우 큼. § 9.7.3 = random trend 까지 추가, Drug × Time 효과 -1.587 로 더욱 강함. 절편-기울기 음의 상관 (-0.47) 으로 ‘baseline 가 심한 환자가 더 빠르게 호전’ 패턴. 환자 EB trajectory 가 그룹 평균을 넘어 광범위한 개인차를 시연. 결론: 모형 선택이 RCT 의 약효 결론을 바꾼다.”

2 § 9.7.2 — Random Intercept Logistic Regression

2.1 모형 정의 — 식 (9.63-9.64)

Level-1 (within-subjects) 모형

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = b_{0i} + b_{1i} \sqrt{\text{Week}}_j \tag{9.63} \]

피험자 \(i\) 의 시점 \(j\) 에서 logit 이 환자별 절편 \(b_{0i}\) 와 환자별 기울기 \(b_{1i}\) 의 1 차식.

Level-2 (between-subjects) 모형 — random intercept만

\[ b_{0i} = \beta_0 + \beta_2 \text{Drug}_i + \upsilon_{0i} \] \[ b_{1i} = \beta_1 + \beta_3 \text{Drug}_i \tag{9.64} \]

랜덤 효과 \(\upsilon_{0i} \sim \mathcal{N}(0, \sigma_v^2)\) — 환자별 절편의 변동만 허용. 기울기는 그룹 평균으로 고정.

직관 — 각 모수의 임상적 의미

식 (9.64) 의 4 개 고정효과:

\(\beta_0\): placebo (\(\text{Drug} = 0\)) 환자의 week 0 logit. 임상적으로 “처치 안 받은 평균 환자의 baseline 중증도”.
\(\beta_1\): placebo 환자의 sqrt(week) 1 단위당 logit 변화. “자연 회복 + 입원 효과” 의 표현.
\(\beta_2\): drug (\(\text{Drug} = 1\)) 환자의 week 0 logit 차이. RCT 가 잘 됐으면 0 에 가까움.
\(\beta_3\): drug 환자의 sqrt(week) 1 단위당 logit 변화의 추가량. 임상적 핵심 — 약효의 차등 효과.

\(\upsilon_{0i}\) 의 역할:

같은 그룹 내에서도 환자마다 baseline 중증도가 다름.
\(\upsilon_{0i} > 0\) 인 환자는 항상 그룹 평균보다 더 심함.
\(\upsilon_{0i} < 0\) 인 환자는 항상 그룹 평균보다 덜 심함.
이 절편 차이가 시간에 따라 일정하게 유지 (절편만 random, 기울기는 fixed).

§ 9.7.3 에서는 기울기까지 랜덤화 — 호전 속도의 환자 차이도 모형화.

Subject-specific 해석의 강조 — 모든 회귀 계수가 “같은 \(\upsilon_{0i}\) 수준의 환자 sub-population” 의 효과. § 9.5 의 식 (9.16) 비축소성에 의해 fixed-effects 와 척도가 다름.

2.2 추정 결과 — 표 9.4

표 9.4 — Random Intercept Logistic Regression 결과 (N = 437)

모수	Estimate	SE	\(Z\)	\(p <\)
Intercept	5.387	0.535	10.07	.001
Drug (0=placebo, 1=drug)	-0.025	0.601	-0.04	.97
Time (\(\sqrt{\text{week}}\))	-1.500	0.228	-6.59	.001
Drug × Time	-1.015	0.274	-3.70	.001
Intercept SD (\(\widehat\sigma_v\))	2.116	0.215

\(-2 \log L = 1249.74\)

직관 — Fixed-Effects 와 Random Intercept 의 비교

모수	Fixed (Table 9.3)	Random Intercept (Table 9.4)	변화
Intercept	3.703	5.387	+1.684 (45 % 증가)
Drug	-0.405	-0.025	거의 0 으로
Time	-1.113	-1.500	-0.387 (35 % 더 큼)
Drug × Time	-0.418 (p=.11)	-1.015 (p<.001)	-0.597 (2.4 배), 유의

관찰 1 — 회귀 계수가 전반적으로 커짐: 식 (9.16) 의 비축소성 때문 (subject-specific scale 이 marginal scale 보다 큼). 모든 효과가 잠재 분산 \(\sigma_v^2 + \pi^2/3\) 으로 scale.

관찰 2 — Drug 주효과가 거의 0: -0.405 → -0.025. Random intercept 가 baseline 차이를 흡수 — 환자 이질성으로 충분히 설명되어 그룹 평균 차이가 사라짐. RCT randomization 이 잘 됐다는 강한 증거.

관찰 3 — Drug × Time 가 유의화: 가장 중요한 변화. Fixed 에서 p = .11 (비유의) → Random intercept 에서 p < .001 (강한 유의). 약효의 차등 효과가 통계적으로 입증됨.

관찰 4 — Intercept SD = 2.12: 환자별 절편의 표준편차가 매우 큼. 그룹 평균 절편이 5.39 인데 환자 표준편차가 2.12 → 환자 사이의 baseline 차이가 약 ±2.12 (logit 척도). odds 척도로는 약 8 배 차이.

2.3 LR 검정 — Random Intercept vs Fixed-Effects

우도비 검정

\(H_0\): \(\sigma_v^2 = 0\) (랜덤 효과 불필요).

\[ \chi^2_1 = -2 \log L_{\text{fixed}} - (-2 \log L_{\text{random}}) = 1362.06 - 1249.74 = 112.3 \]

자유도 1 (\(\sigma_v^2\) 한 개), p < .001.

직관 — “환자 이질성을 무시할 수 없다” 의 강한 증거

\(\chi^2_1 = 112.3\) 는 매우 큰 검정 통계량 (자유도 1 의 .001 quantile 은 약 10.83). p-value 가 사실상 0.

해석:

데이터가 확실히 환자 안에서 상관 — 같은 환자의 시점이 독립이 아님.
Fixed-effects 모형은 이 상관을 무시 → 부적절한 모형.
Random intercept 가 데이터 구조에 훨씬 더 잘 맞음.

LR 검정의 boundary 문제: \(\sigma_v^2 = 0\) 은 모수 공간의 경계 (음수 분산 불가능). 표준 \(\chi^2\) 분포가 약간 부정확 — 실제 분포는 \(0.5 \chi^2_0 + 0.5 \chi^2_1\) (Self & Liang 1987). 그러나 본 사례는 검정 통계량이 너무 커 boundary 보정 후에도 강한 유의.

2.4 ICC — 식 (9.65)

식 (9.65) — Logistic 모형의 ICC

\[ \widehat\rho = \frac{\widehat\sigma_v^2}{\widehat\sigma_v^2 + \pi^2/3} = \frac{2.116^2}{2.116^2 + \pi^2/3} = \frac{4.477}{4.477 + 3.290} = 0.58 \tag{9.65} \]

직관 — ICC = 0.58 의 임상적 의미

§ 9.5.1 에서 본 ICC 의 정의 — 잠재 반응 성향 분산 중 환자 간 분산의 비율. 0.58 은 큰 값:

잠재 반응 성향의 변동 중 58 % 가 환자 차이로 설명.
42 % 만이 시점 잔차 (측정 오차, 환경 변동) 로 설명.
→ 환자 이질성이 매우 강함.

비교 기준:

ICC 범위	해석	모형 권고
0 - 0.05	거의 무시 가능	Fixed-effects 충분
0.05 - 0.20	작지만 유의	Random effects 권고
0.20 - 0.50	보통	Random effects 필수
0.50 - 0.80	큼	Random effects 매우 중요
0.80+	매우 큼	데이터 구조 점검 (의사 반복 가능)

NIMH 데이터의 0.58 은 “큼” 범주 — 정신과 임상시험에서 환자 차이가 시점 변동보다 더 중요한 변동 요인. 이 정보가 단일 환자 수준의 의사결정 (개인 맞춤 약물 등) 에 시사하는 바가 크다.

2.5 ±1 SD Trend Lines — Figure 9.10

Subject-specific 적합 logit (정중도 예시)

식 (9.62) 에 random effect 추가:

\[ P(Y_{ij} = 1) = \frac{1}{1 + \exp[-(5.39 - 0.03 D_i - 1.50 \sqrt{W_j} - 1.02 D_i \sqrt{W_j} + \upsilon_{0i})]} \]

\(\upsilon_{0i} = \pm 2.12\) (즉 \(\pm 1 \widehat\sigma_v\)) 환자에 대해 그룹별로 trend line 그림.

→ 4 개 trend line: placebo (-1 SD), placebo (+1 SD), drug (-1 SD), drug (+1 SD).

직관 — 환자 이질성의 시각화

Figure 9.10 의 패턴:

Placebo (-1 SD): 처치 안 받았지만 baseline 가 낮은 환자 → week 6 에서 P(sick) < 0.5 (즉 well 이 더 많음).
Drug (+1 SD): 처치 받았지만 baseline 가 높은 환자 → week 6 에서 P(sick) 가 여전히 0.5 보다 높을 수 있음.

충격적 결과: “처치 안 받았지만 가벼운 환자가 처치 받았지만 무거운 환자보다 더 빠르게 호전”. 환자 사이의 baseline 차이가 처치 효과보다 클 수 있다는 시각적 증거.

이 발견이 의미:

평균 효과 (그룹 평균 trend) 만 보면 약효가 명확.
그러나 개인 수준에서는 처치 효과가 baseline 차이에 압도될 수 있음.
→ Personalized medicine 의 동기 — 환자 baseline 을 알면 약물 선택을 다르게 할 수 있을지.

§ 9.7.3 의 random trend 모형은 호전 속도의 환자 차이까지 모형화 → 더 풍부한 개인차 분석.

2.6 Empirical Bayes Histograms — Figure 9.11-9.12

환자별 \(\widehat\upsilon_{0i}\) 분포의 패턴

Figure 9.11 (placebo, N=108) 와 Figure 9.12 (drug, N=329) 의 EB 추정 히스토그램.

저자가 본문에서 강조한 두 통계:

Placebo 환자의 약 10 % 가 \(\widehat\upsilon_{0i} \leq -2.12\): -1 SD 이하 — 처치 안 받았지만 baseline 가 매우 낮은 (가벼운) 환자가 표본의 1 / 10. 이 환자들은 placebo 만으로도 회복 가능.
Drug 환자의 약 25 % 가 \(\widehat\upsilon_{0i} > 1.93\): +1 SD 근처 — 약물 받았지만 baseline 가 매우 높은 (심한) 환자가 표본의 1 / 4. 이 환자들은 약효가 미미할 수 있음.

→ “Number Needed to Treat (NNT)” 같은 임상 지표를 이런 분포에서 직접 읽을 수 있다. 평균 효과만 보면 놓치는 환자별 변동.

Bayes Shrinkage 효과 확인: EB 추정값의 분포가 사전 분포 (\(\mathcal{N}(0, 2.12^2)\)) 와 비슷하지만 약간 좁음. 데이터가 적은 환자 (적게 측정된 환자) 가 0 으로 수축되어 분포 폭이 줄어듦.

2.7 Marginalization 4 단계 절차

Subject-specific → Marginal 변환

표 9.4 의 회귀 계수는 subject-specific (랜덤 효과 조건부). Marginal 확률을 얻으려면:

Step 1: 선형 예측자 (랜덤 효과 제외) 계산:

\[ \widehat{y}_i = X_i \widehat\beta \]

Step 2: Marginalization 벡터 계산:

\[ \widehat{s} = \frac{1}{\sigma} \left[ \text{Diag}(\widehat{V}(y_i)) \right]^{1/2} \]

여기서 \(\widehat{V}(y_i) = Z_i \widehat\Sigma Z_i^\top + \sigma^2 I_i\), \(\sigma = 1\) (probit) 또는 \(\sigma = \pi/\sqrt{3}\) (logistic).

Random intercept 만 있으면 \(Z_i = 1_i\), 결과:

\[ \widehat{s} = \sqrt{\widehat\sigma_v^2 / \sigma^2 + 1} \]

Step 3: Element-wise 나눗셈:

\[ \widehat{z}_i = \widehat{y}_i / \widehat{s} \]

Step 4: Marginal 확률:

\[ \widehat{p}_i = \Psi(\widehat{z}_i) \quad \text{(logistic)} \]

직관 — “잠재 분산 보정으로 marginal 환산”

핵심 발상은 식 (9.16) 의 비축소성 — 잠재 분산이 클수록 같은 확률을 만들기 위한 평균값이 작아져야 함. Subject-specific 추정값을 그대로 logistic cdf 에 넣으면 잠재 분산이 \(\sigma_v^2 + \sigma_\epsilon^2\) 인 척도의 답. Marginal 답을 얻으려면 잠재 분산을 \(\sigma_\epsilon^2\) 으로 보정해야 한다.

Step 2 의 \(\widehat{s}\) 가 정확히 이 보정 인자 — \(\widehat{s} = \sqrt{(\sigma_v^2 + \sigma_\epsilon^2) / \sigma_\epsilon^2}\).

Zeger 보정: 본문 권고에 따라 logistic 의 \(\sigma_\epsilon^2 = \pi^2/3\) 대신 \((15\pi/16\sqrt{3})^2 \cdot \pi^2/3\) 사용 → 실제 logistic ↔︎ marginal 변환이 더 정확.

이 4 단계 절차는 SAS PROC IML, R, Python 어느 환경에서도 같은 방식으로 구현 가능. 본문 표 9.5 가 SAS 예시.

Figure 9.13 의 메시지: marginal 적합 곡선이 표 9.2 의 관측 비율과 매우 잘 일치 → 모형이 평균 패턴을 정확히 포착. sqrt(week) 변환의 단일 모수가 두 그룹의 시간 추세를 모두 잘 표현.

3 § 9.7.3 — Random Intercept and Trend Logistic Regression

3.1 모형 정의 — 식 (9.66-9.67)

Level-1 (within-subjects) 모형 — 식 (9.66)

식 (9.63) 와 동일:

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = b_{0i} + b_{1i} \sqrt{\text{Week}}_j \tag{9.66} \]

Level-2 (between-subjects) 모형 — 식 (9.67)

이제 절편과 기울기 모두 랜덤 효과 포함:

\[ b_{0i} = \beta_0 + \beta_2 \text{Drug}_i + \upsilon_{0i} \] \[ b_{1i} = \beta_1 + \beta_3 \text{Drug}_i + \upsilon_{1i} \tag{9.67} \]

랜덤 효과 벡터 \(\upsilon_i = (\upsilon_{0i}, \upsilon_{1i})^\top \sim \mathcal{N}_2(0, \Sigma_v)\) — 이변량 정규.

Probit 버전은 Gibbons & Bock (1987) 가 처음 제안, NIMH 데이터에 Gibbons & Hedeker (1994) 적용. 본 sub-post 는 logistic 버전.

직관 — 환자별 호전 속도까지 다르다는 가정

§ 9.7.2 (random intercept 만) 는 “모든 환자가 같은 속도로 호전, 단 baseline 만 다름” 을 가정. 그러나 임상 현실은 “어떤 환자는 빠르게 호전, 어떤 환자는 느리게 호전” — 호전 속도의 개인차.

식 (9.67) 의 \(\upsilon_{1i}\) 가 이 차이를 모형화:

\(\upsilon_{1i} > 0\): 그룹 평균보다 호전이 느림 (logit 감소가 작음).
\(\upsilon_{1i} < 0\): 그룹 평균보다 호전이 빠름.

\(\upsilon_{0i}\) 와 \(\upsilon_{1i}\) 의 상관 (\(\Sigma_v\) 의 비대각 원소) 도 임상적으로 중요 — “baseline 가 심한 환자는 호전이 빠른가, 느린가?” 의 답.

§ 9.5.2 의 식 (9.24) Cholesky framework 그대로 적용: \(\Sigma_v = TT^\top\) 의 자유 모수 3 개 (\(T_{11}, T_{21}, T_{22}\)) 추정.

3.2 추정 결과 — 표 9.6

표 9.6 — Random Intercept and Trend 결과 (N = 437)

모수	Estimate	SE	\(Z\)	\(p <\)
Intercept	6.025	0.918	6.56	.001
Drug (0=placebo, 1=drug)	0.281	0.761	0.37	.71
Time (\(\sqrt{\text{week}}\))	-1.477	0.451	-3.27	.001
Drug × Time	-1.587	0.479	-3.31	.001
Cholesky 원소
Intercept (\(T_{11}\))	2.726	0.597
Intercept-time covariance (\(T_{21}\))	-0.829	0.352
Time (\(T_{22}\))	1.561	0.248

\(-2 \log L = 1227.43\)

직관 — 세 모형의 효과 추정 비교

모수	Fixed (Table 9.3)	Random Int (Table 9.4)	Random Int+Trend (Table 9.6)
Intercept	3.703	5.387	6.025
Drug	-0.405	-0.025	+0.281
Time	-1.113	-1.500	-1.477
Drug × Time	-0.418 (p=.11)	-1.015 (p<.001)	-1.587 (p<.001)

연속적 변화:

Intercept 가 점점 커짐: 잠재 분산 증가의 결과 (subject-specific scale 이 더 커짐).
Drug 주효과가 비유의 유지: baseline 차이 없음 — randomization 잘 됨.
Time 효과는 비교적 안정 (-1.11 → -1.50 → -1.48): placebo 의 시간 추세는 어떤 모형으로도 비슷.
Drug × Time 효과 점점 강해짐: -0.418 → -1.015 → -1.587. 가장 중요한 효과의 추정이 모형 정교화에 따라 더 명확해짐.

SE 의 변화:

Random intercept (SE 0.274) → Random intercept + trend (SE 0.479): SE 가 커짐. 모형이 더 유연해 추정 불확실성 증가.
그러나 효과 크기 증가가 SE 증가보다 빨라 z-statistic 은 -3.70 → -3.31 로 비슷하게 유지.
p-value 모두 < .001 로 강한 유의.

\(-2 \log L\): 1362.06 → 1249.74 → 1227.43. 모형 정교화에 따라 deviance 가 단조 감소.

3.3 LR 검정 — Random Trend 의 추가 가치

우도비 검정

\(H_0\): \(\sigma_{\upsilon_1}^2 = 0\) 그리고 \(\sigma_{\upsilon_0 \upsilon_1} = 0\) (기울기 분산 + 절편-기울기 공분산 모두 0).

\[ \chi^2_2 = -2 \log L_{\text{int only}} - (-2 \log L_{\text{int+trend}}) = 1249.74 - 1227.43 = 22.31 \]

자유도 2 (\(\sigma_{\upsilon_1}^2\) 와 \(\sigma_{\upsilon_0 \upsilon_1}\)), p < .001.

직관 — 호전 속도의 환자 차이도 무시할 수 없음

\(\chi^2_2 = 22.31\) 는 자유도 2 의 .001 quantile (약 13.82) 보다 훨씬 큼 → 강한 유의.

해석:

환자별 호전 속도가 통계적으로 유의하게 다름.
Random intercept 만으로는 이 차이를 흡수 못함.
→ Random intercept + trend 가 더 적절한 모형.

Boundary 문제: \(\sigma_{\upsilon_1}^2 = 0\) 도 boundary. 정확한 영분포는 \(0.5 \chi^2_1 + 0.5 \chi^2_2\) (혼합 분포). 그러나 검정 통계량이 너무 커 보정 후에도 강한 유의.

모형 선택의 결론: Random intercept + trend 모형이 NIMH 데이터에 가장 적합.

3.4 Cholesky 복원 — 식 (9.68)

식 (9.68) — \(\Sigma_v\) 의 복원

표 9.6 의 Cholesky 원소로부터:

\[ \widehat\Sigma_v = \widehat{T} \widehat{T}^\top = \begin{bmatrix} 2.73 & 0 \\ -0.83 & 1.56 \end{bmatrix} \begin{bmatrix} 2.73 & -0.83 \\ 0 & 1.56 \end{bmatrix} = \begin{bmatrix} 7.43 & -2.27 \\ -2.27 & 3.12 \end{bmatrix} \tag{9.68} \]

분산:

\(\widehat\sigma_{\upsilon_0}^2 = 7.43\) → \(\widehat\sigma_{\upsilon_0} = 2.73\) (절편 SD).
\(\widehat\sigma_{\upsilon_1}^2 = 3.12\) → \(\widehat\sigma_{\upsilon_1} = 1.77\) (기울기 SD).

공분산: \(\widehat\sigma_{\upsilon_0 \upsilon_1} = -2.27\).

상관:

\[ \widehat{r}_{\upsilon_0 \upsilon_1} = \frac{-2.27}{2.73 \cdot 1.77} = -0.47 \]

직관 — 절편-기울기 음의 상관 -0.47 의 임상적 의미

\(\widehat{r} = -0.47\) 은 임상적으로 매우 흥미로운 패턴:

절편이 큰 환자 (\(\upsilon_{0i} > 0\), baseline 가 심함) → 기울기가 음 (\(\upsilon_{1i} < 0\), 빠르게 호전).
절편이 작은 환자 (\(\upsilon_{0i} < 0\), baseline 가 가벼움) → 기울기가 양 (\(\upsilon_{1i} > 0\), 천천히 호전).

해석:

천장 효과 (ceiling effect): baseline 가 매우 심한 환자는 더 호전될 여지가 많다 (logit 척도에서). 가벼운 환자는 이미 거의 well 이라 추가 호전 여지가 적다.
약효의 표적: 약물이 가장 심한 환자에게 가장 큰 효과 — 임상시험의 자연스러운 결과.
회귀 향한 평균 (regression to the mean) 의 표현: 극단값에서 시작한 환자가 시간이 지나면서 평균으로 돌아가는 통계적 현상.

그룹 평균 + 분산 의 해석:

Placebo intercept 평균: \(\beta_0 = 6.03\), drug intercept 평균: \(\beta_0 + \beta_2 = 6.03 + 0.28 = 6.31\) (둘이 비슷, RCT 정합).
Placebo slope 평균: \(\beta_1 = -1.48\), drug slope 평균: \(\beta_1 + \beta_3 = -1.48 - 1.59 = -3.07\) (drug 가 약 2 배 빠른 호전).
Intercept SD = 2.73 (그룹 평균보다 큼! → 같은 그룹 내 환자 차이가 그룹 평균 차이를 압도).
Slope SD = 1.77 (이것도 큼 → 호전 속도 차이가 매우 다양).

→ “그룹 평균” 만 보고하는 것이 환자 변동을 얼마나 가리는지 보여줌. Mixed-effects 모형이 이 정보를 명시적으로 제공.

3.5 Empirical Bayes 절편-기울기 산점도 — Figure 9.14

환자별 \((\widehat\upsilon_{0i}, \widehat\upsilon_{1i})\) 의 분포

Figure 9.14 의 산점도 (placebo 와 drug 환자 별도 마커) 가 보여주는 패턴:

음의 상관 시각적 확인: 좌상→우하 방향의 점 분포 (Cholesky 복원 결과 -0.47 과 일치).
그룹 분리: drug 환자가 평균적으로 더 음의 기울기 영역에 분포 (빠른 호전).
개별 변동: 그룹 안에서도 광범위한 변동 — placebo 중에 매우 음의 기울기 (호전 빠름) 환자, drug 중에 양의 기울기 (악화) 환자 모두 존재.

음의 공분산의 의미 (저자 본문 인용):

“subjects with more negative intercepts (more likely to be well at baseline) have more positive slopes (less likely to improve over time).”

— 즉, baseline 에서 이미 가벼운 환자 (절편이 작음) 는 추가 호전 여지가 적어 기울기가 더 양수. 천장/바닥 효과의 자연스러운 결과.

이 산점도가 임상적으로 가치 있는 이유: 개별 환자의 trajectory 를 예측 가능. 새 환자의 baseline 만 보면 호전 속도를 어느 정도 예측 가능 (음의 상관 활용).

3.6 환자별 Trend Lines — Figure 9.15-9.16

식 — 개인별 trajectory

각 환자의 fitted probability:

\[ \widehat{P}(Y_{ij} = 1) = \frac{1}{1 + \exp[-(6.03 + 0.28 D_i - 1.48 \sqrt{W_j} - 1.59 D_i \sqrt{W_j} + \widehat\upsilon_{0i} + \widehat\upsilon_{1i} \sqrt{W_j})]} \]

각 환자의 EB 추정 \(\widehat\upsilon_{0i}, \widehat\upsilon_{1i}\) 대입.

직관 — “그룹 평균” 너머의 광범위한 trajectory

Figure 9.15 (placebo, N=108) 와 Figure 9.16 (drug, N=329) 의 trend line 묶음:

Placebo 그룹의 패턴 (Figure 9.15):

대부분 환자가 천천히 호전 (양의 기울기 또는 약한 음의 기울기).
일부 환자는 매우 빠르게 호전 (가파른 음의 기울기) — 약물 없이도 회복.
일부 환자는 끝까지 sick — 자연 회복 안 됨.

Drug 그룹의 패턴 (Figure 9.16):

대부분 환자가 가파르게 호전.
그러나 일부 환자는 끝까지 sick — 약물에 반응 안 함.

임상적 함의:

Personalized Medicine의 정량적 토대: 그룹 평균만 보면 “약물이 효과 있다” 만 알 수 있지만, individual trajectory 를 보면 어떤 환자가 약효 받고 누가 안 받는지 패턴이 드러남.
Treatment Failure 의 식별: Drug 그룹 중 끝까지 sick 인 환자 = 다른 약물 필요. 일찍 식별할 수 있으면 전환 결정 빠름.
Spontaneous Recovery 의 인지: Placebo 중 빠르게 호전한 환자 = 자연 회복 가능군. 이들에게 약물 처방의 비용/이익 재고.

이 분석이 RCT 결과의 임상 적용에서 핵심. “평균은 거짓말” — 환자 분포가 진짜 정보.

3.7 Marginal 적합 — Figure 9.17

Random intercept 와 random trend 의 marginal 차이

Figure 9.17 (random intercept + trend 모형의 marginal) 와 Figure 9.13 (random intercept 만의 marginal) 비교:

두 그림 모두 marginal 곡선이 관측 비율과 잘 일치.
그러나 random trend 모형은 환자 trajectory 의 광범위한 변동을 반영.

Marginalization 절차는 동일 (4 단계) 이지만, \(\widehat\Sigma_v\) 가 random intercept + trend 에서 다음으로 변경:

\[ \widehat{V}(y_i) = Z_i \widehat{T} \widehat{T}^\top Z_i^\top + \sigma^2 I_i \]

여기서 \(Z_i\) 가 시점 수 × 2 행렬 (\(Z_i = (1, \sqrt{W_j})\) 의 stacked form), \(\widehat{T} \widehat{T}^\top\) 가 식 (9.68) 의 \(\widehat\Sigma_v\).

저자의 결론 인용:

“the model fits the marginal proportions well. This isn’t surprising given that the previous random intercept model also fit these marginal proportions well.”

→ Marginal 적합 만으로는 두 모형을 구분 못함. 두 모형의 차이는 individual trajectory 와 SE 추정에서만 드러남. 모형 선택은 LR 검정과 임상적 해석으로 결정해야 한다.

4 세 모형의 종합 비교

비교 표

항목	Fixed (§9.7.1)	Random Int (§9.7.2)	Random Int+Trend (§9.7.3)
가정	시점 독립	환자 baseline 차이	+ 호전 속도 차이
자유 모수	4	5	7
Drug × Time	-0.418 (p=.11)	-1.015 (p<.001)	-1.587 (p<.001)
Drug × Time SE	0.262	0.274	0.479
Drug × Time Z	-1.60	-3.70	-3.31
환자 SD	—	2.116	2.726 (절편), 1.77 (기울기)
ICC	—	0.58	— (다차원)
\(-2 \log L\)	1362.06	1249.74	1227.43
LR vs 직전 모형	—	\(\chi^2_1 = 112.3\), p<.001	\(\chi^2_2 = 22.31\), p<.001
임상적 결론	약효 비유의	약효 유의 (subject-specific)	약효 유의 + 환자 trajectory 광범위

모형 선택의 교훈

이 NIMH 사례는 종단 GLMM 의 본질적 가치 를 시연:

모형이 결론을 바꿈: Fixed-effects 만 봤으면 “약효 입증 못함” 결론. Mixed-effects 로 정확한 SE 추정 후 강한 약효 발견.
환자 이질성이 RCT 의 핵심 정보: ICC = 0.58 같은 큰 값은 그룹 평균만 보면 놓침.
모형 정교화의 점진적 가치: Random intercept → Random intercept + trend 로 갈수록 효과 추정이 정확해짐. LR 검정으로 정량 평가.
Marginal vs Subject-specific 의 양립: 두 정보 모두 가치 있음. Group 비교는 marginal, 개인 의사결정은 subject-specific.
그래프의 위력: Empirical Bayes trajectory 가 전체 환자의 다양성을 시각화 — “평균” 너머의 정보 전달.

실무 권고: 종단 이항 데이터를 분석할 때 항상 세 모형 모두 적합:

Baseline (fixed-effects): 단순 비교용.
Random intercept: 환자 이질성 검정.
Random intercept + trend: 가장 일반적. LR 검정으로 random trend 의 가치 평가.

이 세 모형이 모두 같은 결론을 주면 결과가 강건. 다른 결론을 주면 어느 모형이 데이터에 맞는지 LR 검정 + 임상적 해석으로 결정.

5 응용 분야

분야	활용	NIMH 와의 평행
정신과 RCT	항우울제·항정신병 약물의 시간별 효과	본 시연
항암 임상	tumor response 추적 + 환자 trajectory	Empirical Bayes 산점도로 responder 식별
심부전 관리	NYHA class 호전 여부 + 환자별 진행	절편-기울기 음의 상관
재활 의학	운동 능력 회복 + 환자 변동	Random trend 가 재활 속도 차이 모형화
디지털 치료제	앱 사용 후 PHQ-9 호전 + 사용자 trajectory	개인 맞춤 알림 시기 결정
교육 RCT	학습 중재 효과 + 학생별 학습 곡선	학생별 절편 + 기울기 차이

6 코드 예시

6.1 Step 1: Random Intercept 모형 (R lme4 등가)

import numpy as np
import pandas as pd
import statsmodels.api as sm


# 앞 09-6 의 시뮬레이션 함수 재사용 + random intercept 추가
def simulate_nimh_with_random(seed: int = 2026) -> pd.DataFrame:
    """NIMH 데이터 + 환자 random intercept 추가 시뮬레이션."""
    rng = np.random.default_rng(seed)
    attrition = {
        "placebo": {0: 107, 1: 105, 2: 5, 3: 87, 4: 2, 5: 2, 6: 70},
        "drug":    {0: 327, 1: 321, 2: 9, 3: 287, 4: 9, 5: 7, 6: 265}
    }

    rows = []
    subject_id = 0
    for group, weeks in attrition.items():
        drug = 1 if group == "drug" else 0
        N = weeks[0]
        for i in range(N):
            # 환자 random intercept (table 9.4: sigma_v = 2.116)
            upsilon = rng.normal(0, 2.116)
            for week, count in weeks.items():
                if rng.random() < (count / N):
                    sqrt_w = np.sqrt(week)
                    # Table 9.4 추정값 적용
                    eta = 5.387 - 0.025 * drug - 1.500 * sqrt_w \
                          - 1.015 * drug * sqrt_w + upsilon
                    p = 1.0 / (1.0 + np.exp(-eta))
                    y = rng.binomial(1, p)
                    rows.append({
                        "subject": subject_id, "drug": drug,
                        "week": week, "sqrt_week": sqrt_w, "y": y
                    })
            subject_id += 1
    return pd.DataFrame(rows)


df = simulate_nimh_with_random()

# Random intercept 적합 — statsmodels variational Bayes
model = sm.BinomialBayesMixedGLM.from_formula(
    "y ~ drug + sqrt_week + drug:sqrt_week",
    random={"a": "0 + C(subject)"}, data=df)
result = model.fit_vb()

print("표 9.4 (Random Intercept) 와 비교:")
print(f"  Fixed effects: {result.fe_mean.round(3).tolist()}")
print(f"  sigma_v:       {np.sqrt(result.vc_mean[0]):.3f} (table: 2.116)")

R 등가 코드

library(lme4)
fit_int <- glmer(y ~ drug + sqrt_week + drug:sqrt_week + (1 | subject),
                 data = df, family = binomial, nAGQ = 20)
summary(fit_int)

# ICC 계산
sigma2_v <- as.numeric(VarCorr(fit_int)$subject)
icc <- sigma2_v / (sigma2_v + pi^2 / 3)
cat("ICC =", round(icc, 3), "\n")  # 표 9.4 의 0.58

R lme4::glmer(nAGQ = 20) 가 adaptive Gauss-Hermite quadrature 로 정확. statsmodels VB 는 빠르나 정확도 약간 낮음. 학습은 R 권고.

6.2 Step 2: Random Intercept + Trend 모형

# Note: statsmodels 는 random intercept + slope 의 동시 적합이 제한적.
# 본격 분석은 R lme4 또는 statsmodels.MixedLM (정규만) 또는 PyMC

# R 등가 코드 (권장):
r_code = """
library(lme4)
fit_full <- glmer(y ~ drug + sqrt_week + drug:sqrt_week
                       + (1 + sqrt_week | subject),
                  data = df, family = binomial,
                  control = glmerControl(optimizer = "bobyqa"))
summary(fit_full)

# Cholesky 추정 + Sigma_v 복원
VarCorr(fit_full)
# Sigma_v 의 직접 추출
Sigma_v <- VarCorr(fit_full)$subject
print(Sigma_v)

# 상관계수
attr(Sigma_v, "correlation")  # -0.47 예상

# LR test vs random intercept only
fit_int_only <- glmer(y ~ drug + sqrt_week + drug:sqrt_week + (1 | subject),
                      data = df, family = binomial)
anova(fit_int_only, fit_full)  # Chi^2 = 22.31, df = 2 예상
"""
print(r_code)

추정 안정성

Random intercept + trend 는 random intercept 만 보다 수렴이 어려울 수 있음:

모수 수 증가로 우도 표면이 평평한 영역 발생.
Optimizer 선택 (bobyqa, Nelder_Mead, nlminbwrap) 이 결과에 영향.
nAGQ = 1 (Laplace) 가 기본, nAGQ > 1 은 random intercept + slope 에서 지원 안 됨 (R lme4 제약).

수렴 안 되면:

Optimizer 변경.
glmmTMB 패키지 시도 (Laplace + AD).
Bayesian (brms + Stan) 으로 우회.

NIMH 의 경우 수렴이 잘 되어 표 9.6 결과 산출 가능.

6.3 Step 3: Cholesky 분해와 Σ_v 복원 (수치 검증)

import numpy as np


# 표 9.6 의 Cholesky 추정값
T_11 = 2.726   # intercept SD (Cholesky)
T_21 = -0.829  # intercept-time covariance (Cholesky)
T_22 = 1.561   # time SD (Cholesky)

T = np.array([[T_11, 0.0],
              [T_21, T_22]])

# 식 (9.68) — Sigma_v 복원
Sigma_v = T @ T.T

print("Cholesky factor T:")
print(T)
print("\nReconstructed Sigma_v (식 9.68):")
print(Sigma_v)
print(f"\n비교 — 표 9.6 의 보고값:")
print(f"  Var(intercept) = {Sigma_v[0,0]:.2f} (표: 7.43)")
print(f"  Cov(int, slope) = {Sigma_v[0,1]:+.2f} (표: -2.27)")
print(f"  Var(slope) = {Sigma_v[1,1]:.2f} (표: 3.12)")

# 표준편차와 상관계수
sd_int = np.sqrt(Sigma_v[0, 0])
sd_slope = np.sqrt(Sigma_v[1, 1])
corr = Sigma_v[0, 1] / (sd_int * sd_slope)

print(f"\nSD(intercept) = {sd_int:.2f} (표: 2.73)")
print(f"SD(slope)     = {sd_slope:.2f} (표: 1.77)")
print(f"Cor(int, slope) = {corr:+.2f} (표: -0.47)")

# 양정치 검증
eigvals = np.linalg.eigvalsh(Sigma_v)
print(f"\nEigenvalues: {eigvals}  → 양정치? {np.all(eigvals > 0)}")

결과 해석

식 (9.68) 의 모든 숫자가 정확히 재현됨. Cholesky factor 로부터 분산-공분산 행렬을 복원하는 절차가 § 9.5.2 의 일반 framework 의 직접 응용.

실무에서는 R lme4 등이 자동으로 VarCorr() 로 \(\Sigma_v\) 와 상관계수를 출력. 수동 계산은 모형 이해 + 검증 목적.

음수 상관계수 -0.47 가 양정치 보존 과 양립하는지 확인: \(\det(\Sigma_v) = 7.43 \cdot 3.12 - (-2.27)^2 = 23.18 - 5.15 = 18.03 > 0\) → 양정치. Cholesky framework 가 자동으로 보장.

6.4 Step 4: 세 모형의 -2 log L 비교 (LR test 직접 계산)

import numpy as np
from scipy import stats


# 표 9.3, 9.4, 9.6 의 -2 log L
neg2logL_fixed = 1362.06
neg2logL_int = 1249.74
neg2logL_int_trend = 1227.43

# LR test 1: Random Int vs Fixed (df = 1, sigma_v^2)
chi2_1 = neg2logL_fixed - neg2logL_int
p_1 = 1 - stats.chi2.cdf(chi2_1, df=1)
print(f"LR test (Random Int vs Fixed):")
print(f"  Chi^2 = {chi2_1:.2f}, df = 1, p = {p_1:.2e}")
print(f"  (표 9.4 의 X^2 = 112.3 와 일치)")

# LR test 2: Random Int+Trend vs Random Int (df = 2, sigma_v1^2 + cov)
chi2_2 = neg2logL_int - neg2logL_int_trend
p_2 = 1 - stats.chi2.cdf(chi2_2, df=2)
print(f"\nLR test (Random Int+Trend vs Random Int):")
print(f"  Chi^2 = {chi2_2:.2f}, df = 2, p = {p_2:.2e}")
print(f"  (표 9.6 의 X^2 = 22.31 와 일치)")

# Boundary 보정 (Self & Liang 1987)
# Random Int vs Fixed: 0.5 * chi2_0 + 0.5 * chi2_1
p_1_bdy = 0.5 * (1 - stats.chi2.cdf(chi2_1, df=1))
print(f"\nBoundary 보정 후 p (LR 1): {p_1_bdy:.2e}")
print(f"   → 너무 작아 결론 변화 없음")

# Random Int+Trend vs Random Int: 0.5 * chi2_1 + 0.5 * chi2_2
p_2_bdy = 0.5 * (1 - stats.chi2.cdf(chi2_2, df=1)) \
          + 0.5 * (1 - stats.chi2.cdf(chi2_2, df=2))
print(f"Boundary 보정 후 p (LR 2): {p_2_bdy:.2e}")
print(f"   → 결론 변화 없음 (강한 유의 유지)")

Boundary 보정의 실무적 의미

표준 LR 검정 (\(\chi^2_k\)) 은 모수 공간 안의 점에 대한 가설용. 분산 = 0 가설은 모수 공간 경계에 위치.

Self & Liang (1987) 결과: \(H_0: \sigma^2 = 0\) 의 정확한 영분포는 \(\frac{1}{2}\chi^2_{k-1} + \frac{1}{2}\chi^2_k\) 의 혼합. p-value 가 표준 \(\chi^2\) 보다 작아짐 (즉 더 쉽게 거부).

NIMH 데이터의 검정 통계량 (112.3, 22.31) 이 너무 커서 보정이 결론에 영향 없음. 그러나 검정 통계량이 임계값 근처일 때는 보정 필수.

R lme4 의 anova() 는 표준 \(\chi^2\) 를 보고 — 실제로는 약간 보수적. RLRsim 패키지가 정확한 분포 제공.

7 관련 주제

선행 지식

§ 9.7~9.7.1 — NIMH 데이터 소개 + Fixed-effects baseline (직전 sub-post)
§ 9.5.1-9.5.3 — Mixed-effects 모형 정의, ICC, Cholesky framework
§ 9.5.4-9.5.5 — Multilevel form (식 9.29-9.31)
§ 9.6 추정 — Marginal MLE, Empirical Bayes (식 9.48), Gauss-Hermite quadrature

후속 주제

§ 9.8 — Summary
Ch.10 GLMM 순서형 — IMPS 의 원래 7 점 척도 (이항화 없이)
Ch.14 결측 데이터 — NIMH attrition 의 informative dropout 검증
§ 4 정규 종단 MRM — 같은 framework 의 정규 버전

관련 개념

Gibbons & Bock (1987) — Random intercept + trend probit 모형 원전
Gibbons & Hedeker (1994) — NIMH 데이터의 probit 분석 (logit 의 대응)
Self & Liang (1987) — LR 검정의 boundary 분포
Verbeke & Molenberghs (2000) — 종단 데이터의 random effects 모형 종합
GEE (Ch.8) — Marginal vs subject-specific 의 다른 한쪽
mm-06 GLMM 이진 — AI Agent 비즈니스 예시 직관