1 들어가며 — 13-0 의 깊이

Ch.13 Overview (13-0) 에서 3-level data 의 큰 그림과 식 (13.1) 의 정의를 다뤘다. 본 sub-post 는 § 13.1 의 수학적 깊이 + § 13.1.1 의 NIMH illustration center clustering 분석.

한 줄 요약

“§ 13.1 = 식 (13.1) 의 행렬 표현 해부 — \(Z_i\) 가 cluster intercept (모든 row 1) + subject blocks 의 hierarchical 구조, \(\Sigma_i\) 가 block-diagonal (cluster + subject blocks). EAP estimator + score + EM + Fisher scoring (§ 4.5 의 직접 확장). § 13.1.1 NIMH illustration = 9 treatment centers + placebo/chlorpromazine + IMPS Item 79. 2-level 모형 (Table 13.2): Drug × Time -0.564 (p<.001), chlorpromazine 6 주 호전이 placebo 의 약 2.6 배. 3-level 확장 (Table 13.4): center random effect 추가, \(\widehat\sigma_\gamma^2 = 0.039\) (비유의), fixed effects 거의 같음 — 본 사례에서 cluster 효과 작음. Table 13.5 variance decomposition: baseline subject 32% → week 6 74%, center 1.66-4.39% (Donner 1982 일치). 그러나 일반 권고: 3-level 적합 후 cluster variance 검정 — 작으면 2-level 로 reduce.”

2 § 13.1 — 식 (13.1) 의 깊이

2.1 행렬 표현 해부

식 (13.1) — 3-level linear

\[ y_i = X_i \beta + Z_i \upsilon_i^* + \varepsilon_i \tag{13.1} \]

표기 (13-0 복습):

\(i = 1, \ldots, N\): cluster.
\(j = 1, \ldots, n_i\): subject in cluster \(i\).
\(k = 1, \ldots, n_{ij}\): observation for subject \(j\).
\(N_i = \sum_j n_{ij}\): cluster \(i\) 총 관측.
\(R_i = n_i \times r + 1\): cluster \(i\) 의 random effects 차원.

\(Z_i\) Matrix 의 Block 구조

\(Z_i\) (\(N_i \times R_i\)) 의 정확한 형태:

\[ Z_i = \begin{bmatrix} 1 & Z_{i1} & 0 & 0 & \cdots & 0 \\ 1 & 0 & Z_{i2} & 0 & \cdots & 0 \\ 1 & 0 & 0 & Z_{i3} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 0 & 0 & 0 & \cdots & Z_{i n_i} \end{bmatrix} \]

각 row block 이 한 subject 의 관측들. Column 구조:

첫 column: cluster random intercept (\(\gamma_i\)) 의 design — 모든 row 1.
나머지 columns: subject 별 design (\(Z_{ij}\), \(r\) columns), block-diagonal.

→ 같은 subject 의 관측이 같은 \(\gamma_i\) + 같은 \(\upsilon_{ij}\) 영향받음. → 다른 subject (같은 cluster) 는 같은 \(\gamma_i\) 만 공유, \(\upsilon_{ij}\) 는 별도.

예시 — 한 cluster 에 3 subjects, subject 마다 \(r = 2\) random effects (intercept + slope), 시점 마다 4 관측:

\(N_i = 3 \times 4 = 12\) rows.
\(R_i = 3 \times 2 + 1 = 7\) columns.
\(Z_i\) 가 \(12 \times 7\) matrix:
- Subject 1 의 4 row 가 columns (1, 2, 3) 만 non-zero.
- Subject 2: columns (1, 4, 5).
- Subject 3: columns (1, 6, 7).
- Column 1 은 모두 1 (cluster intercept).

이 sparse block-like 구조가 계산 효율성 의 핵심 — 같은 subject 안에서만 random effects 공유.

2.2 \(\Sigma_i\) Block-Diagonal 구조

Random effects 분산-공분산

\(\upsilon_i^*\) 의 분산-공분산 \(\Sigma_i\) (\(R_i \times R_i\)):

\[ \Sigma_i = \begin{bmatrix} \sigma_\gamma^2 & 0 & 0 & \cdots & 0 \\ 0 & \Sigma_\upsilon & 0 & \cdots & 0 \\ 0 & 0 & \Sigma_\upsilon & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \Sigma_\upsilon \end{bmatrix} \]

표기:

첫 block (1×1): \(\sigma_\gamma^2\) — cluster random intercept 분산.
나머지 blocks (\(r \times r\)): \(\Sigma_\upsilon\) — 모든 subject 가 같은 분산 (i.i.d. 가정).
비대각 blocks: 0 (cluster vs subject, subject vs subject 독립).

직관 — Block-Diagonal 의 가정

가정 1 — Cluster effect 와 subject effects 독립:

\(\gamma_i\) 와 \(\upsilon_{ij}\) 의 covariance = 0.
“Cluster 의 평균 효과” 와 “환자별 변동” 이 별도 메커니즘.

가정 2 — Subject effects 독립 (within cluster):

\(\upsilon_{ij}\) 와 \(\upsilon_{ij'}\) (\(j \neq j'\)) 의 covariance = 0.
같은 cluster 의 환자들이 conditional independent (cluster effect 조건부).

가정 3 — Subject effects 동질:

모든 \(\upsilon_{ij} \sim \mathcal{N}_r(0, \Sigma_\upsilon)\) (같은 분산).
환자별 분산 차이 없음.

가정 위반 시 처리:

Cluster-subject 상관: 더 일반 모형 (cross-level random effects).
Subject 간 추가 상관: 가족 효과 등 추가 random effect.
환자별 다른 분산: scaling terms (§ 10.2.2 framework).

block-diagonal 의 계산 가치:

행렬 inverse 가 block 별 inverse 의 합.
Determinant 가 block determinant 의 곱.
→ 매우 큰 cluster (\(n_i = 100\)) 에서도 효율적 계산.

2.3 Joint Distribution

\(y_i\) 와 \(\upsilon^*\) 의 joint multivariate normal

\[ \begin{bmatrix} y_i \\ \upsilon_i^* \end{bmatrix} \sim \mathcal{N}\left(\begin{bmatrix} X_i \beta \\ 0 \end{bmatrix}, \begin{bmatrix} Z_i \Sigma_i Z_i^\top + \sigma_\varepsilon^2 I_i & Z_i \Sigma_i \\ \Sigma_i Z_i^\top & \Sigma_i \end{bmatrix}\right) \]

표기:

\(V(y_i) = Z_i \Sigma_i Z_i^\top + \sigma_\varepsilon^2 I_i\) — marginal 분산.
\(\text{Cov}(y_i, \upsilon_i^*) = Z_i \Sigma_i\).
\(V(\upsilon_i^*) = \Sigma_i\).

직관 — Marginal 분산의 의미

\(V(y_i) = Z_i \Sigma_i Z_i^\top + \sigma_\varepsilon^2 I_i\) 의 구조:

\(Z_i \Sigma_i Z_i^\top\): random effects 의 marginal 기여.
- 같은 subject 의 관측이 \(\sigma_\gamma^2 + \text{subj effects}\) 로 covariance.
- 다른 subject (같은 cluster) 의 관측이 \(\sigma_\gamma^2\) 만 covariance.
- 다른 cluster 의 관측이 covariance 0.
\(\sigma_\varepsilon^2 I_i\): residual 의 대각.

→ 자연스러운 hierarchical covariance structure.

ICC 의 정확한 도출:

같은 subject, 다른 시점 (subject intercept만 random):

\[ \text{Cov}(y_{ijk}, y_{ijk'}) = \sigma_\gamma^2 + \sigma_{\upsilon_0}^2 \]

같은 cluster, 다른 subject:

\[ \text{Cov}(y_{ijk}, y_{ij'k}) = \sigma_\gamma^2 \]

다른 cluster:

\[ \text{Cov}(y_{ijk}, y_{i'j'k'}) = 0 \]

분산:

\[ V(y_{ijk}) = \sigma_\gamma^2 + \sigma_{\upsilon_0}^2 + \sigma_\varepsilon^2 \]

Cluster ICC:

\[ \text{ICC}_{\text{cluster}} = \frac{\sigma_\gamma^2}{\sigma_\gamma^2 + \sigma_{\upsilon_0}^2 + \sigma_\varepsilon^2} \]

Subject ICC (조건부, cluster 같음):

\[ \text{ICC}_{\text{subject}} = \frac{\sigma_\gamma^2 + \sigma_{\upsilon_0}^2}{\sigma_\gamma^2 + \sigma_{\upsilon_0}^2 + \sigma_\varepsilon^2} \]

→ 13-0 의 ICC 분해의 정확한 도출.

2.4 EAP Estimator + Score

EAP estimator (식)

Empirical Bayes (사후 평균):

\[ \bar\upsilon^* = [Z_i^\top (\sigma_\varepsilon^2 I_i)^{-1} Z_i + \Sigma_i^{-1}]^{-1} Z_i^\top (\sigma_\varepsilon^2 I_i)^{-1} (y_i - X_i \beta) \]

사후 분산:

\[ \Sigma_{\upsilon^* \mid y_i} = [Z_i^\top (\sigma_\varepsilon^2 I_i)^{-1} Z_i + \Sigma_i^{-1}]^{-1} \]

Score 도출

Marginal log-likelihood \(\log L = \sum_i \log h(y_i)\).

\(\beta\) 에 대한 score:

\[ \frac{\partial \log L}{\partial \beta} = \sigma_\varepsilon^{-2} \sum_i X_i^\top \mu_i \]

\(\Sigma^*\) (분산 모수 vector) 에 대한 score:

\[ \frac{\partial \log L}{\partial \Sigma^*} = \frac{1}{2} \sum_i D_i^\top G_i^\top \text{vec}[\Sigma_i^{-1} (\Sigma_{\upsilon^* \mid y_i} + \bar\upsilon^* (\bar\upsilon^*)^\top - \Sigma_i) \Sigma_i^{-1}] \]

\(\sigma_\varepsilon^2\) 에 대한 score:

\[ \frac{\partial \log L}{\partial \sigma_\varepsilon^2} = \frac{1}{2 \sigma_\varepsilon^4} \sum_i \{-n_i \sigma^2 + \mu_i^\top \mu_i + \text{tr}[\Sigma_{\upsilon^* \mid y_i} Z_i^\top Z_i]\} \]

with \(\mu_i = y_i - Z_i \bar\upsilon^* - X_i \beta\).

직관 — § 4.5 의 직접 확장

저자 본문 인용:

“Note that these derivatives are functionally the same as those for the 2-level model. Thus, the same approach using the EM algorithm and Fisher scoring solution can be used.”

→ § 4.5 의 toolkit 그대로 + 행렬 차원 증가.

EM 알고리즘:

E-step: \(\bar\upsilon^*\) 와 \(\Sigma_{\upsilon^* \mid y_i}\) 계산 (위 식).
M-step: \(\beta, \Sigma^*, \sigma_\varepsilon^2\) score 식으로 갱신.
수렴까지 반복.

Fisher scoring:

정보 행렬 = \(E[-\partial^2 \log L / \partial \theta \partial \theta^\top]\).
Newton-Raphson update: \(\theta^{new} = \theta^{old} + I^{-1} \nabla \log L\).

계산 효율:

2-level: \(N\) subject 별 행렬 연산.
3-level: \(N\) cluster 별 행렬 연산 (cluster 안 subjects 결합).
→ cluster size 가 클 때 더 큰 행렬 inversion (block-diagonal 활용으로 효율).

Software:

R lme4::lmer(y ~ x + (1 | cluster) + (1 | cluster:subject)): REML/ML.
R nlme::lme(random = ~1 | cluster/subject): nested syntax.
SAS PROC MIXED with two RANDOM statements.
모두 같은 EM/Fisher scoring 알고리즘.

3 § 13.1.1 — NIMH Center Clustering Illustration

3.1 데이터 — IMPS Item 79

NIMH Schizophrenia Collaborative Study (Ch.9 와 같은 데이터의 continuous version)

Source: NIMH (Lorr & Klett 1966).

Outcome: IMPS Item 79 (“severity of illness”), 7 점 척도:

Normal, not at all ill
Borderline mentally ill
Mildly ill
Moderately ill
Markedly ill
Severely ill
Among the most extremely ill

Design:

9 treatment centers participated.
4 conditions: placebo + 3 drugs (chlorpromazine, fluphenazine, thioridazine).
본 분석: placebo + chlorpromazine 만 (3 약물 효과 비슷).
Time: 7 weeks (0, 1, 2, 3, 4, 5, 6).
Time transformation: square root of week (linearize).

Table 13.1 — 시점별 표본

Treatment	Week 0	Week 1	Week 2	Week 3	Week 4	Week 5	Week 6
Placebo	110	108	5	89	2	2	72
Chlorpromazine	110	108	3	96	4	5	87

관찰: 매우 unbalanced — week 2, 4, 5 에 표본 매우 적음.

직관 — Unbalanced Data 의 처리

NIMH 데이터의 attrition + irregular measurement:

Week 0, 1, 3, 6: 주요 measurement 시점 (~70-110 명).
Week 2, 4, 5: irregular (~2-9 명).

Unbalanced data 의 mixed-effects 처리:

같은 환자의 시점 수 다양 (1 ~ 7).
같은 시점의 환자 수 매우 다양.
→ Mixed-effects 가 자연 처리 (variable cluster size).
Standard ANOVA 또는 repeated measures 는 처리 어려움.

Sqrt(week) 변환의 동기:

IMPS rating 의 시간 추세가 비선형 (early week 빠른 호전 → late week 안정화).
\(\sqrt{\text{week}}\) 변환이 선형화.
→ 단일 slope parameter 로 시간 효과 모형화.

§ 9.7 (Ch.9) 의 NIMH 이항 분석에서 본 것과 같은 변환.

3.2 2-Level Model — Table 13.2

2-Level Mixed-Effects 결과 (Table 13.2)

모형:

Level-1 (within-subject): \(y_{ijk} = b_{0i} + b_{1i} \sqrt{\text{week}_k} + \varepsilon_{ijk}\)
Level-2 (between-subject): \(b_{0i} = \beta_0 + \beta_2 \text{Drug}_i + \upsilon_{0i}\), \(b_{1i} = \beta_1 + \beta_3 \text{Drug}_i + \upsilon_{1i}\)

Parameter	Estimate	SE	\(p <\)
Placebo intercept (\(\beta_0\))	5.345	0.085	.001
Placebo vs Chlorpromazine (\(\beta_2\))	0.029	0.120	.812
Placebo slope (\(\beta_1\))	-0.343	0.066	.001
Drug × Time (\(\beta_3\))	-0.564	0.091	.001
Subject intercept SD²	0.342	0.082
Subject slope SD²	0.224	0.044
Subject intercept-slope cov	0.007	0.047
Error SD²	0.571	0.042

직관 — 임상적 해석

Placebo intercept (\(\widehat\beta_0 = 5.345\)): baseline severity = “between severely (6) and markedly ill (5)”. 강한 유의성.

Drug 주효과 (\(\widehat\beta_2 = 0.029\), p=.812): chlorpromazine baseline 5.345 + 0.029 = 5.374 — placebo 와 거의 같음 (RCT randomization 정합).

Placebo slope (\(\widehat\beta_1 = -0.343\), p<.001): sqrt(week) 1 단위 당 0.343 unit 호전. 6 주 후 호전:

\[ \Delta_{\text{placebo}} = -0.343 \times \sqrt{6} = -0.84 \text{ units} \]

→ 베이스라인 5.345 → week 6 = 5.345 - 0.84 = 4.50 (moderately to markedly ill).

Drug × Time (\(\widehat\beta_3 = -0.564\), p<.001): chlorpromazine 의 추가 호전. 총 chlorpromazine slope:

\[ \widehat\beta_1 + \widehat\beta_3 = -0.343 - 0.564 = -0.907 \]

6 주 후 호전:

\[ \Delta_{\text{drug}} = -0.907 \times \sqrt{6} = -2.22 \text{ units} \]

→ Chlorpromazine: 5.374 - 2.22 = 3.15 (mildly ill).

약물의 임상적 가치:

Chlorpromazine 의 추가 호전: 2.22 - 0.84 = 1.38 units.
“Severely → moderately/mildly” 로의 결정적 변화.
강한 통계적 + 임상적 유의성.

환자 변동성:

Subject intercept SD² = 0.342, SD = 0.585.
±2 SD: 환자별 baseline 이 ±1.17 unit (일부 환자가 평균보다 매우 심각/덜 심각).
Subject slope SD² = 0.224, SD = 0.473.
±2 SD: 환자별 호전 속도가 ±0.95 unit/sqrt(week).

절편-기울기 상관: \(r = 0.007 / \sqrt{0.342 \times 0.224} = 0.024\) — 매우 작음.

→ “처음 심각한 환자가 더 빠르게 호전” 의 패턴 없음 (Ch.9 § 9.7.3 의 이항 분석 -0.47 와 다름 — continuous version 의 metric 차이).

3.3 3-Level Extension — Table 13.3-13.4

9 Centers, Variable Sample Sizes (Table 13.3)

Treatment	Center 1	2	3	4	5	6	7	8	9	Total
Placebo	13	20	13	15	13	7	10	10	6	107
Chlorpromazine	9	22	8	18	15	9	10	12	7	110

→ Center 별 6-22 환자 — unbalanced.

Table 13.4 — 3-Level Mixed-Effects 결과

3-Level 모형: 식 (13.1) 에 \(\gamma_i \sim \mathcal{N}(0, \sigma_\gamma^2)\) 추가.

Parameter	Estimate	SE	\(p <\)
Placebo intercept (\(\beta_0\))	5.339	0.106	.001
Drug (\(\beta_2\))	0.038	0.116	.741
Placebo slope (\(\beta_1\))	-0.341	0.066	.001
Drug × Time (\(\beta_3\))	-0.566	0.091	.001
Subject intercept SD²	0.285	0.078
Subject slope SD²	0.225	0.044
Subject intercept-slope cov	0.025	0.045
Center SD² (\(\sigma_\gamma^2\))	0.039	0.031	(NS)
Error SD²	0.570	0.040

직관 — 2-Level vs 3-Level 비교

Parameter	2-Level	3-Level	차이
\(\beta_0\)	5.345	5.339	거의 같음
\(\beta_2\)	0.029	0.038	거의 같음
\(\beta_1\)	-0.343	-0.341	거의 같음
\(\beta_3\)	-0.564	-0.566	거의 같음
Subject intercept SD²	0.342	0.285	감소
Subject slope SD²	0.224	0.225	같음
Center SD²	—	0.039	추가

핵심 관찰:

Fixed effects 거의 동일 — 본 사례에서 cluster 효과가 작아 추정 변화 미미.
Subject intercept SD² 감소 (0.342 → 0.285) — 일부 변동이 center variance 로 흡수.
Center SD² = 0.039 — 작음. SD = 0.198 (< 1/5 unit).
Center SD² 비유의 (SE = 0.031) — Wald test 에 의하면 NS.

임상적 해석:

Center 효과의 SD = 0.198 unit:
- 가장 strict center vs 가장 lenient center 의 평균 차이가 약 ±0.4 unit.
- 임상적으로 미미.
본 NIMH 데이터에서는 9 centers 가 비슷한 의료 표준 + protocol.
→ 2-level 분석으로 충분.

저자 본문 인용:

“Inspection of Table 13.4 reveals very similar results to those obtained earlier ignoring the center effect. The variance attributable to centers is not quite statistically significant.”

→ 본 사례에서 cluster effect 무시해도 결과 큰 차이 없음.

그러나 일반적 권고:

3-level 데이터면 일단 3-level 모형 적합.
Cluster variance 검정.
작으면 2-level 으로 reduce 정당.
크면 3-level 유지 + 정확한 SE.

3.4 Variance Components 시점별 — Table 13.5

시점별 분산 비율 (% of total)

Timepoint	Error	Subject	Center
Baseline	63.72	31.89	4.39
Week 1	48.73	47.91	3.36
Week 2	40.27	56.96	2.77
Week 3	34.41	63.22	2.37
Week 4	30.08	67.85	2.07
Week 5	26.74	71.42	1.84
Week 6	24.08	74.26	1.66

직관 — Subject Variance 의 시간 증가

시점에 따른 패턴:

Error (level-1, 잔차): 일정 (절대값) — 시간 무관 가정. 비율은 감소.
Subject (level-2): 시간 따라 증가 — slope variance 의 효과.
Center (level-3): 일정 (절대값) — 시간 무관 가정.

Subject 분산 증가의 메커니즘:

같은 환자의 시점 \(k\) 에서 잠재 variance:

\[ V_{\text{subject}}(t) = \sigma_{\upsilon_0}^2 + 2 t \sigma_{\upsilon_{01}} + t^2 \sigma_{\upsilon_1}^2 \]

(t = sqrt(week)).

\(t = 0\) (baseline): \(V = \sigma_{\upsilon_0}^2 = 0.285\).
\(t = \sqrt{6}\) (week 6): \(V = 0.285 + 2 \times \sqrt{6} \times 0.025 + 6 \times 0.225 = 0.285 + 0.122 + 1.35 = 1.76\).

→ Slope variance 의 효과로 subject variance 가 시간에 따라 ~6 배 증가.

시점별 ICC 의 시각화:

Baseline: 환자별 차이 32% — 환자가 처음에는 비슷.
Week 6: 환자별 차이 74% — 시간이 지남에 따라 환자별 호전 속도 차이 누적.

Center variance 의 일정함:

Center intercept 만 random (slope 없음).
→ Center 효과가 시간 무관.
대안: random center slope 추가 — 시간 따라 center 효과 커짐 가능.
본 분석에서는 단순화 (center variance 자체가 작음).

Donner (1982) 와의 일치:

Center ICC 1.66-4.39% — 다른 cluster studies 와 비슷한 범위.
의료 multi-center trials 의 일반적 cluster effect.
→ 모형 설계 단계에서 이 정도 cluster ICC 예상 가능.

4 Random vs Fixed Center 의 결정

두 framework 의 비교

저자 본문 인용:

“Treatment of the center effects as a random term in the model means that the specific centers used in the study are considered to be a representative sample from a larger population of potential centers. Conversely, if interest is only in making inferences about the specific centers of a dataset (e.g., is there a difference between center A and center B?), the center could then be regarded as a ‘fixed’ and not a ‘random’ effect in the model.”

Random center:

가정: 9 centers 가 더 큰 populationcenter 의 sample.
추정: center variance (\(\sigma_\gamma^2\)) — 모집단의 변동.
추론: 다른 center 에서도 유사 결과 기대.
사용 시점: external validity, 정책 (다른 hospital 도입 가능성).

Fixed center:

가정: 9 centers 자체에 대한 추론.
추정: 8 dummy variables (center 1-8 vs reference).
추론: 이 specific centers 의 차이만.
사용 시점: internal validity, 특정 hospital 비교.

직관 — 권고 절차

Random center 권고 시점:

External generalization 이 목적 — “다른 hospital 에서도 비슷한 효과?”
Center 가 sample 로 인식: 9 centers 가 가능한 모든 center 의 부분.
Cluster-level covariate 분석: center size, urban/rural 등 covariate 와 상호작용.
자료가 unbalanced: random effects 가 불균형 처리 자연.

Fixed center 권고 시점:

Specific centers 의 비교 — “Center A 가 B 보다 더 나은가?”
Center 가 population: 모든 가능한 center 가 데이터에 포함.
Center 수 매우 적음 (5 미만): random effects 추정 불안정.

NIMH 의 9 centers:

9 centers 가 미국 전체 정신과 hospital 의 sample.
Random effects 권고 (저자의 선택).
장기 정책 (다른 center 에서 약물 도입) 평가.

Fixed centers 의 한계:

9 centers × 추가 dummy 8 개 = 모수 8 개 추가.
Cluster-level covariate (center size) 와 confounded.
Marginal effect 만 추정 (center-specific 정책 분석 불가).

현대 표준 — Random center:

Multi-center clinical trials 의 표준.
ICH E9 (Statistical Principles for Clinical Trials) 권고.
Mixed-effects framework 의 기본 가정.

5 응용 분야

분야	Cluster (level 3)	Subject (level 2)	Observation (level 1)
Multi-center 임상시험	Treatment center	Patient	Visit
Hospital 비교	Hospital	Patient	Daily measure
학교 효과	School	Student	Annual test
가족 연구	Family	Member	Repeated measure
다국가 RCT	Country	Patient	Time
Spatial 연구	Region	County	Year
회사 효과	Company	Employee	Performance

6 코드 예시

6.1 Step 1: NIMH-like 3-Level 데이터 시뮬레이션

library(lme4)
library(dplyr)


# NIMH-like 3-level 데이터 시뮬레이션
set.seed(2026)
n_centers <- 9
center_sizes <- c(13, 20, 13, 15, 13, 7, 10, 10, 6)  # Table 13.3 placebo
                                                    # (chlorpromazine 도 비슷)
n_total_per_group <- sum(center_sizes)

# 모수 (Table 13.4 와 비슷)
beta_0 <- 5.339
beta_drug <- 0.038
beta_time <- -0.341
beta_drug_time <- -0.566

sigma_v0 <- sqrt(0.285)  # subject intercept SD
sigma_v1 <- sqrt(0.225)  # subject slope SD
sigma_gamma <- sqrt(0.039)  # center SD
sigma_e <- sqrt(0.570)  # residual SD

# Center random effects
gamma_i <- rnorm(n_centers, 0, sigma_gamma)

# 데이터 생성
df <- data.frame()
subject_id_counter <- 0
for (i in 1:n_centers) {
  for (group in c("placebo", "chlorpromazine")) {
    drug_indicator <- ifelse(group == "chlorpromazine", 1, 0)
    n_subj <- center_sizes[i]  # placebo, chlorpromazine 비슷
    for (j in 1:n_subj) {
      subject_id_counter <- subject_id_counter + 1
      # Subject random effects
      u <- mvtnorm::rmvnorm(1, mean = c(0, 0),
                              sigma = matrix(c(0.285, 0.025, 0.025, 0.225), 2))
      v0 <- u[1]
      v1 <- u[2]

      # Time points (NIMH 의 4 주요 시점만 단순화)
      for (week in c(0, 1, 3, 6)) {
        sqrt_week <- sqrt(week)
        eta <- beta_0 + beta_drug * drug_indicator +
               beta_time * sqrt_week +
               beta_drug_time * drug_indicator * sqrt_week +
               gamma_i[i] + v0 + v1 * sqrt_week
        y <- eta + rnorm(1, 0, sigma_e)
        df <- rbind(df, data.frame(
          center = i, subject = subject_id_counter,
          group = group, drug = drug_indicator,
          week = week, sqrt_week = sqrt_week, y = y
        ))
      }
    }
  }
}

cat("전체 관측:", nrow(df), "\n")
cat("Subjects:", length(unique(df$subject)), "\n")
cat("Centers:", length(unique(df$center)), "\n")

시뮬레이션 검증

논문의 추정값으로 데이터 생성 — 다시 적합하면 추정값 회복 가능.

검증 포인트:

환자 수 ~ 217 (Table 13.3 의 placebo 107 + chlorpromazine 110).
각 환자가 4 시점 (단순화).
9 centers, center 별 unbalanced.

6.2 Step 2: 2-Level vs 3-Level 적합 비교

# 2-Level mixed-effects (Table 13.2 재현)
fit_2level <- lmer(y ~ drug * sqrt_week + (sqrt_week | subject),
                   data = df, REML = TRUE)
summary(fit_2level)

# 3-Level mixed-effects (Table 13.4 재현)
fit_3level <- lmer(y ~ drug * sqrt_week + (sqrt_week | subject) +
                   (1 | center),
                   data = df, REML = TRUE)
summary(fit_3level)

# Fixed effects 비교
cat("\nFixed effects 비교:\n")
print(round(rbind(
  "2-level" = fixef(fit_2level),
  "3-level" = fixef(fit_3level)
), 3))

# Variance components 비교
cat("\nVariance components 비교:\n")
cat("2-level:\n")
print(VarCorr(fit_2level))
cat("3-level:\n")
print(VarCorr(fit_3level))

# LR test (cluster effect 의 유의성)
lr_stat <- 2 * (logLik(fit_3level) - logLik(fit_2level))
cat("\nLR test (3-level vs 2-level):\n")
cat("  chi^2 =", round(lr_stat, 2), "\n")
cat("  p (boundary 보정 0.5 * chi^2_1) =",
    format.pval(0.5 * (1 - pchisq(lr_stat, 1)), digits = 3), "\n")

2-Level vs 3-Level 의 검증

예상 결과 (Table 13.2 vs 13.4 재현):

Fixed effects 거의 동일 (드물게 0.01-0.05 차이).
Subject intercept SD² 약간 감소 (0.342 → 0.285).
Center SD² 추가 (~0.039).

LR test:

\(H_0: \sigma_\gamma^2 = 0\) (3-level 불필요).
Boundary 문제 — Self-Liang 1987 의 0.5 * \(\chi^2_1\) 분포.
p-value 가 작으면 (< .05) cluster effect 유의 → 3-level 권고.

시뮬레이션의 한계:

정확한 Table 13.2/13.4 재현 어려움 (representative sample 의 변동).
패턴 (방향성, 크기) 만 비슷하면 검증 성공.

6.3 Step 3: Variance Components 시점별 분해 (Table 13.5)

# Variance components 시점별 분해 (Table 13.5 재현)
weeks_to_analyze <- c(0, 1, 2, 3, 4, 5, 6)

# 추정 분산 추출
sigma_e2 <- sigma(fit_3level)^2
sigma_g2 <- as.numeric(VarCorr(fit_3level)$center)
varcorr_subj <- VarCorr(fit_3level)$subject
sigma_v0_2 <- varcorr_subj[1, 1]
sigma_v1_2 <- varcorr_subj[2, 2]
sigma_v01 <- varcorr_subj[1, 2]

# 시점별 subject variance: V(t) = sigma_v0^2 + 2*t*sigma_v01 + t^2*sigma_v1^2
# (t = sqrt(week))
result <- data.frame(timepoint = character(), error = numeric(),
                     subject = numeric(), center = numeric())

for (week in weeks_to_analyze) {
  t <- sqrt(week)
  v_subject <- sigma_v0_2 + 2 * t * sigma_v01 + t^2 * sigma_v1_2
  v_total <- sigma_e2 + v_subject + sigma_g2
  result <- rbind(result, data.frame(
    timepoint = ifelse(week == 0, "Baseline", paste("Week", week)),
    error = round(100 * sigma_e2 / v_total, 2),
    subject = round(100 * v_subject / v_total, 2),
    center = round(100 * sigma_g2 / v_total, 2)
  ))
}

print(result)

Table 13.5 의 직접 계산

위 코드가 § 13.1.1 의 Table 13.5 의 정확한 산출 방식.

핵심 관찰:

Baseline: error 비율 가장 큼 (subject slope 영향 없음).
Week 6: subject 비율 가장 큼 (slope variance 누적).
Center 비율: 일정한 절대값, 비율은 시간에 따라 감소 (total variance 증가 때문).

시각화 권고:

library(tidyr)
library(ggplot2)

result_long <- pivot_longer(result, cols = c(error, subject, center),
                             names_to = "component", values_to = "percent")

ggplot(result_long, aes(x = timepoint, y = percent, fill = component)) +
  geom_bar(stat = "identity", position = "stack") +
  scale_fill_manual(values = c("error" = "lightgray",
                                 "subject" = "steelblue",
                                 "center" = "tomato")) +
  labs(x = "Time", y = "% of Total Variance",
       title = "Variance Components Across Time")

→ 시간에 따른 variance decomposition 의 시각화 — 정책·임상 보고서에 활용.

6.4 Step 4: Empirical Bayes Center Effects

# Center 별 EB random effects
ranef_center <- ranef(fit_3level, condVar = TRUE)$center

# Caterpillar plot
library(ggplot2)
center_df <- data.frame(
  center = 1:nrow(ranef_center),
  estimate = ranef_center[, 1],
  se = sqrt(attr(ranef(fit_3level, condVar = TRUE)$center, "postVar")[1, 1, ])
)
center_df$rank <- order(center_df$estimate)
center_df <- center_df[order(center_df$estimate), ]
center_df$x <- 1:nrow(center_df)

ggplot(center_df, aes(x = x, y = estimate)) +
  geom_point(size = 3, color = "steelblue") +
  geom_errorbar(aes(ymin = estimate - 1.96 * se,
                    ymax = estimate + 1.96 * se),
                width = 0.3, color = "steelblue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "Center Rank", y = "EB Estimate (gamma_i)",
       title = "Empirical Bayes Center Effects (with 95% CI)")

# Center 별 fitted intercept (overall + center effect)
overall_intercept <- fixef(fit_3level)["(Intercept)"]
center_df$fitted_intercept <- overall_intercept + center_df$estimate
print(center_df[, c("center", "estimate", "se", "fitted_intercept")])

Caterpillar Plot 의 정책 활용

Center ranking:

가장 높은 estimate: 평균보다 baseline 자살률 ↑ — 의료 quality 점검.
가장 낮은 estimate: baseline 자살률 ↓ — best practice 학습.
95% CI 가 0 을 포함하지 않으면 유의하게 다름.

Shrinkage 효과:

작은 center (n=6, 7): CI 가 넓음 — 추정 불안정.
큰 center (n=20, 22): CI 가 좁음 — 추정 정확.
→ EB 가 작은 center 를 모집단 평균으로 수축.

임상 의사결정:

Outlier center 식별 → 추가 조사.
Best practice center → 다른 center 에 transfer.
Center 별 자원 배분 결정.

§ 12.4.2 (Hospital mortality, Thomas 1992) 와 동일 framework — Hospital ranking 의 mixed-effects 표준.

7 관련 주제

선행 지식

Ch.13 Overview — 3-level data 의 큰 그림
Ch.4 정규 종단 MRM — 2-level linear (3-level 의 토대)
§ 4.5 MRM 추정 — EM + Fisher scoring (3-level 에 직접 적용)
§ 9.7 NIMH 이항 분석 — 같은 NIMH 데이터의 이항 버전
§ 12.4.2 Empirical Bayes — Center ranking 의 표준

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

§ 13.2 ~ 13.2.4 — Nonlinear 3-level (probit, logistic, ordinal, nominal, count)
§ 13.3 — Ch.13 Summary
Ch.14 — 결측 데이터

관련 개념

Lorr & Klett (1966) — IMPS 척도 원전
Gibbons et al. (1988) — NIMH 분석 (3 약물 효과 비슷 발견)
Donner (1982) — Cluster ICC 의 일반 범위 (1-5%)
Donner & Klar (2000) — Cluster randomized trials
Murray (1998) — Group-randomized trials 표준
Jacobs et al. (1989), Siddiqui et al. (1996) — 추가 cluster 사례
Blair et al. (1983) — Type I error inflation by ignored ICC
Self & Liang (1987) — Boundary LR test
Raudenbush & Bryk (2002) — Multilevel models 표준
ICH E9 — Multi-center clinical trials 통계 표준