§ 13.2.4.2 ~ 13.2.4.3 — Three-Level Nominal 과 Count 의 깊이

§13.2.4.2 Nominal (Hedeker 2003): 식 (13.30) conditional likelihood · 식 (13.31-13.32) multinomial probabilities · 식 (13.33) 범주별 random effects · 식 (13.34-13.35) score · §13.2.4.3 Count (Siddiqui 1996): 식 (13.36) λ with offset · 식 (13.37-13.38) Poisson conditional · 식 (13.39-13.40) likelihood + score

Hedeker & Gibbons (2006) Ch.13 §13.2.4.2 + §13.2.4.3 의 깊이 있는 풀이. 13-4 sub-post 의 ordinal 확장에 이어 nominal + count 의 3-level 일반화. §13.2.4.2 Nominal (Hedeker 2003): reference cell formulation + multinomial logit. 식 (13.30) conditional likelihood = \(\prod \prod p_{ijc}^{d_{ijkc}}\). 식 (13.31-13.32) multinomial probabilities (reference cell \(c = 1\)). 식 (13.33) 결정적 차이 — ordinal 과 달리 모든 covariate 효과 \(\beta_c\) + random effects \(T_c, \sigma_{c(3)}\) 가 범주별 다름. 식 (13.34) score = \((d_{ijkc} - p_{ijkc}) \cdot \partial z / \partial \eta_c\) — Ch.11 식 (11.10) 의 3-level 일반화. §13.2.4.3 Count (Siddiqui 1996): 식 (13.36) \(\lambda_{ijk} = t_{ijk} \exp(z_{ijk})\)time/exposure offset \(t_{ijk}\) 포함. 식 (13.37) response model 동일 (probit/ordinal/nominal 와 같은 형태). 식 (13.38-13.39) Poisson conditional probability + likelihood. 식 (13.40) score = \(\sum (y_{ijk} - \lambda_{ijk}) \cdot \partial z / \partial \eta\) — Ch.12 식 (12.22) 의 3-level 일반화. § 11 / § 12 의 2-level 의 직접 확장.

Statistics
저자

Kwangmin Kim

공개

2026년 05월 06일

1 들어가며 — 13-4 의 자연 확장

§ 13.2.4 ~ 13.2.4.1 (13-4) 에서 일반 framework + ordinal 확장을 다뤘다. 본 sub-post 는 그 framework 의 nominal + count 응답 형태로의 자연 확장 — 응답 형태별 conditional likelihood \(\ell_{ij}\) 만 변경.

한 줄 요약

“§13.2.4.2 Nominal (Hedeker 2003) = ordinal (§ 13.2.4.1) 의 두 제약 (회귀 계수 c-무관 + random effect 분산 c-공통) 을 모두 풀어 모든 covariate + random effects 가 범주별 다름 (식 13.33). 식 (13.30) conditional likelihood + 식 (13.31-13.32) multinomial probabilities (reference cell). 식 (13.34) score \((d - p) \cdot \partial z / \partial \eta_c\) — Ch.11 식 (11.10) 의 3-level 일반화. §13.2.4.3 Count (Siddiqui 1996): 식 (13.36) \(\lambda_{ijk} = t_{ijk} \exp(z_{ijk})\)time/exposure offset \(t_{ijk}\) 포함, rate 모형. 식 (13.40) score \((y - \lambda) \cdot \partial z / \partial \eta\) — Ch.12 식 (12.22) 의 3-level 일반화. 두 모형 모두 § 13.2.4 의 통합 framework 직접 적용 — score 식만 응답별 다르고 적분 분해 + Fisher scoring 동일.”

2 § 13.2.4.2 — Nominal Outcomes

2.1 발전사 + 표기

Hedeker (2003) — Mixed-Effects Multinomial 정립

저자 본문 인용:

“As described in Hedeker [2003], the nominal mixed model is typically written in terms of the logistic representation. Also, it’s common to use the first category as a reference cell and to estimate the covariate effects, which vary across categories, relative to this reference cell.”

관행:

  • Logistic representation: probit 보다 logistic 흔함 (OR 해석).
  • Reference cell: 첫 범주 (\(c = 1\)) 가 reference, 다른 범주 (\(c = 2, \ldots, C\)) 가 비교 대상.
  • 범주 가변 효과: 모든 covariate 효과가 reference 대비 다름.
직관 — Ordinal vs Nominal 의 결정적 차이

Ordinal (§ 13.2.4.1):

  • 모든 cumulative logit 에 같은 회귀 계수 \(\beta\).
  • 모든 cluster 가 같은 random effects 분산.
  • → Proportional odds 가정.

Nominal (§ 13.2.4.2):

  • 각 contrast 마다 별도 회귀 계수 \(\beta_c\) (\(c = 2, \ldots, C\)).
  • 각 contrast 마다 별도 random effects 분산 (\(T_c, \sigma_{c(3)}\)).
  • → 모수 폭발 + 가장 일반적.

Ch.10 → Ch.11 의 평행 (2-level):

  • Ordinal Ch.10: \(\beta\) 한 셋 + \(\sigma\) 한 개.
  • Nominal Ch.11: \(\beta_c\) \(C-1\) 셋 + \(\sigma_c\) \(C-1\) 개.

§ 13.2.4.1 → § 13.2.4.2 의 평행 (3-level):

  • 같은 패턴, 다만 random effects 가 cluster + subject 두 수준.

모수 차원: Ordinal 의 \(C\) 배.

2.2 식 (13.30-13.32) — 모형 정의

식 (13.30) — Conditional likelihood

\[ \ell_{ij}(\theta) = \prod_{k=1}^{n_{ij}} \prod_{c=1}^C (p_{ijc})^{d_{ijkc}} \tag{13.30} \]

표기:

  • \(d_{ijkc} = 1\) if \(y_{ijk} = c\), 0 otherwise (각 관측에서 \(C\) 개 indicator 중 하나만 1).
  • \(p_{ijc}\): 관측 (\(i, j, k\)) 의 범주 \(c\) 응답 확률, random effects 조건부.
식 (13.31-13.32) — Multinomial probabilities

Reference cell (\(c = 1\)) 표기:

\[ p_{ijc} = \frac{\exp(z_{ijkc})}{1 + \sum_{h=2}^C \exp(z_{ijkh})} \quad (c = 2, 3, \ldots, C) \tag{13.31} \]

\[ p_{ij1} = \frac{1}{1 + \sum_{h=2}^C \exp(z_{ijkh})} \tag{13.32} \]

→ Ch.11 식 (11.1-11.2) 의 직접 확장.

직관 — Reference Cell 의 정당성

§ 11 의 식 (11.1-11.2) 와 동일한 발상:

  • \(c = 1\): \(\beta_1 = 0, T_1 = 0, \sigma_{1(3)} = 0\) 으로 fixed.
  • 다른 범주: reference 대비 효과.
  • 결과: 모수 식별 가능.

Reference 선택:

  • 보통 가장 흔한 범주 또는 “control” 범주.
  • 임상적 해석 용이.
  • 적합도, 예측은 reference 선택 무관.

3-Level 의 새로운 점:

  • \(\beta_c, T_c, \sigma_{c(3)}\) 모두 범주별.
  • 모수 폭발: \(p \times (C-1)\) 회귀 + \(r(r+1)/2 \times (C-1)\) subject random + \((C-1)\) cluster random.
  • 4 범주 + 5 covariate + 2 subject random effects = \(5 \times 3 + 3 \times 3 + 3 = 27\) 모수.

명목 3-level 모형의 추정 어려움 — 큰 표본 + 적은 random effects 권고.

2.3 식 (13.33) — Response Model

범주별 separate response model

\[ z_{ijkc} = \sigma_{c(3)} \theta_{0i} + z_{ijk}^\top T_c \theta_{ij} + x_{ijk}^\top \beta_c \tag{13.33} \]

각 contrast \(c\) 가 별도 회귀 + 별도 random effects:

  • \(\sigma_{c(3)}\): 범주 \(c\) 의 cluster random intercept SD.
  • \(T_c\): 범주 \(c\) 의 subject random effects Cholesky.
  • \(\beta_c\): 범주 \(c\) 의 fixed effects.
직관 — 범주별 잠재 효용 framework

§ 11 의 extremal concept (Bock 1972) 의 3-level 적용:

  • 각 범주 \(c\) 에 별도 잠재 효용 \(y_{ijkc}\).
  • 잠재 효용 = cluster effect + subject effect + fixed + 잡음.
  • 응답 = 최대 잠재 효용의 범주.

3-Level 의 random effects 차원:

  • Cluster: \(C - 1\) 개 (범주별).
  • Subject: \(C - 1\)\(T_c\) matrix.
  • 적분 차원이 폭발 — adaptive Gauss-Hermite 또는 MCMC 필수.

임상 응용 예시 (의료 이용 — 3 범주: no care / outpatient / inpatient):

  • \(\sigma_{2(3)}\): hospital 의 outpatient utilization 변동.
  • \(\sigma_{3(3)}\): hospital 의 inpatient utilization 변동.
  • 두 효과가 별도 — hospital 마다 다른 mix.

두 차원 hospital quality 분석:

  • Outpatient capacity vs inpatient capacity.
  • 정책 결정 — 각 service 별 자원 배분.

§ 11.3.1 (organ transplantation, 11-3) 의 3-level 확장 — region × hospital × patient 같은 구조.

2.4 식 (13.34-13.35) — Score

식 (13.34) — 범주별 score

\[ \frac{\partial \log \ell_{ij}(\theta)}{\partial \eta_c} = \sum_{k=1}^{n_{ij}} (d_{ijkc} - p_{ijkc}) \frac{\partial z_{ijkc}}{\partial \eta_c} \tag{13.34} \]

→ “잔차 × 모수 미분” 의 자연 형태.

식 (13.35) — 모수별 미분

\[ \frac{\partial z_{ijkc}}{\partial \beta_c} = x_{ijk}, \quad \frac{\partial z_{ijkc}}{\partial v(T_c)} = J_r (\theta \otimes z_{ijk}), \quad \frac{\partial z_{ijkc}}{\partial \sigma_{c(3)}} = \theta_{(3)} \tag{13.35} \]

직관 — Ch.11 식 (11.10) 의 3-Level 일반화

Ch.11 의 nominal 2-level score (식 11.10):

\[ \frac{\partial \log L}{\partial \Delta} = \sum_i h^{-1}(Y_i) \int [\sum_j D(y_{ij} - P_{ij}) \otimes \partial \Delta] \ell_i g(\theta) d\theta \]

§ 13.2.4.2 의 3-level score (식 13.34):

\[ \frac{\partial \log \ell_{ij}}{\partial \eta_c} = \sum_k (d_{ijkc} - p_{ijkc}) \frac{\partial z_{ijkc}}{\partial \eta_c} \]

같은 형태, 다른 차원:

  • 잔차 (관측 indicator - 예측 확률).
  • 모수 미분.
  • 곱의 합.

3-level 의 추가 적분 (식 13.22-13.23):

  • Subject score (식 13.34) 가 외부 + 내부 적분의 안에 들어감.
  • 외부: cluster \(\theta_{(3)}\).
  • 내부: subject \(\theta_{(2)}\).
  • → § 13.2.4 의 일반 framework 직접 적용.

모수 vector 의 차원:

  • 각 범주 \(c\) 마다 별도 score.
  • 총 모수: \(\sum_c (\dim(\beta_c) + \dim(v(T_c)) + 1)\).
  • 큰 차원에서 추정 어려움 — Bayesian (brms) 권고.

3 § 13.2.4.3 — Count Outcomes

3.1 식 (13.36-13.37) — 모형 정의

Siddiqui (1996) — 2-level Mixed-Effects Poisson 의 3-level 확장

저자 본문 인용:

“Siddiqui [1996] describes the statistical development of the 2-level mixed-effects Poisson regression model. Here, we will present the key computational features for the extension to three levels.”

식 (13.36) — Expected count with offset:

\[ \lambda_{ijk} = t_{ijk} \exp(z_{ijk}) \tag{13.36} \]

표기:

  • \(t_{ijk}\): 관측 (\(i, j, k\)) 의 follow-up time / exposure.
  • \(z_{ijk}\) = 식 (13.37) 의 response model.

식 (13.37) — Response model:

\[ z_{ijk} = \sigma_{(3)} \theta_{0i} + z_{ijk}^\top T \theta_{ij} + x_{ijk}^\top \beta \tag{13.37} \]

→ Probit (식 13.21) 와 동일 형태.

직관 — Time/Exposure Offset

식 (13.36) 의 \(t_{ijk}\) 의 결정적 역할:

  • \(\lambda\) 가 단위 시간 당 사건 수 (rate) 가 아니라 총 관측 동안 사건 수.
  • \(t_{ijk}\) 가 관측 \(i, j, k\) 의 노출 시간 (또는 인구).
  • → log link 형태: \(\log \lambda = \log t + z = \log t + \sigma_{(3)}\theta_{0i} + z'T\theta + x'\beta\).

Offset 의 자연성:

  • 의료 이용 횟수: \(t\) = follow-up days.
  • 자살률: \(t\) = population (county-year).
  • 사고: \(t\) = vehicle-miles traveled.
  • 결함: \(t\) = production volume.

없으면 (\(t_{ijk} = 1\) 가정):

  • 모든 관측이 같은 노출 — 단순 count 모형.

3-Level 의 의미:

  • Cluster level (예: hospital), subject level (예: department), observation level (예: month).
  • 각 수준의 random effects 가 unobserved heterogeneity 흡수.
  • \(t_{ijk}\) 가 관측별 노출 처리.

§ 12.5 (자살률) 의 framework 가 정확히 같은 구조 — 다만 § 12.5 는 2-level (county random effect 만), § 13.2.4.3 은 3-level (cluster + subject random effects).

3.2 식 (13.38-13.39) — Conditional Probability + Likelihood

식 (13.38) — Poisson conditional probability

\[ P(Y_{ijk} = y_{ijk} \mid \theta) = \exp(-\lambda_{ijk}) \frac{\lambda_{ijk}^{y_{ijk}}}{y_{ijk}!} \tag{13.38} \]

표준 Poisson density.

식 (13.39) — Conditional likelihood

\[ \ell_{ij}(\theta) = \prod_{k=1}^{n_{ij}} \exp(-\lambda_{ijk}) \frac{\lambda_{ijk}^{y_{ijk}}}{y_{ijk}!} \]

전개:

\[ = \exp\left[-\sum_k \lambda_{ijk} + \sum_k y_{ijk} \log t_{ijk} + \sum_k y_{ijk} z_{ijk} - \sum_k \log(y_{ijk}!)\right] \tag{13.39} \]

직관 — Poisson 곱의 정리

식 (13.39) 의 단계별 정리:

  1. Original 곱: \(\prod_k e^{-\lambda} \lambda^y / y!\).

  2. Log 변환: \(\sum_k [-\lambda_{ijk} + y_{ijk} \log \lambda_{ijk} - \log y_{ijk}!]\).

  3. \(\lambda_{ijk} = t_{ijk} e^{z_{ijk}}\) 대입:

\[ \log \lambda_{ijk} = \log t_{ijk} + z_{ijk} \]

  1. 정리:

\[ \sum_k [-t_{ijk} e^{z_{ijk}} + y_{ijk} (\log t_{ijk} + z_{ijk}) - \log y_{ijk}!] \]

  1. 재배열 (식 13.39 형태):

\[ -\sum_k \lambda_{ijk} + \sum_k y_{ijk} \log t_{ijk} + \sum_k y_{ijk} z_{ijk} - \sum_k \log y_{ijk}! \]

관찰:

  • \(\log t_{ijk}\)\(\log y_{ijk}!\) 는 모수 무관 (constants for given data).
  • → 추정에는 \(-\sum \lambda + \sum y_{ijk} z_{ijk}\) 만 영향.

§ 12.4.1 (12-3) 의 식 (12.21) 과 동일 형태 — 2-level Poisson 의 직접 확장.

3.3 식 (13.40) — Score

식 (13.40) — Score

\[ \frac{\partial \log \ell_{ij}(\theta)}{\partial \eta} = \sum_{k=1}^{n_{ij}} (y_{ijk} - \lambda_{ijk}) \frac{\partial z_{ijk}}{\partial \eta} \tag{13.40} \]

with:

\[ \frac{\partial z_{ijk}}{\partial \beta} = x_{ijk}, \quad \frac{\partial z_{ijk}}{\partial v(T)} = J_r (\theta_{(2)} \otimes z_{ijk}), \quad \frac{\partial z_{ijk}}{\partial \sigma_{(3)}} = \theta_{(3)} \]

직관 — Ch.12 식 (12.22) 의 3-Level 일반화

Ch.12 의 2-level Poisson score (식 12.22):

\[ \frac{\partial \log \ell}{\partial \beta} = \sum_j (y_{ij} - \lambda_{ij}) x_{ij} \]

§ 13.2.4.3 의 3-level score (식 13.40):

\[ \frac{\partial \log \ell_{ij}}{\partial \eta} = \sum_k (y_{ijk} - \lambda_{ijk}) \frac{\partial z_{ijk}}{\partial \eta} \]

같은 “잔차 × 공변량” 형태, 차이는:

  1. 합의 범위: 2-level 은 같은 환자의 시점 합, 3-level 은 같은 subject 의 관측 합.
  2. 모수 미분: 3-level 에 cluster random effect \(\theta_{(3)}\) 추가.
  3. 적분: 식 (13.22-13.23) 의 외부 + 내부 적분 안에 들어감.

Score 의 핵심 항 \((y_{ijk} - \lambda_{ijk})\):

  • 잔차 = 관측 - 예측.
  • 표준 GLM 의 자연 형태 (canonical link).
  • Poisson + log link 의 가장 단순한 score.

Hessian (식 13.40 의 미분):

\[ \frac{\partial^2 \log \ell}{\partial \eta \partial \eta^\top} = -\sum_k \lambda_{ijk} \frac{\partial z_{ijk}}{\partial \eta} \frac{\partial z_{ijk}}{\partial \eta^\top} \]

→ 음정치 (Newton-Raphson 안정).

§ 12 의 toolkit 그대로 적용 — adaptive Gauss-Hermite + Fisher scoring.

4 § 11 / § 12 와의 비교 (2-Level → 3-Level)

통합 framework 의 평행
응답 2-Level (Ch.) 3-Level (§ 13.2.4)
Binary Ch.9 (probit + logit) § 13.2.1-13.2.2
Ordinal Ch.10 (cumulative) § 13.2.4.1
Nominal Ch.11 (multinomial) § 13.2.4.2
Count Ch.12 (Poisson + ZIP) § 13.2.4.3
직관 — 4 가지 일반화의 공통 구조

모든 일반화가 같은 패턴:

  1. Conditional likelihood \(\ell_{ij}\): 응답 형태별 다름.
    • Binary: Bernoulli 곱.
    • Ordinal: Cumulative 차이의 곱.
    • Nominal: Multinomial 곱.
    • Count: Poisson 곱.
  2. Score 구조: 모두 “잔차 × 모수 미분” 형태.
    • Binary: \((y - \Phi(z)) \cdot \partial z / \partial \eta\) (with weight).
    • Ordinal: Cumulative 차이의 미분 (with Kronecker delta for thresholds).
    • Nominal: \((d - p) \cdot \partial z / \partial \eta_c\) (범주별).
    • Count: \((y - \lambda) \cdot \partial z / \partial \eta\).
  3. Random effects 구조 (식 13.21 또는 13.33 형태):
    • \(z = \sigma_{(3)}\theta_{0i} + z'T\theta_{ij} + x'\beta\) (or \(\beta_c\)).
    • 3-level: cluster + subject + observation.
  4. 적분 분해 (식 13.17-13.19):
    • 모든 응답에 동일.
    • 차원이 \(r + 1\) 으로 폭발 회피.
  5. Fisher scoring (식 13.24-13.25):
    • 모든 응답에 동일.
    • BHHH-like empirical information.

§ 13.2.4 의 통합 framework 가 모든 응답 형태의 일반화를 지원.

임상 응용 분야 매핑:

응답 3-Level 응용
Binary Multi-center 임상시험 (TVSFP, 13-3)
Ordinal School × Class × Student 학업 등급
Nominal Hospital × Department × Patient 진료 형태
Count Region × County × Year 의 자살률 (12-5 의 3-level 확장)

5 응용 분야

분야 응답 형태 3-Level 구조
학업 등급 (3+ 범주) Ordinal School × Class × Student
의료 이용 형태 Nominal Hospital × Department × Patient
통근 수단 Nominal Region × Family × Member
County × Year 자살 Count State × County × Year
Hospital 결함 횟수 Count System × Hospital × Month
Patient × Day 의료 사고 Count Hospital × Patient × Day
Family × Member × Year 의료 이용 Count Region × Family × Member

→ “Three-level + categorical/count” 의 모든 사례.

6 코드 예시

6.1 Step 1: 3-Level Nominal 시뮬레이션

library(mclogit)
library(dplyr)


# 3-level nominal 시뮬레이션 (3 범주, reference = 1)
set.seed(2026)
n_clusters <- 15
n_subjects_per_cluster <- 8
n_times_per_subject <- 4

# 모수 (3 범주 = reference + 2 contrast)
beta_2 <- c(0.5, 0.3)  # 범주 2 의 회귀 계수 (intercept, x)
beta_3 <- c(-0.5, 0.5)  # 범주 3 의 회귀 계수
sigma_2_cluster <- 0.4
sigma_3_cluster <- 0.6  # 범주 3 의 cluster 변동 더 큼
sigma_2_subject <- 0.5
sigma_3_subject <- 0.7

# Random effects (각 범주별)
cluster_effects_2 <- rnorm(n_clusters, 0, sigma_2_cluster)
cluster_effects_3 <- rnorm(n_clusters, 0, sigma_3_cluster)

df <- data.frame()
subject_id <- 0
for (i in 1:n_clusters) {
  for (j in 1:n_subjects_per_cluster) {
    subject_id <- subject_id + 1
    subj_eff_2 <- rnorm(1, 0, sigma_2_subject)
    subj_eff_3 <- rnorm(1, 0, sigma_3_subject)

    for (k in 1:n_times_per_subject) {
      x <- rnorm(1)
      # 식 (13.33) — z_ijkc for c = 2, 3 (reference = 1)
      z_2 <- beta_2[1] + beta_2[2] * x + cluster_effects_2[i] + subj_eff_2
      z_3 <- beta_3[1] + beta_3[2] * x + cluster_effects_3[i] + subj_eff_3

      # 식 (13.31-13.32) — Probabilities
      denom <- 1 + exp(z_2) + exp(z_3)
      p_1 <- 1 / denom
      p_2 <- exp(z_2) / denom
      p_3 <- exp(z_3) / denom

      y <- sample(1:3, 1, prob = c(p_1, p_2, p_3))
      df <- rbind(df, data.frame(
        cluster = i, subject = subject_id,
        x = x, y = factor(y)
      ))
    }
  }
}

cat("전체 관측:", nrow(df), "\n")
cat("응답 분포:\n")
print(table(df$y))
시뮬레이션의 검증
  • 3 범주 응답 분포 확인 (각 범주가 합리적 비율).
  • Cluster + subject random effects 가 두 contrast (c=2, c=3) 에 별도 영향.
  • → 식 (13.33) 의 정확한 generative process.

6.2 Step 2: 3-Level Nominal 적합 (R mclogit::mblogit)

# 3-level nominal — mblogit 가 random effects 지원
fit_3level_nominal <- mblogit(y ~ x,
                               random = ~1 | subject,  # subject random effect
                               data = df)
summary(fit_3level_nominal)

# Cluster random effect 직접 지원 약함 — 단순화 또는 brms 우회
# Bayesian 우회 권고

# Bayesian 적합 (brms)
library(brms)
fit_brm_nominal <- brm(y ~ x + (1 | subject) + (1 | cluster),
                       data = df, family = categorical(),
                       chains = 2, iter = 2000, refresh = 0)
summary(fit_brm_nominal)
R 의 3-Level Nominal 적합

Frequentist 한계:

  • mclogit::mblogit: random intercept + slope 지원, 3-level 직접 어려움.
  • nnet::multinom: random effects 무.

Bayesian 권고:

  • brms::brm(family = categorical()): 3-level 직접 지원.
  • 각 contrast 별 random effects 추정.
  • Posterior predictive checks.

대안:

  • SAS PROC NLMIXED: 직접 코딩.
  • Hedeker MIXNO: nominal mixed-effects 표준 도구 (단, 옛 software).

현대 표준 (2026):

  • brms 또는 Stan 직접.
  • \(C-1\) 개 contrast 별 random effects → 모수 폭발 → strong prior 또는 simplification.

6.3 Step 3: 3-Level Count 시뮬레이션

library(lme4)


# 3-level Poisson 시뮬레이션 (식 13.36)
set.seed(2026)
n_regions <- 10  # cluster (level 3)
n_counties_per_region <- 8  # subject (level 2)
n_years <- 5  # observation (level 1)

# 모수
beta_intercept <- 0.5
beta_drug <- -0.3  # antidepressant 효과
sigma_region <- 0.3
sigma_county <- 0.4

# Random effects
region_effects <- rnorm(n_regions, 0, sigma_region)

df <- data.frame()
county_id <- 0
for (r in 1:n_regions) {
  for (c in 1:n_counties_per_region) {
    county_id <- county_id + 1
    county_effect <- rnorm(1, 0, sigma_county)
    population <- round(runif(1, 10000, 1000000))  # county population

    for (y in 1:n_years) {
      drug_use <- rnorm(1)  # antidepressant log-pills per person
      # 식 (13.36-13.37): lambda = t * exp(z), z = beta + random
      z <- beta_intercept + beta_drug * drug_use +
           region_effects[r] + county_effect
      lambda <- population * exp(z) / 100000  # rate per 100,000 * population
      count <- rpois(1, lambda)
      df <- rbind(df, data.frame(
        region = r, county = county_id, year = y,
        drug_use = drug_use, population = population,
        count = count
      ))
    }
  }
}

cat("전체 관측:", nrow(df), "\n")
cat("Counties:", length(unique(df$county)), "\n")
cat("Regions:", length(unique(df$region)), "\n")
자살률 데이터 시뮬레이션

3-level 자살률 데이터의 자연 구조:

  • Level 3: Region (state).
  • Level 2: County within region.
  • Level 1: County-year (시간).
  • Outcome: Suicide count.
  • Offset: Population.

§ 12.5 (Gibbons 2005) 의 framework + 추가 region clustering.

\(t_{ijk}\) (offset):

  • \(\log t = \log(\text{population}/100000)\) — log per 100,000 rate.
  • 회귀 식: \(\log \lambda = \log t + z\).

6.4 Step 4: 3-Level Poisson 적합 (R lme4::glmer)

# 3-level Poisson with offset
fit_3level_poisson <- glmer(count ~ drug_use + offset(log(population/100000)) +
                            (1 | region) + (1 | county),
                            data = df, family = poisson,
                            control = glmerControl(optimizer = "bobyqa"))
summary(fit_3level_poisson)

# 추정값 비교
cat("\n추정값 (진짜 vs 추정):\n")
cat("  beta_intercept (진짜 0.5):",
    round(fixef(fit_3level_poisson)["(Intercept)"], 3), "\n")
cat("  beta_drug (진짜 -0.3):",
    round(fixef(fit_3level_poisson)["drug_use"], 3), "\n")

# Random effects
varcomp <- VarCorr(fit_3level_poisson)
cat("\nRandom effects SD:\n")
cat("  Region SD (진짜 0.3):",
    round(sqrt(as.numeric(varcomp$region)), 3), "\n")
cat("  County SD (진짜 0.4):",
    round(sqrt(as.numeric(varcomp$county)), 3), "\n")

# Empirical Bayes — county 별 random effects (자살률 ranking 용도)
ranef_county <- ranef(fit_3level_poisson)$county
cat("\n10 outlier counties (highest random effect):\n")
top_10 <- order(ranef_county[, 1], decreasing = TRUE)[1:10]
print(round(ranef_county[top_10, 1], 3))
정책 활용 — County Ranking

3-level Poisson 의 EB random effects:

  • County 별 자살률 추정 (covariate 보정 후).
  • Ranking → 정책 우선순위.

§ 12.5 (자살률, 2-level) vs § 13.2.4.3 (3-level) 의 차이:

  • 2-level: county random effect 만 (state clustering 무시).
  • 3-level: state random effect 추가 → state 효과 정확 분리.

State 효과의 임상적 의미:

  • 같은 state 의 county 들이 비슷한 정책, 의료 시스템.
  • → State clustering 무시 시 county random effect 추정 부정확.
  • 정확한 ranking 위해 3-level 권고.

Convergence:

  • 3-level Poisson 은 큰 표본 + 적은 random effects 면 안정.
  • Convergence 어려우면: nAGQ = 1 (default), bobyqa optimizer.
  • 또는 glmmTMB (더 안정) 또는 brms (Bayesian).

7 관련 주제

선행 지식

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

  • § 13.3 — Ch.13 Summary
  • Ch.14 — 결측 데이터

관련 개념

  • Hedeker (2003) — Mixed-Effects multinomial 정립
  • Siddiqui (1996) — 2-level mixed-effects Poisson 의 토대
  • Bock (1972) — Multinomial logit + extremal concept
  • Magnus (1988) — Elimination matrix
  • Stroud & Sechrest (1966) — Gauss-Hermite quadrature
  • McCullagh & Nelder (1989) — GLM canonical link
  • Ch.11.3.1 Organ Transplant — Nominal 의 competing risks 응용
  • § 12.5 Suicide Rates — Count 의 정책 응용

Subscribe

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