1 들어가며 — 랜덤 절편이 추가되면 무엇이 달라지는가

Ch.9 Overview 에서 GLMM 이항의 전체 구도를, § 9.2-9.3 에서 단일 수준 로지스틱·프로빗을, § 9.4 에서 잠재 변수 framework 를 다뤘다. 본 포스트는 그 재료들을 결합해 mixed-effects 이항 모형의 기본 구조 를 완성한다.

§ 9.5 본문은 식 (9.13) 의 랜덤 절편 추가에서 출발해 § 9.5.5 까지 다섯 개 소절로 펼쳐진다. 본 포스트는 그중 모형 정의 (9.5 intro) + ICC (9.5.1) + 다중 랜덤 효과 (9.5.2) + 이질 분산 (9.5.3) 의 네 단계를 다룬다. 다수준 표현 (9.5.4) 과 response function 형태 (9.5.5) 는 뒤이은 sub-post 에서 별도로 정리한다.

한 줄 요약

“랜덤 절편 한 개를 로지스틱에 추가하면 잠재 변수 분산이 \(\pi^2/3\) 에서 \(\sigma_v^2 + \pi^2/3\) 로 늘어난다 (9.5 intro). 이 분산 비율이 ICC 다 (9.5.1). 절편이 여러 개로 늘어나면 Cholesky 로 안정 추정 (9.5.2). 그룹별로 다른 분산이 필요하면 design vector 를 dummy code 로 바꾸면 끝 — 쌍둥이 tetrachoric, IRT 2PL 까지 같은 틀 (9.5.3).”

2 § 9.5 — Mixed-Effects Logistic 모형 정의

2.1 왜 랜덤 효과가 필요한가

군집 (clustered) 이항 데이터의 문제

Single-level 로지스틱 (§ 9.2) 은 관측이 독립 임을 가정한다. 그러나 종단 데이터에서:

같은 환자의 반복 측정 → 같은 환자 안에서 양의 상관
같은 클리닉/학교의 환자들 → 같은 클러스터 안에서 양의 상관

독립 가정이 깨지면 standard error 가 과소평가되고 검정의 type I error 가 부풀려진다. 추정 자체는 일치 (consistent) 라 해도, 추론은 잘못된다.

직관 — “환자별 응답 성향” 의 모형화

같은 처치를 받아도 어떤 환자는 호전 잘하고 어떤 환자는 잘 안 한다. 측정 못한 개인 특성 (유전, 생활 습관, 동기 부여 등) 이 응답 확률에 일관되게 영향을 미친다. 이 개인별 잠재 반응 성향 을 모형에 명시하면:

같은 환자의 반복 측정이 자동으로 상관된다 (같은 \(\upsilon_i\) 를 공유하므로).
처치 효과 \(\beta\) 는 “이 개인 특성을 통제한 상태에서” 의 효과로 해석된다 (subject-specific).

이 발상이 § 4 의 정규 종단 MRM 과 정확히 동일하다. § 9 의 차이는 link 함수가 비선형 (logit) 이라는 것뿐이다. 비선형이 만드는 결과 — subject-specific vs marginal 의 비축소성 (식 9.16) — 가 § 9.5 의 가장 중요한 통찰이다.

2.2 식 (9.13) — 랜덤 절편 로지스틱의 출발

표기

\(i = 1, \ldots, N\): level-2 단위 (피험자, 클러스터)
\(j = 1, \ldots, n_i\): level-1 단위 (피험자 내 반복 측정)
\(Y_{ij} \in \{0, 1\}\): \(i\) 번째 피험자의 \(j\) 번째 이항 반응
\(p_{ij} = P(Y_{ij} = 1)\)
\(x_{ij}\): \((p+1) \times 1\) 공변량 벡터 (절편 1 포함)
\(\beta\): \((p+1) \times 1\) 회귀 계수 벡터
\(\upsilon_i\): 피험자 \(i\) 의 랜덤 절편, \(\upsilon_i \sim \mathcal{N}(0, \sigma_v^2)\)

식 (9.13)

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = x_{ij}^\top \beta + \upsilon_i \tag{9.13} \]

단일 수준 로지스틱 (식 9.1) 에 피험자별 절편 \(\upsilon_i\) 한 개를 더한 형태.

직관 — “log-odds 가 환자별로 평행 이동”

식 (9.13) 의 우변에서 \(\upsilon_i\) 는 처치 효과 \(\beta\) 와 무관한 환자 고유의 절편 이동. 그래프로 그리면:

\(\upsilon_i > 0\) 인 환자는 모든 시점·모든 공변량 값에서 log-odds 가 위로 평행 이동 → 응답 확률이 일관되게 높다.
\(\upsilon_i < 0\) 인 환자는 일관되게 낮다.

기울기 \(\beta\) 는 환자 사이에서 동일 — 처치 효과는 누구에게나 같은 방향, 다만 출발점이 다르다.

§ 9.5.2 에서 기울기까지 환자별로 다르게 만드는 모형 (랜덤 기울기) 으로 확장되지만, 출발은 절편 하나가 가장 자연스럽다.

2.3 식 (9.14) — 표준화 형태

식 (9.14)

\(\upsilon_i = \sigma_v \theta_i\) 로 변환 (\(\theta_i \sim \mathcal{N}(0, 1)\)):

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = x_{ij}^\top \beta + \sigma_v \theta_i \tag{9.14} \]

왜 표준화하나 — 추정과 해석의 두 이점

추정 측면: Gauss-Hermite 구적법 (§ 9.6) 은 표준 정규 \(\theta_i\) 에 대한 적분을 미리 계산된 노드·가중치로 근사한다. 모형이 \(\sigma_v\) 까지 직접 분리해 두면 적분 변수가 항상 표준 정규로 고정되어 알고리즘이 단순해진다.

해석 측면: \(\sigma_v\) 가 회귀식에 명시적으로 등장한다. 즉 \(\beta\) 와 같은 척도 (log-odds) 위에 있다. 예를 들어 \(\sigma_v = 0.7\) 이면 “환자 간 표준편차가 log-odds 0.7 만큼 — 평균 환자에서 ±1 SD 환자로 옮겨가면 odds 가 \(e^{0.7} \approx 2\) 배” 식으로 직접 읽을 수 있다.

2.4 식 (9.15) — 잠재 변수 표현

식 (9.15)

§ 9.4 의 잠재 변수 framework 로 다시 쓰면:

\[ y_{ij} = x_{ij}^\top \beta + \sigma_v \theta_i + \epsilon_{ij} \tag{9.15} \]

\(Y_{ij} = 1 \iff y_{ij} > 0\), \(\epsilon_{ij}\) 의 분포가 logistic 이면 logit, 정규면 probit.

핵심 통찰 — 잠재 분산이 두 조각으로

식 (9.15) 의 우변에서 무작위 부분은 \(\sigma_v \theta_i + \epsilon_{ij}\) 두 개. 따라서 잠재 변수의 조건부 분산 (공변량 \(x\) 를 고정한 상태) 은:

\[ V(y_{ij} \mid x_{ij}) = \sigma_v^2 + \sigma_\epsilon^2 \]

여기서:

\(\sigma_v^2\): 피험자 간 (between-subject) 분산. 환자 고유 성향의 이질성.
\(\sigma_\epsilon^2\): 피험자 내 (within-subject) 잔차 분산. logistic 이면 \(\pi^2/3 \approx 3.29\), probit 이면 1 로 고정 (잠재 분산의 식별 불가능성 때문, § 9.4 참조).

Single-level 모형 (식 9.11) 에서는 잠재 분산이 \(\sigma_\epsilon^2\) 한 조각이었다. 랜덤 효과를 추가하면 그 위에 \(\sigma_v^2\) 한 조각이 더 쌓인다. 이 단순한 사실 — 분산이 두 조각으로 분해된다는 것 — 가 § 9.5 의 거의 모든 결과를 만들어낸다.

2.5 식 (9.16) — Subject-Specific vs Marginal 의 비축소성

식 (9.16) — 회귀 계수 scale 차이

같은 데이터에 fixed-effects (또는 GEE) 와 mixed-effects 를 적합하면 회귀 계수가 다음 비율로 차이난다:

\[ \beta_M \approx \sqrt{\frac{\sigma_v^2 + \sigma_\epsilon^2}{\sigma_\epsilon^2}} \, \beta_F \tag{9.16} \]

\(\beta_F\): fixed-effects (또는 GEE marginal) 계수, \(\beta_M\): mixed-effects (subject-specific) 계수.

직관 — “분산이 커지면 계수도 커져야 같은 확률”

같은 응답 확률 \(P(Y = 1) = 0.7\) 을 만들려면, 잠재 변수 분산이 클수록 평균값을 더 위로 끌어올려야 한다 (분산이 크면 분포가 더 넓게 퍼져 임계값을 넘기기 어려워지므로).

수식으로 보면, 잠재 변수가 분산 \(V\) 를 가지면 \(P(Y = 1) = P(y > 0) = \Phi(\mu / \sqrt{V})\) (probit). 같은 확률을 유지하려면 \(\mu / \sqrt{V}\) 가 일정해야 하므로:

\[ \frac{\beta_M}{\sqrt{\sigma_v^2 + \sigma_\epsilon^2}} = \frac{\beta_F}{\sqrt{\sigma_\epsilon^2}} \]

→ \(\beta_M = \sqrt{(\sigma_v^2 + \sigma_\epsilon^2)/\sigma_\epsilon^2} \, \beta_F\).

비선형 link 의 본질적 결과 — log 의 평균이 평균의 log 와 같지 않다는 사실 (Jensen 부등식의 변종) 의 회귀 계수 버전. 정규 (정체) 모형에서는 \(V\) 가 어떤 값이든 회귀 계수가 변하지 않으므로 이 차이가 존재하지 않는다.

Zeger 보정 — 실용적 권고

Zeger et al. (1988) 는 logistic 모형에서 \(\sigma_\epsilon^2 = \pi^2/3\) 대신 \((15/16)^2 \pi^2 / 3\) 를 쓰는 것이 두 추정값을 일치시키는 데 더 정확함을 시뮬레이션으로 보였다. 이론적 정확값과 경험적 보정값이 약간 다른 이유는 logistic cdf 와 정규 cdf 의 꼬리 모양 차이.

실무에서 GEE 결과와 GLMM 결과를 비교할 때 이 보정 인자를 곱해 두면 둘이 잘 일치하는 것을 확인할 수 있다.

3 § 9.5.1 — Intraclass Correlation (ICC)

3.1 정의 — 잠재 분산의 분해

ICC 공식 — Logistic 모형

\[ \widehat{\rho}_v = \frac{\widehat{\sigma}_v^2}{\widehat{\sigma}_v^2 + \pi^2 / 3} \]

Probit 모형은 분모의 \(\pi^2/3\) 을 \(1\) 로 대체:

\[ \widehat{\rho}_v = \frac{\widehat{\sigma}_v^2}{\widehat{\sigma}_v^2 + 1} \]

직관 — “설명 안 된 분산 중 피험자 간이 차지하는 비율”

ICC 는 잠재 반응 성향의 총 분산 중 피험자 간 분산이 차지하는 비율. 식 (9.15) 의 분산 분해 \(V(y) = \sigma_v^2 + \sigma_\epsilon^2\) 에서 첫 조각의 비중.

해석:

\(\rho_v = 0\): 피험자 간 차이가 없음. 모든 변동은 피험자 내 잔차. → 단일 수준 로지스틱이면 충분.
\(\rho_v = 1\): 피험자 내 변동이 없음. 같은 환자의 반복 측정이 완벽 상관. → 종단 데이터로 다룰 의미가 없음 (모든 측정이 동일).
중간값 (예: 0.3-0.5): 의료 연구에서 흔한 범위. 환자 차이가 무시할 수 없으나 처치 효과 추정의 여지도 있음.

ICC 는 단순히 “랜덤 효과를 넣어야 하는가” 를 정량화하는 진단량 — ICC 가 0 에 가까우면 GLMM 과 single-level 의 결과가 거의 같지만, ICC 가 0.2 이상이면 두 결과가 크게 다르다.

3.2 왜 분모가 \(\pi^2 / 3\) 인가 — Logistic 분포의 분산

표준 logistic 의 분산

표준 logistic 분포 (location 0, scale 1) 의 분산은 \(\pi^2 / 3 \approx 3.29\). 유도:

\[ V(\epsilon) = \int_{-\infty}^{\infty} z^2 \cdot \frac{e^z}{(1 + e^z)^2} \, dz = \frac{\pi^2}{3} \]

비슷하게 표준 정규는 분산 1, 표준 Gumbel 은 \(\pi^2 / 6\) ≈ 1.645 (complementary log-log). 이 차이가 § 9.5.5 에서 link 함수별 모형의 scale 차이를 결정한다.

핵심 — 잠재 분산은 식별 불가능 (§ 9.4) 이라 어떤 값으로 고정해야 한다. 관행적으로 logistic 은 \(\pi^2/3\), probit 은 \(1\) 을 쓴다. 이 고정 때문에 ICC 분모가 link 별로 다른 상수로 결정된다.

3.3 실무 해석 예시

예시 — 정신과 임상시험

HAMD 우울증 척도가 7 이하인지 (관해 여부) 를 5 회 반복 측정. 랜덤 절편 로지스틱 적합 결과 \(\widehat{\sigma}_v = 1.4\).

\[ \widehat{\rho}_v = \frac{1.4^2}{1.4^2 + \pi^2 / 3} = \frac{1.96}{1.96 + 3.29} \approx 0.37 \]

해석: “관해 여부의 환자 간 변동이 잠재 반응 성향 총 변동의 약 37 % — 즉 처치 효과를 정확히 보려면 환자 차이를 통제해야 함. GEE 마진 모형으로는 이 차이가 보이지 않음.”

만약 같은 모형을 GEE 로 적합하면, 처치 효과 계수가 약 \(1/\sqrt{1 + 1.96/3.29} = 1/\sqrt{1.60} \approx 0.79\) 배 작게 나온다 (식 9.16 의 역방향) — subject-specific 효과의 약 80 % 가 marginal 효과.

4 § 9.5.2 — More General Mixed-Effects Models

4.1 식 (9.24) — 다중 랜덤 효과로의 확장

다중 랜덤 효과의 표기

\(z_{ij}\): \(r \times 1\) 랜덤 효과 design 벡터 (절편 1 + 시간 + 기타 시간 가변 공변량 등). 보통 \(z_{ij}\) 의 첫 원소는 1 (랜덤 절편).
\(\upsilon_i \sim \mathcal{N}_r(0, \Sigma_v)\): \(r\) 차원 다변량 정규, 평균 0, 분산-공분산 행렬 \(\Sigma_v\).
\(T\): \(TT^\top = \Sigma_v\) 인 \(r \times r\) 하삼각 행렬 (Cholesky 분해).
\(\upsilon_i = T \theta_i\) 로 표준화 (\(\theta_i \sim \mathcal{N}_r(0, I)\)).

식 (9.24)

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = x_{ij}^\top \beta + z_{ij}^\top T \theta_i \tag{9.24} \]

직관 — “랜덤 절편 + 랜덤 기울기” 의 자연스러운 일반화

단일 랜덤 절편 (식 9.14) 은 \(r = 1\), \(z_{ij} = 1\) (스칼라), \(T = \sigma_v\) 의 특수 경우. 랜덤 기울기를 추가하려면 \(r = 2\), \(z_{ij} = (1, t_{ij})^\top\), 그리고 \(T\) 는 \(2 \times 2\) 하삼각 행렬:

\[ T = \begin{bmatrix} \sigma_{v_0} & 0 \\ \tau_{10} & \tau_{11} \end{bmatrix} \]

이때 \(TT^\top\) 가 분산-공분산 행렬 \(\Sigma_v\):

\[ \Sigma_v = \begin{bmatrix} \sigma_{v_0}^2 & \sigma_{v_0} \tau_{10} \\ \sigma_{v_0} \tau_{10} & \tau_{10}^2 + \tau_{11}^2 \end{bmatrix} \]

대각: 절편·기울기 분산. 비대각: 절편·기울기 공분산.

→ 환자별로 출발점 \(\upsilon_{0i}\) 와 변화 속도 \(\upsilon_{1i}\) 가 다른 모형. 정규 MRM (§ 4.3) 과 구조는 동일, 차이는 link 가 비선형이라는 점뿐이다.

4.2 왜 Cholesky 분해를 쓰나

직관 — 행렬의 “제곱근” 으로 양정치성 보장

\(\Sigma_v\) 를 직접 추정하면 두 가지 문제가 생긴다:

문제 1 — 양정치성 위반: 분산-공분산 행렬은 양정치 (positive definite) 여야 한다. 그런데 자유 추정으로 얻은 행렬이 비양정치가 되어 모형 자체가 무효화될 수 있다. 특히 분산 추정값이 음수에 가까운 경우.

문제 2 — Boundary 문제: 분산이 0 에 가까울 때 (\(\sigma_v^2 \to 0\)), 우도 표면이 가장자리에서 평평해져 수렴이 느리고 표준오차 추정이 불안정.

Cholesky 분해는 두 문제를 동시에 해결:

양정치성 자동 보장: \(T\) 가 어떤 값이든 \(TT^\top\) 는 항상 양반정치 (positive semi-definite). 대각이 양수이면 양정치.
Boundary 처리: \(T\) 의 대각 원소가 0 에 가까워져도 (즉 분산이 0 에 가까워져도) 우도 표면에서 매끄럽게 다룰 수 있다. 분산이 진짜 0 일 때 수치적으로 안정.

비유하자면 행렬의 “제곱근” — 음수가 안 되도록 제곱근을 추정하면 제곱은 자동으로 비음수가 된다. 이 발상의 행렬 버전이 Cholesky.

Cholesky 분해의 행렬 형태

\(\Sigma_v\) 가 \(r \times r\) 양정치 대칭 행렬이면 유일한 하삼각 \(T\) 가 존재해 \(\Sigma_v = TT^\top\):

\[ \begin{bmatrix} \sigma_{11} & \cdots & \sigma_{1r} \\ \vdots & \ddots & \vdots \\ \sigma_{r1} & \cdots & \sigma_{rr} \end{bmatrix} = \begin{bmatrix} t_{11} & 0 & 0 \\ \vdots & \ddots & 0 \\ t_{r1} & \cdots & t_{rr} \end{bmatrix} \begin{bmatrix} t_{11} & \cdots & t_{r1} \\ 0 & \ddots & \vdots \\ 0 & 0 & t_{rr} \end{bmatrix} \]

원소 수: \(\Sigma_v\) 의 자유 모수 \(r(r+1)/2\) 개 = \(T\) 의 자유 모수 \(r(r+1)/2\) 개. 1:1 대응이라 정보 손실 없음.

추정 절차:

데이터에 대해 우도를 \(T\) 의 함수로 표현 (\(\Sigma_v\) 대신).
\(T\) 의 원소를 ML/REML 로 추정.
필요하면 \(\widehat\Sigma_v = \widehat{T} \widehat{T}^\top\) 로 환원.

5 § 9.5.3 — Heterogeneous Variance Terms

5.1 핵심 — Design Vector 만 바꾸면 그룹별 분산이 만들어진다

§ 9.5.2 와 § 9.5.3 의 차이 — 같은 식 (9.24) 의 두 해석

§ 9.5.2 의 다중 랜덤 효과는 여러 종류의 랜덤 효과 (절편 + 기울기 + …) 를 같은 피험자에게 부여한다. \(z_{ij}\) 의 원소는 1, \(t_{ij}\), \(t_{ij}^2\) 등 시간/공변량 함수.

§ 9.5.3 의 이질 분산은 하나의 랜덤 효과 (절편) 인데 그룹별로 다른 분산 을 허용한다. \(z_{ij}\) 의 원소는 그룹 dummy code (예: MZ/DZ 표시 변수). 모형 구조 (식 9.24) 는 같지만 design vector 의 의미가 다르다.

이 발상이 강력한 이유 — 구현 측면에서 새로운 코드를 짤 필요 없이 같은 GLMM 알고리즘을 그대로 쓸 수 있음. design vector 만 바꾸면 끝. § 9.5.3 의 두 응용 (쌍둥이, IRT) 모두 이 trick 의 사례.

5.2 응용 1 — 쌍둥이 연구 (식 9.25)

쌍둥이 데이터 구조

\(i = 1, \ldots, N\): 쌍둥이 쌍 (level-2 unit).
\(j = 1, 2\): 쌍 안의 두 개체 (level-1 unit). \(j\) 는 단지 두 개체 구분용.
\(MZ_i\), \(DZ_i\): 쌍이 일란성 (MZ) 인지 이란성 (DZ) 인지의 dummy code. (\(MZ_i = 1 \Rightarrow DZ_i = 0\), 반대도 같음.)

design vector \(z_{ij} = (MZ_i, DZ_i)^\top\).

식 (9.25)

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = x_{ij}^\top \beta + \begin{bmatrix} MZ_i & DZ_i \end{bmatrix} \begin{bmatrix} \sigma_{v(MZ)} \\ \sigma_{v(DZ)} \end{bmatrix} \theta_i \tag{9.25} \]

\(\theta_i \sim \mathcal{N}(0, 1)\) — 단일 표준 정규 효과. 그러나 dummy code 와 \(T\) 벡터의 곱이 그룹별로 다른 분산을 만들어낸다.

직관 — “같은 잠재 효과, 그룹별 다른 폭”

\(\theta_i\) 는 표준 정규 한 개. MZ 쌍이면 \(\sigma_{v(MZ)}\) 만큼 증폭, DZ 쌍이면 \(\sigma_{v(DZ)}\) 만큼 증폭. 결과적으로:

MZ 쌍 안에서 두 개체의 잠재 변수 공통 부분 = \(\sigma_{v(MZ)} \theta_i\) → 분산 \(\sigma_{v(MZ)}^2\).
DZ 쌍 안에서 두 개체의 잠재 변수 공통 부분 = \(\sigma_{v(DZ)} \theta_i\) → 분산 \(\sigma_{v(DZ)}^2\).

같은 모형 구조, 하나의 \(T\) 벡터, 두 가지 그룹별 분산 — design vector 가 mode 를 결정.

유전학적 해석: MZ 쌍은 유전자 100 % 공유, DZ 쌍은 평균 50 % 공유. 잠재 형질의 유전적 기여가 크면 MZ 쌍 안 상관이 DZ 쌍 안 상관보다 크게 나타난다 (즉 \(\sigma_{v(MZ)} > \sigma_{v(DZ)}\)). 두 분산을 비교하면 유전성 (heritability) 의 단서를 얻는다.

5.3 Tetrachoric Correlation — Probit 형태에서

Probit + intercept-only → tetrachoric correlation

식 (9.25) 를 probit link 로 쓰고 (logit 대신), 공변량을 절편 하나로 단순화하면 그룹별 ICC 가:

\[ ICC_{MZ} = \frac{\sigma_{v(MZ)}^2}{\sigma_{v(MZ)}^2 + 1}, \quad ICC_{DZ} = \frac{\sigma_{v(DZ)}^2}{\sigma_{v(DZ)}^2 + 1} \]

이 값들이 쌍둥이 쌍 안의 tetrachoric correlation — 두 이항 변수가 잠재 정규 변수에서 절단된 것이라 가정할 때의 잠재 정규 변수 사이 상관계수.

직관 — 왜 probit 형태에서 자연스러운가

Tetrachoric correlation 은 정의상 잠재 정규 변수의 상관. probit 모형이 정규 잠재 변수를 가정하므로 (logistic 잠재 변수 대신), tetrachoric 해석이 직접 맞는다. logistic 형태에서는 잠재 분포가 logistic 이라 polychoric/tetrachoric 의 정확한 정의에서 벗어난다.

유전학·가족 연구에서 probit 을 선호하는 이유: 형질 유전 모형이 정규 잠재 변수를 가정하는 전통이 있고, tetrachoric 으로 추정된 상관이 곧바로 유전성 모수 (heritability) 의 입력이 된다.

공변량을 추가하면 “조정된 (adjusted) tetrachoric correlation” 이 된다 — 다른 변수를 통제한 상태에서의 잠재 변수 상관.

5.4 응용 2 — IRT 2PL 모형 (식 9.27-9.28)

IRT 데이터 구조

\(i = 1, \ldots, N\): 응시자 (level-2 unit).
\(j = 1, \ldots, m\): 시험 문항 (level-1 unit). 각 응시자가 \(m\) 개 문항에 응답.
\(Y_{ij} \in \{0, 1\}\): 정답 여부.
\(\theta_i\): 응시자의 잠재 능력 (latent ability), \(\mathcal{N}(0, 1)\).

문항 모수:

\(a_j\): discrimination parameter (변별도) — 능력 차이를 얼마나 잘 구분하는가.
\(b_j\): difficulty parameter (난이도) — 능력 평균 응시자가 정답 확률 0.5 가 되는 능력 수준.

식 (9.26-9.28) — 2PL 의 세 가지 등가 표현

원형 (식 9.26):

\[ P(Y_{ij} = 1 \mid \theta_i) = \frac{1}{1 + \exp[-a_j(\theta_i - b_j)]} \tag{9.26} \]

\(c_j = -a_j b_j\) 로 재표기 (식 9.27):

\[ P(Y_{ij} = 1 \mid \theta_i) = \frac{1}{1 + \exp[-(c_j + a_j \theta_i)]} \tag{9.27} \]

logit 형태로 (식 9.28):

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = c_j + a_j \theta_i \tag{9.28} \]

직관 — IRT 가 mixed-effects logistic 의 특수 경우

식 (9.28) 을 식 (9.24) 와 비교:

\(x_{ij}^\top \beta\) 자리에 문항별 절편 \(c_j\) — design 으로 표현하면 \(X\) 가 문항 dummy 행렬, \(\beta = (c_1, \ldots, c_m)^\top\).
\(z_{ij}^\top T \theta_i\) 자리에 \(a_j \theta_i\) — design 으로 표현하면 \(z_{ij}\) 가 \(m \times 1\) 문항 dummy 벡터 (현재 문항 위치만 1), \(T = (a_1, a_2, \ldots, a_m)^\top\) (\(m \times 1\) 벡터).

→ 2PL = “각 문항에 별도의 randomeffect 분산을 허용한 mixed-effects logistic”. 문항별 dummy 가 design vector 에 들어가 § 9.5.3 의 이질 분산 framework 와 정확히 같은 형태.

응시자 능력 \(\theta_i\) 는 표준 정규 한 개 (단일 차원), 문항별 변별도 \(a_j\) 가 그 능력의 “확대 배율”. 변별도가 큰 문항은 능력 차이를 크게 증폭, 작은 문항은 미미하게 증폭.

문항 행렬 표현:

\[ \begin{bmatrix} l_{1i} \\ l_{2i} \\ l_{3i} \\ l_{4i} \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ c_4 \end{bmatrix} + \begin{bmatrix} a_1 & 0 & 0 & 0 \\ 0 & a_2 & 0 & 0 \\ 0 & 0 & a_3 & 0 \\ 0 & 0 & 0 & a_4 \end{bmatrix} \theta_i \cdot \mathbf{1} \]

(여기서 각 행이 해당 문항의 logit. 두 번째 행렬의 비대각 0 이 문항 dummy 의 작동.)

Rasch 모형 — 1PL 특수 경우

2PL 에서 \(a_1 = a_2 = \cdots = a_m = a\) 로 모든 문항의 변별도를 같게 제약하면:

\[ \log \left[ \frac{p_{ij}}{1 - p_{ij}} \right] = c_j + a \theta_i \]

이는 랜덤 절편 logistic 회귀모형 + 문항 dummy 공변량 — 가장 단순한 GLMM 이항. Rasch (1960) 가 제안한 1PL (one-parameter logistic) 모형.

→ Rasch (1PL) ⊂ 2PL ⊂ 일반 mixed-effects logistic. 같은 모형 가족, 제약 정도만 다름.

통합 시각의 가치 — 왜 이 framework 가 중요한가

쌍둥이 분석과 IRT 는 외형상 매우 다른 응용 — 한쪽은 유전학, 한쪽은 교육 측정. 그러나 식 (9.24) 의 mixed-effects logistic 으로 보면 둘 다 “design vector 가 dummy code 인 GLMM 의 특수 경우”:

쌍둥이: design vector = 쌍 type (MZ/DZ) dummy.
IRT 2PL: design vector = 문항 dummy.

이 통합이 실무적으로 가치 있는 이유:

공통 알고리즘: 같은 GLMM 추정 코드 (Gauss-Hermite 적분 등) 로 두 응용 모두 처리.
공변량 자유 추가: 전통적 IRT 는 응시자/문항 모수만 다루지만, mixed-effects 표현에서는 응시자 특성 (성별, 학년) 이나 문항 특성 (주제, 형식) 을 공변량으로 자연스럽게 추가 가능. 차별 기능 분석 (DIF, differential item functioning) — 같은 능력의 다른 그룹이 같은 문항을 다르게 푸는지 검토 — 가 회귀 계수 검정으로 단순화.
모형 비교: 1PL vs 2PL, MZ-DZ 동등 분산 vs 이질 분산 — 모두 일반 LR 검정으로 처리.

§ 9.5.3 의 짧은 한 페이지가 사실상 “mixed-effects logistic 의 표현력이 어디까지 닿는가” 의 시연이다. design vector 의 자유로운 선택이 모형 가족을 폭발적으로 확장한다는 것이 핵심 메시지.

6 응용 분야

분야	활용	구체적 예시
임상 종단	처치별 치유 확률 추적	HAMD ≤ 7 관해 여부의 환자별 종단 변화, 약물 vs 위약 비교
교육 평가	응시자 능력·문항 난이도 동시 추정	표준화 시험의 IRT 2PL/Rasch 적합
행동 유전학	형질의 유전적 기여 분리	쌍둥이 데이터의 MZ/DZ tetrachoric correlation 비교로 heritability 추정
보건 서비스	클리닉 간 진료 결과 차이	다기관 임상시험에서 클리닉 ICC 측정
A/B 테스트	사용자별 전환 성향 모형화	같은 사용자의 반복 노출에서 사용자 간 ICC 측정, GEE 와 비교

7 코드 예시

7.1 Step 1: ICC 직접 계산 (numpy)

import numpy as np


def icc_logistic(sigma_v: float) -> float:
    """Logistic mixed-effects 의 ICC.
    분모는 sigma_v^2 + pi^2/3.
    """
    return sigma_v ** 2 / (sigma_v ** 2 + np.pi ** 2 / 3)


def icc_probit(sigma_v: float) -> float:
    """Probit mixed-effects 의 ICC. 분모는 sigma_v^2 + 1."""
    return sigma_v ** 2 / (sigma_v ** 2 + 1)


# 다양한 sigma_v 에서 두 ICC 비교
for sigma_v in [0.5, 1.0, 1.4, 2.0, 3.0]:
    rho_logit = icc_logistic(sigma_v)
    rho_probit = icc_probit(sigma_v)
    print(f"sigma_v = {sigma_v:.1f}: "
          f"ICC_logit = {rho_logit:.3f}, "
          f"ICC_probit = {rho_probit:.3f}")

결과 해석 포인트

같은 \(\sigma_v\) 에서 probit 의 ICC 가 logit 보다 항상 큰 이유 — 분모의 잔차 분산이 logit (\(\pi^2/3 \approx 3.29\)) 보다 probit (\(1\)) 이 작기 때문. 즉 probit 척도에서 같은 분산이 ICC 비중으로는 더 크게 보인다.

이 차이는 scale 의 차이일 뿐 실질적 차이가 아님 — 두 모형이 같은 데이터에 적합되면 ICC 가 다르게 나오지만, 잠재 변수 분산을 같은 척도로 환산하면 동일한 결론.

7.2 Step 2: Cholesky 분해의 작동 확인

import numpy as np


def build_cov_from_cholesky(t11: float, t21: float, t22: float) -> np.ndarray:
    """2x2 Cholesky factor 로부터 분산-공분산 행렬 복원.

    T = [[t11, 0], [t21, t22]]
    Sigma = T @ T.T
    """
    T = np.array([[t11, 0.0], [t21, t22]])
    return T @ T.T


def cov_to_correlation(Sigma: np.ndarray) -> float:
    """2x2 분산-공분산에서 상관계수 추출."""
    return Sigma[0, 1] / np.sqrt(Sigma[0, 0] * Sigma[1, 1])


# 예: 절편 SD 1.0, 기울기 SD 0.6, 상관 -0.3 인 분산 구조를 만들고 싶다면
# t11 = 1.0
# t21 = sigma_intercept * sigma_slope * rho / sigma_intercept = sigma_slope * rho
# t22 = sigma_slope * sqrt(1 - rho^2)
sigma_int, sigma_slope, rho = 1.0, 0.6, -0.3
t11 = sigma_int
t21 = sigma_slope * rho
t22 = sigma_slope * np.sqrt(1 - rho ** 2)

Sigma = build_cov_from_cholesky(t11, t21, t22)
print("Reconstructed Sigma_v:")
print(Sigma)
print(f"  Var(intercept) = {Sigma[0,0]:.3f} (expected {sigma_int**2:.3f})")
print(f"  Var(slope)     = {Sigma[1,1]:.3f} (expected {sigma_slope**2:.3f})")
print(f"  Cor(int,slope) = {cov_to_correlation(Sigma):+.3f} (expected {rho:+.3f})")

# 양정치 확인 — 고유값이 모두 양수
eigvals = np.linalg.eigvalsh(Sigma)
print(f"\nEigenvalues: {eigvals}  → all positive? {np.all(eigvals > 0)}")

검증 포인트

\(T\) 의 자유 추정 \((t_{11}, t_{21}, t_{22})\) 가 어떤 실수값이든 — 심지어 \(t_{21}\) 이 음수여도 — \(\Sigma_v = TT^\top\) 는 항상 양정치. 이것이 직접 추정 (분산·공분산을 따로) 대비 가장 큰 장점.

만약 \(t_{22} \to 0\) 이면 기울기 분산이 0 에 접근 — 모형이 “기울기 변동 없음” 으로 수렴. boundary 에서도 \(\Sigma_v\) 가 정의되어 ML 표면이 안정.

7.3 Step 3: 쌍둥이 이질 분산 시뮬레이션 (numpy)

import numpy as np


def simulate_twin_data(n_pairs: int, sigma_mz: float, sigma_dz: float,
                       beta0: float, p_mz: float = 0.5,
                       seed: int = 2026) -> dict:
    """이질 분산 mixed-effects logistic 으로 쌍둥이 이항 데이터 생성.

    각 쌍에 단일 잠재 효과 theta_i ~ N(0,1).
    MZ 쌍은 sigma_mz, DZ 쌍은 sigma_dz 로 증폭.
    쌍 안의 두 개체는 같은 잠재 효과를 공유 (상관 발생).
    """
    rng = np.random.default_rng(seed)

    # 쌍 type 결정
    is_mz = rng.binomial(1, p_mz, size=n_pairs).astype(bool)

    # 각 쌍의 잠재 효과
    theta = rng.normal(0, 1, size=n_pairs)
    sigma_pair = np.where(is_mz, sigma_mz, sigma_dz)
    upsilon = sigma_pair * theta

    # 쌍 안의 두 개체 (j=1,2) — 같은 upsilon 공유
    Y = np.zeros((n_pairs, 2), dtype=int)
    for j in range(2):
        eta = beta0 + upsilon  # 절편만 있는 모형
        p = 1.0 / (1.0 + np.exp(-eta))
        Y[:, j] = rng.binomial(1, p)

    return {"Y": Y, "is_mz": is_mz, "upsilon": upsilon}


# 시뮬레이션 — MZ 쌍 분산이 DZ 보다 큼 (유전 기여 큼)
data = simulate_twin_data(n_pairs=2000, sigma_mz=1.5, sigma_dz=0.8,
                          beta0=0.0, p_mz=0.5)

# 쌍 안 일치율 (concordance) 계산
concordance_mz = np.mean(data["Y"][data["is_mz"], 0] == data["Y"][data["is_mz"], 1])
concordance_dz = np.mean(data["Y"][~data["is_mz"], 0] == data["Y"][~data["is_mz"], 1])

print(f"MZ 쌍 일치율: {concordance_mz:.3f}")
print(f"DZ 쌍 일치율: {concordance_dz:.3f}")
print(f"\n예상 ICC (probit 척도):")
print(f"  ICC_MZ = {1.5**2 / (1.5**2 + 1):.3f}")
print(f"  ICC_DZ = {0.8**2 / (0.8**2 + 1):.3f}")

결과 해석

MZ 쌍 분산이 더 크면 (유전 기여가 더 크면) 쌍 안의 두 개체가 같은 응답을 할 확률 (일치율) 이 더 높다. 이 일치율의 차이가 heritability 의 단서.

실무에서는 일치율 차이를 직접 분석하지 않고, mixed-effects logistic 으로 \(\sigma_{v(MZ)}, \sigma_{v(DZ)}\) 를 추정해 ICC 비교 또는 LR 검정 (\(H_0: \sigma_{v(MZ)} = \sigma_{v(DZ)}\)) 으로 통계적 추론.

7.4 Step 4: statsmodels 실무 적합 — 종단 데이터 GLMM

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


# 종단 이항 데이터 시뮬레이션 (시간 × 처치 효과 + 환자 랜덤 절편)
def simulate_longitudinal_binary(n_subjects: int, n_times: int,
                                 beta0: float, beta_treat: float,
                                 beta_time: float, sigma_v: float,
                                 seed: int = 2026) -> pd.DataFrame:
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n_subjects):
        treat = rng.binomial(1, 0.5)
        upsilon = rng.normal(0, sigma_v)  # 환자별 랜덤 절편
        for j in range(n_times):
            t = j  # 0, 1, ..., n_times-1
            eta = beta0 + beta_treat * treat + beta_time * t + upsilon
            p = 1.0 / (1.0 + np.exp(-eta))
            y = rng.binomial(1, p)
            rows.append({"subject": i, "time": t, "treat": treat, "y": y})
    return pd.DataFrame(rows)


df = simulate_longitudinal_binary(n_subjects=300, n_times=5,
                                  beta0=-1.0, beta_treat=0.8, beta_time=0.3,
                                  sigma_v=1.2)

# Mixed-effects logistic — BinomialBayesMixedGLM 또는 GEE 대안
# statsmodels 의 BinomialBayesMixedGLM 으로 적합
formula = "y ~ treat + time"
random = {"a": "0 + C(subject)"}  # subject별 랜덤 절편

model = sm.BinomialBayesMixedGLM.from_formula(formula, random, df)
result = model.fit_vb()  # variational Bayes

print(result.summary())
print(f"\n추정된 ICC = sigma_v^2 / (sigma_v^2 + pi^2/3) = "
      f"{result.vc_mean[0]**2 / (result.vc_mean[0]**2 + np.pi**2/3):.3f}")

실무 노트 — 추정 방법 선택

statsmodels 의 BinomialBayesMixedGLM 은 변분 근사 (variational Bayes) 기반. 정확한 Gauss-Hermite 구적법 (§ 9.6) 은 R 의 lme4::glmer 또는 MASS::glmmPQL 가 더 표준적.

R 코드 등가:

library(lme4)
fit <- glmer(y ~ treat + time + (1 | subject),
             data = df, family = binomial)
summary(fit)

# ICC 계산
sigma2_v <- as.numeric(VarCorr(fit)$subject)
icc <- sigma2_v / (sigma2_v + pi^2 / 3)

추정 방법 비교:

방법	정확도	속도	적합 패키지
Gauss-Hermite quadrature	높음	느림	R `lme4::glmer(nAGQ > 1)`
Adaptive quadrature	매우 높음	보통	SAS NLMIXED, Stata `meqrlogit`
Laplace approximation	보통	빠름	R `lme4::glmer` (default)
MCMC (full Bayes)	매우 높음	매우 느림	R `brms`, `rstanarm`
Variational Bayes	보통	빠름	Python `statsmodels`

랜덤 효과가 1-2 차원이면 quadrature 가 표준. 차원이 늘어나면 (3 이상) Laplace 또는 MCMC.

8 관련 주제

선행 지식

Ch.9 Overview — 이항 GLMM 의 전체 framework
§ 9.2-9.3 — Single-level 로지스틱·프로빗 회귀
§ 9.4 — Threshold concept (잠재 변수 framework)
§ 4.2-4.3 — 정규 종단 MRM 의 랜덤 절편 모형 (구조적 등가물)

후속 주제

§ 9.5.4 — Multilevel Representation (level-1 / level-2 분해 표기)
§ 9.5.5 — Response Functions (logit, probit, complementary log-log)
§ 9.6 — Marginal MLE 추정과 Gauss-Hermite quadrature
§ 9.11 — Subject-specific vs population-averaged 효과의 정량 비교
Ch.10 GLMM 순서형 — 누적 logit 으로의 자연 확장
Ch.11 GLMM 명목 — 기준 범주 logit 으로의 확장

관련 개념

GEE (Ch.8) — Marginal vs subject-specific 비교의 다른 한쪽
mm-06 GLMM 이진 — AI Agent 비즈니스 예시 직관
IRT 2PL/Rasch — 교육 측정에서의 mixed-effects logistic 응용
Tetrachoric correlation — 행동 유전학·가족 연구