Kwangmin Kim - 혼합효과 회귀모형 — 연속 반응변수 개요

1 왜 혼합효과 회귀모형인가

ANOVA(구형성 가정)와 MANOVA(완전 데이터 요구)는 종단 데이터 분석에서 두 가지 근본적 한계를 가진다.

한계	ANOVA	MANOVA	MRM
결측 피험자	목록삭제 또는 제한적 처리	완전 제외	모든 피험자 포함
분산-공분산 구조	구형성(compound symmetry) 강요	자유 구조, 단 모수 폭발	유연한 중간 구조
불균등 시점	불가	불가	피험자마다 다른 시점 허용
시변 공변량	불가	불가	포함 가능
개인별 추세	불가 (집단 평균만)	불가	경험 베이즈로 추정

혼합효과 회귀모형(Mixed-effects Regression Model, MRM)은 이 모든 한계를 동시에 해결한다. 랜덤 피험자 효과를 회귀모형에 포함시켜 관측값의 상관 구조를 설명하고, 피험자별 개인 추세를 추정한다.

MRM은 다양한 이름으로 불린다: 랜덤효과 모형(Laird & Ware, 1982), 계층 선형 모형(Raudenbush & Bryk, 2002), 다수준 모형(Goldstein, 1995), 혼합 모형(Wolfinger, 1993). 핵심은 동일하다 — 개인 수준의 랜덤 효과를 포함한 회귀 모형 (Hedeker & Gibbons, 2006, Ch.4).

2 출발점: 단순 선형회귀

MRM을 이해하는 가장 쉬운 방법은 익숙한 단순 선형회귀에서 시작하여 단계별로 확장하는 것이다.

2.1 독립 관측 가정의 한계 (식 4.1)

\[ y_{ij} = \beta_0 + \beta_1 t_{ij} + \varepsilon_{ij} \]

\(i = 1, \ldots, N\): 피험자
\(j = 1, \ldots, n_i\): 피험자 \(i\) 의 관측 횟수
\(\varepsilon_{ij} \overset{iid}{\sim} N(0, \sigma^2)\): 독립 동분산 오차

문제: 종단 데이터에서 같은 피험자의 반복 관측은 독립이 아니다. 또한 이 모형은 모든 피험자가 동일한 기울기 \(\beta_1\) 을 갖는다고 가정한다 — 개인 간 변화 이질성을 완전히 무시한다.

독립 가정 위반의 결과

OLS로 추정하면: - 표준오차가 과소 추정 → 1종 오류 증가 - 같은 피험자의 관측값이 주는 “새로운 정보”를 과대평가

MRM은 개인 랜덤 효과를 추가하여 이 상관을 명시적으로 모형화한다.

3 랜덤 절편 모형 (Random Intercept MRM)

3.1 모형 구조 (식 4.2)

단순 선형회귀에 피험자별 랜덤 절편 \(v_{0i}\) 를 추가한다:

\[ y_{ij} = \beta_0 + \beta_1 t_{ij} + v_{0i} + \varepsilon_{ij} \]

\(v_{0i} \sim N(0, \sigma^2_v)\): 피험자 \(i\) 의 기저 수준 편차 (랜덤 효과)
\(\varepsilon_{ij} \sim N(0, \sigma^2)\): 조건부 독립 잔차 (랜덤 효과 조건부)

3.2 위계적 표현 (수준 1, 수준 2)

같은 모형을 위계적으로 분해하면 구조가 더 명확해진다 (식 4.3~4.4):

수준 1 (피험자 내, Within-subject):

\[ y_{ij} = b_{0i} + b_{1i} t_{ij} + \varepsilon_{ij} \]

수준 2 (피험자 간, Between-subject):

\[ b_{0i} = \beta_0 + v_{0i}, \quad b_{1i} = \beta_1 \]

직관: 수준-1은 각 피험자 내에서 시간에 따른 변화를 기술한다. 수준-2는 피험자별 절편(\(b_{0i}\))이 집단 평균 \(\beta_0\) 주위에 분포하되, 기울기(\(b_{1i}\))는 모든 사람이 \(\beta_1\) 로 동일하다고 가정한다. 두 수준을 합치면 단일 방정식 (식 4.2)가 된다.

시각적 이해: 모든 피험자의 추세선이 서로 평행하다. 출발점(절편)은 다르지만, 기울기는 동일하다. 차이는 상수 \(v_{0i}\) 만큼 일정하게 유지된다.

3.3 복합대칭과 ICC

랜덤 절편 모형이 함의하는 공분산 구조 (식 4.5):

\[ V(y_{ij}) = \sigma^2_v + \sigma^2 \quad \text{(분산, 시점 불변)} \]

\[ C(y_{ij}, y_{ij'}) = \sigma^2_v, \quad j \neq j' \quad \text{(공분산, 모든 시차 동일)} \]

이는 복합대칭(compound symmetry) 구조다. 급내상관계수(ICC)는:

\[ \text{ICC} = \frac{\sigma^2_v}{\sigma^2_v + \sigma^2} \]

해석: ICC는 총 변동 중 피험자 수준에 귀속되는 비율이다. ICC가 클수록 같은 피험자의 반복 측정이 높은 상관을 가지며, 독립 가정을 무시하는 분석의 왜곡이 커진다.

4 랜덤 절편 + 기울기 모형 (Random Intercept and Trend MRM)

4.1 왜 기울기도 랜덤이어야 하는가

현실의 종단 데이터에서 모든 피험자가 동일한 속도로 변화한다고 가정하는 것은 과도하게 단순하다.

일부 환자는 치료에 빠르게 반응하고, 다른 환자는 반응이 없다.
분산이 시간에 따라 증가하는 경우(치료 반응의 이질성), 랜덤 절편만으로는 이를 포착할 수 없다.

4.2 모형 구조 (식 4.6)

수준-2 모형을 확장하여 기울기도 개인별로 다르게 한다:

\[ b_{0i} = \beta_0 + v_{0i}, \quad b_{1i} = \beta_1 + v_{1i} \]

합치면:

\[ y_{ij} = \underbrace{(\beta_0 + v_{0i})}_{\text{개인 절편}} + \underbrace{(\beta_1 + v_{1i})}_{\text{개인 기울기}} t_{ij} + \varepsilon_{ij} \]

랜덤 효과의 결합 분포:

\[ \begin{bmatrix} v_{0i} \\ v_{1i} \end{bmatrix} \sim \mathcal{N}\left(\mathbf{0},\; \Sigma_v = \begin{bmatrix} \sigma^2_{v_0} & \sigma_{v_0 v_1} \\ \sigma_{v_0 v_1} & \sigma^2_{v_1} \end{bmatrix}\right) \]

모수	의미
\(\sigma^2_{v_0}\)	절편 이질성 — 개인별 기저 수준이 얼마나 퍼져 있는가
\(\sigma^2_{v_1}\)	기울기 이질성 — 개인별 변화 속도가 얼마나 다른가
\(\sigma_{v_0 v_1}\)	절편-기울기 공분산 — 초기값이 높은 사람이 더 빠르게/느리게 변하는가

절편-기울기 공분산의 직관

\(\sigma_{v_0 v_1} < 0\): 초기값이 높은 사람이 더 빠르게 감소 (천장 효과 또는 치료 반응)

\(\sigma_{v_0 v_1} > 0\): 초기값이 높은 사람이 더 빠르게 증가 (성장 가속)

\(\sigma_{v_0 v_1} = 0\): 초기값과 변화 속도가 무관

4.3 시간 코딩의 중요성 (§4.4.2)

절편 관련 모수(\(\beta_0\), \(v_{0i}\), \(\sigma^2_{v_0}\))의 의미는 시간 변수 \(t\) 의 영점(origin) 에 의존한다.

코딩 방식	예시 (\(n=6\))	절편 의미	주의사항
기저선 코딩	\(t = 0,1,2,3,4,5\)	기저선(\(t=0\)) 수준	가장 직관적
중심화 코딩	\(t = -2.5, -1.5, \ldots, 2.5\)	연구 중간점 수준	\(\sigma^2_{v_0}\) 값 달라짐
종료점 코딩	\(t = -5, -4, \ldots, 0\)	최종 시점 수준	역방향 해석 필요

핵심: 기울기 추정값 \(\hat{\beta}_1\) 과 기울기 분산 \(\hat{\sigma}^2_{v_1}\) 은 코딩에 무관하게 동일하다. 그러나 절편 추정값과 절편 분산, 그리고 절편-기울기 공분산의 부호까지 달라질 수 있다.

5 공변량 포함: 집단 효과와 교호작용

5.1 시불변 공변량의 수준-2 투입 (식 4.7)

진단명(DX: 내인성=1 vs 비내인성=0)처럼 피험자 수준 공변량은 수준-2 모형에 투입한다:

\[ b_{0i} = \beta_0 + \beta_2 DX_i + v_{0i} \] \[ b_{1i} = \beta_1 + \beta_3 DX_i + v_{1i} \]

\(\beta_0\): 비내인성 집단의 기저선 평균
\(\beta_1\): 비내인성 집단의 주당 변화율
\(\beta_2\): 내인성 집단의 기저선 차이
\(\beta_3\): 집단 × 시간 교호작용 — 두 집단의 변화율 차이

합치면 단일 방정식:

\[ y_{ij} = \beta_0 + \beta_1 t_{ij} + \beta_2 DX_i + \beta_3 (DX_i \cdot t_{ij}) + v_{0i} + v_{1i} t_{ij} + \varepsilon_{ij} \]

직관: 행렬 표현에서 이것은 그냥 \(X_i\) 에 열을 추가하는 것뿐이다. 수준-1/수준-2 분리는 해석의 편의를 위한 것이며 실제로는 하나의 모형이 추정된다.

6 행렬 표현 (Matrix Formulation)

6.1 일반 MRM (식 4.8)

\[ \underbrace{y_i}_{n_i \times 1} = \underbrace{X_i}_{n_i \times p} \underbrace{\beta}_{p \times 1} + \underbrace{Z_i}_{n_i \times r} \underbrace{v_i}_{r \times 1} + \underbrace{\varepsilon_i}_{n_i \times 1} \]

기호	의미	차원
\(y_i\)	피험자 \(i\) 의 반응 벡터	\(n_i \times 1\)
\(X_i\)	고정효과 설계 행렬	\(n_i \times p\)
\(\beta\)	고정효과 계수 벡터	\(p \times 1\)
\(Z_i\)	랜덤효과 설계 행렬	\(n_i \times r\)
\(v_i\)	개인 랜덤 효과 벡터	\(r \times 1\)
\(\varepsilon_i\)	잔차 벡터	\(n_i \times 1\)

랜덤 절편 + 기울기 모형에서 (\(r=2\)):

\[ X_i = Z_i = \begin{bmatrix} 1 & t_{i1} \\ 1 & t_{i2} \\ \vdots & \vdots \\ 1 & t_{in_i} \end{bmatrix}, \quad \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}, \quad v_i = \begin{bmatrix} v_{0i} \\ v_{1i} \end{bmatrix} \]

6.2 분포 가정

\[ \varepsilon_i \sim \mathcal{N}(0, \sigma^2 I_{n_i}), \quad v_i \sim \mathcal{N}(0, \Sigma_v) \]

이로부터 \(y_i\) 의 주변 분포:

\[ y_i \sim \mathcal{N}(X_i \beta, \; \underbrace{Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i}}_{V_i}) \]

\(V_i = Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i}\) 가 핵심이다. 이것이 반복 측정의 분산-공분산 행렬을 결정한다. 랜덤 효과 구조 \(\Sigma_v\) 와 설계 행렬 \(Z_i\) 의 조합으로 복합대칭부터 비구조적까지 다양한 공분산 패턴을 표현할 수 있다 (식 4.12).

공분산 구조의 유연성

랜덤 절편만 (\(Z_i = \mathbf{1}\)):

\[V_i = \sigma^2_v \mathbf{1}\mathbf{1}' + \sigma^2 I \quad \text{(복합대칭)}\]

랜덤 절편 + 기울기 (\(n_i=3\), \(t=0,1,2\)):

\[V_i = \sigma^2 I + Z_i \Sigma_v Z_i' = \begin{bmatrix} \sigma^2_{v_0}+\sigma^2 & \sigma^2_{v_0}+\sigma_{v_0v_1} & \sigma^2_{v_0}+2\sigma_{v_0v_1} \\ \cdots & \sigma^2_{v_0}+2\sigma_{v_0v_1}+\sigma^2_{v_1}+\sigma^2 & \cdots \\ \cdots & \cdots & \sigma^2_{v_0}+4\sigma_{v_0v_1}+4\sigma^2_{v_1}+\sigma^2 \end{bmatrix}\]

분산과 공분산이 시점에 따라 달라진다 — 복합대칭을 넘는 유연성이다.

7 경험 베이즈 추정 (Empirical Bayes)

7.1 개인별 랜덤 효과 추정 (식 4.10)

모형 모수 \(\hat{\beta}\), \(\hat{\Sigma}_v\), \(\hat{\sigma}^2\) 을 추정한 후, 개인별 랜덤 효과 \(v_i\) 를 추정한다:

\[ \hat{v}_i = \left[Z_i'(\sigma^2 I)^{-1}Z_i + \Sigma_v^{-1}\right]^{-1} Z_i'(\sigma^2 I)^{-1}(y_i - X_i\hat{\beta}) \]

직관: 이 수식은 두 정보를 결합한다.

데이터 기반 추정 \(Z_i'(\sigma^2 I)^{-1}(y_i - X_i\hat{\beta})\): 피험자 \(i\) 의 잔차로부터 개인 효과를 직접 추정
모집단 사전분포 \(\Sigma_v^{-1}\): 랜덤 효과가 모집단 분포 \(N(0, \Sigma_v)\) 를 따른다는 정보

관측 횟수 \(n_i\) 가 많을수록 데이터 기반 추정에 더 많은 가중치가 부여된다. \(n_i\) 가 작으면 모집단 평균 쪽으로 수축(shrinkage)된다 — 이것이 “경험 베이즈” 추정의 핵심이다.

개인별 성장 곡선 추정값:

\[ \hat{b}_{0i} = \hat{\beta}_0 + \hat{v}_{0i}, \quad \hat{b}_{1i} = \hat{\beta}_1 + \hat{v}_{1i} \]

이를 통해 집단 평균 추세 외에 개인별 추세를 추정할 수 있다 — ANOVA/MANOVA로는 불가능한 기능이다.

8 추정: ML vs REML

8.1 최대가능도(ML) 추정

MRM 모수 \((\beta, \Sigma_v, \sigma^2)\) 는 반응 벡터의 주변 가능도를 최대화하여 추정한다:

\[ L(\beta, \Sigma_v, \sigma^2) = \prod_{i=1}^{N} f(y_i \mid \beta, \Sigma_v, \sigma^2) \]

여기서 \(y_i \sim \mathcal{N}(X_i\beta, Z_i\Sigma_v Z_i' + \sigma^2 I)\) 이므로 이 가능도는 다변량 정규 밀도의 곱이다.

8.2 REML (제한 최대가능도)

ML은 고정효과를 추정할 때 발생하는 자유도 손실을 고려하지 않아 분산 성분을 과소 추정하는 경향이 있다 (소표본에서 더 뚜렷). REML은 고정효과를 소거한 잔차의 가능도를 최대화하여 분산 성분의 편향을 교정한다.

비교	ML	REML
목적	모든 모수 동시 추정	분산 성분 편향 교정
분산 추정	과소 추정 (소표본)	불편 추정
모형 비교	고정효과 차이 포함한 nested model 비교 가능	같은 고정효과 구조 내에서만 비교 가능
기본값	—	대부분의 소프트웨어 기본

8.3 가설 검정

고정효과 \(\beta\): Wald 검정 (\(Z = \hat{\beta}/\text{SE}\)), 표준 정규 분포 참조

분산 성분 \(\sigma^2_v\): Wald 검정 부적합 — 분산은 0 이상이어야 하므로 경계값 문제 발생. 대신 우도비 검정(LRT)를 사용하며, 경계 문제로 인해 얻은 \(p\)-값을 2로 나눠 보수적 보정을 적용한다.

중첩 모형 비교: LRT 통계량 \(= -2\log L_{\text{제한}} - (-2\log L_{\text{확장}})\), \(\chi^2\) 분포(자유도 = 추가 모수 수)로 평가

9 Reisby 정신과 데이터 예시

9.1 연구 설계

Reisby et al.(1977)의 연구는 우울증 입원 환자 66명에게 이미프라민(imipramine) 투여 후 해밀턴 우울 척도(HDRS)를 6주간 추적 측정했다. HDRS 점수가 낮을수록 우울 증상이 적다.

피험자: N=66 (내인성 37명, 비내인성 29명)
시점: Week 0(기저선), 1, 2, 3, 4, 5
완전 데이터 보유 피험자: 46명 → MANOVA는 20명 제외

9.2 시점별 요약 통계 (Table 4.1)

	주 0	주 1	주 2	주 3	주 4	주 5
평균	23.44	21.84	18.31	16.42	13.62	11.95
표준편차	4.53	4.70	5.49	6.42	6.97	7.22
\(n\)	61	63	65	65	63	58

주목: 평균은 감소(치료 효과), 표준편차는 증가(개인 반응 이질성). 이 패턴은 복합대칭을 명백히 위반하며, 랜덤 기울기 모형이 필요함을 시사한다.

9.3 랜덤 절편 모형 결과 (Table 4.3)

\[ y_{ij} = \beta_0 + \beta_1 \cdot \text{Week}_{ij} + v_{0i} + \varepsilon_{ij} \]

모수	추정값	SE	Z	\(p\)
\(\hat{\beta}_0\) (절편)	23.55	0.64	36.80	<.0001
\(\hat{\beta}_1\) (기울기)	−2.38	0.14	−17.00	<.0001
\(\hat{\sigma}^2_{v_0}\)	16.15	3.41	—	—
\(\hat{\sigma}^2\)	19.04	1.53	—	—

\(-2\log L = 2285.19\)

해석: - 환자들은 평균적으로 기저선에서 HDRS=23.55이고, 주당 2.38점씩 감소한다. - ICC = 16.15/(16.15+19.04) = 0.46 — 설명되지 않은 분산의 46%가 개인 수준 - 주별 표준편차가 증가하는 패턴이 이 모형에서는 포착되지 않는다 (복합대칭 한계)

9.4 랜덤 절편 + 기울기 모형 결과 (Table 4.5)

모수	추정값	SE
\(\hat{\beta}_0\)	23.58	0.55
\(\hat{\beta}_1\)	−2.38	0.21
\(\hat{\sigma}^2_{v_0}\) (절편 분산)	12.63	3.53
\(\hat{\sigma}_{v_0 v_1}\) (절편-기울기 공분산)	−1.42	1.04
\(\hat{\sigma}^2_{v_1}\) (기울기 분산)	2.08	0.52
\(\hat{\sigma}^2\) (잔차 분산)	12.22	1.12

\(-2\log L = 2219.04\)

우도비 검정 (랜덤 절편 vs 랜덤 절편+기울기):

\[ \Delta(-2\log L) = 2285.19 - 2219.04 = 66.15, \quad df=2, \quad p \ll .001 \]

복합대칭 가정이 강하게 기각된다. 랜덤 기울기 모형이 유의하게 더 적합하다.

기울기 이질성 해석: - 집단 평균 기울기: \(\hat{\beta}_1 = -2.38\) - 기울기 표준편차: \(\sqrt{2.08} = 1.44\) - 95% 범위: \(-2.38 \pm 1.96 \times 1.44 = [-5.20, \; +0.44]\)

양수 기울기가 범위에 포함된다 — 일부 환자는 시간이 지나도 호전되지 않거나 악화될 수 있다.

절편-기울기 공분산 해석: - \(\hat{\sigma}_{v_0 v_1} = -1.42\): 상관으로 변환하면 \(r = -1.42/\sqrt{12.63 \times 2.08} = -0.28\) - 해석: 기저선 HDRS가 높은 환자(더 우울)일수록 더 빠르게 호전되는 경향

9.5 진단 집단 비교 결과 (Table 4.8)

모수	추정값	SE	Z	\(p\)
\(\hat{\beta}_0\) (NE 절편)	22.48	0.79	28.30	<.0001
\(\hat{\beta}_1\) (NE 기울기)	−2.37	0.31	−7.59	<.0001
\(\hat{\beta}_2\) (E 절편 차이)	1.99	1.07	1.86	.063
\(\hat{\beta}_3\) (E 기울기 차이)	−0.03	0.42	−0.06	.95

우도비 검정: \(\chi^2(2) = 4.10\), \(p > .05\) → 두 집단 간 유의한 차이 없음.

내인성 집단이 기저선에서 약 2점 높은 경향이 있지만(경계적), 변화 속도는 두 집단이 동일하다. 즉, 진단 유형이 약물 치료 반응 속도를 예측하지 못한다.

10 R 코드

library(nlme)
library(lme4)

# ----------------------------------------------------------------
# Reisby 데이터 구조 시뮬레이션 (실제 데이터 대신)
# ----------------------------------------------------------------
set.seed(42)
N <- 66
n_weeks <- 6
weeks <- 0:5

# 집단 설정 (내인성 37명, 비내인성 29명)
group <- c(rep(1, 37), rep(0, 29))

# 모집단 모수 (Table 4.5 기반)
beta0 <- 23.58; beta1 <- -2.38
sigma_v0 <- sqrt(12.63); sigma_v1 <- sqrt(2.08)
rho_v    <- -1.42 / (sigma_v0 * sigma_v1)  # 상관
sigma_e  <- sqrt(12.22)

# 랜덤 효과 생성 (이변량 정규)
Sigma_v <- matrix(c(12.63, -1.42, -1.42, 2.08), 2)
library(MASS)
random_effects <- mvrnorm(N, mu=c(0,0), Sigma=Sigma_v)
v0 <- random_effects[,1]
v1 <- random_effects[,2]

# 데이터 생성 (불완전 데이터 포함)
df_list <- lapply(1:N, function(i) {
  # 일부 피험자는 관측 누락
  obs_weeks <- weeks[sample(c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE),
                             6, replace=FALSE,
                             prob=c(0.95,0.95,0.99,0.99,0.97,0.88))]
  y <- (beta0 + v0[i]) + (beta1 + v1[i]) * obs_weeks +
       rnorm(length(obs_weeks), 0, sigma_e)
  data.frame(
    id    = i,
    week  = obs_weeks,
    hdrs  = y,
    dx    = group[i]
  )
})
df <- do.call(rbind, df_list)
df$id <- factor(df$id)

cat("총 관측 수:", nrow(df), "\n")
cat("피험자 수:", length(unique(df$id)), "\n")

# ----------------------------------------------------------------
# 모형 1: 랜덤 절편 모형
# ----------------------------------------------------------------
m1 <- lmer(hdrs ~ week + (1 | id), data=df, REML=FALSE)
cat("\n--- 모형 1: 랜덤 절편 ---\n")
print(summary(m1))

# ICC 계산
var_re <- as.data.frame(VarCorr(m1))
icc <- var_re$vcov[1] / sum(var_re$vcov)
cat(sprintf("ICC = %.3f\n", icc))

# ----------------------------------------------------------------
# 모형 2: 랜덤 절편 + 기울기 모형
# ----------------------------------------------------------------
m2 <- lmer(hdrs ~ week + (week | id), data=df, REML=FALSE)
cat("\n--- 모형 2: 랜덤 절편 + 기울기 ---\n")
print(summary(m2))

# ----------------------------------------------------------------
# 모형 비교 (우도비 검정)
# ----------------------------------------------------------------
cat("\n--- 모형 비교 ---\n")
print(anova(m1, m2))

# ----------------------------------------------------------------
# 모형 3: 집단(진단) 효과 추가
# ----------------------------------------------------------------
m3 <- lmer(hdrs ~ week * dx + (week | id), data=df, REML=FALSE)
cat("\n--- 모형 3: 집단 × 시간 교호작용 ---\n")
print(summary(m3))

# 모형 2 vs 3 비교
cat("\n--- 진단 효과 유의성 검정 ---\n")
print(anova(m2, m3))

# ----------------------------------------------------------------
# 경험 베이즈 추정 (개인별 절편·기울기)
# ----------------------------------------------------------------
eb <- ranef(m2)$id
cat("\n--- 경험 베이즈 추정값 (첫 5명) ---\n")
print(head(eb, 5))

# 개인별 성장 곡선
fixef_m2 <- fixef(m2)
eb_full <- ranef(m2)$id
b0_hat <- fixef_m2["(Intercept)"] + eb_full[,"(Intercept)"]
b1_hat <- fixef_m2["week"] + eb_full[,"week"]

cat("\n개인 기울기 요약:\n")
print(summary(b1_hat))
cat(sprintf("기울기 범위: [%.2f, %.2f]\n", min(b1_hat), max(b1_hat)))

11 Python 코드

import numpy as np
import pandas as pd
from scipy import stats

# statsmodels MixedLM
import statsmodels.formula.api as smf

# ----------------------------------------------------------------
# Reisby 데이터 시뮬레이션
# ----------------------------------------------------------------
np.random.seed(42)
N = 66
weeks_all = np.arange(6)

# 집단 (내인성=1, 비내인성=0)
dx = np.array([1]*37 + [0]*29)

# 모집단 모수
beta0, beta1 = 23.58, -2.38
Sigma_v = np.array([[12.63, -1.42], [-1.42, 2.08]])
sigma_e = np.sqrt(12.22)

# 랜덤 효과 생성
ve = np.random.multivariate_normal([0, 0], Sigma_v, N)
v0, v1 = ve[:, 0], ve[:, 1]

# Long-format 데이터 생성
rows = []
for i in range(N):
    # 일부 관측 누락
    obs_mask = np.random.binomial(1, [0.95,0.95,0.99,0.99,0.97,0.88])
    for j, w in enumerate(weeks_all):
        if obs_mask[j]:
            y = (beta0 + v0[i]) + (beta1 + v1[i]) * w + \
                np.random.normal(0, sigma_e)
            rows.append({'id': i, 'week': w, 'hdrs': y, 'dx': dx[i]})

df = pd.DataFrame(rows)
print(f"총 관측 수: {len(df)}")
print(f"피험자 수: {df['id'].nunique()}")

# ----------------------------------------------------------------
# 모형 1: 랜덤 절편 모형
# ----------------------------------------------------------------
m1 = smf.mixedlm("hdrs ~ week", df, groups=df["id"])
r1 = m1.fit()
print("\n--- 모형 1: 랜덤 절편 ---")
print(r1.summary())

# ICC
var_re1 = r1.cov_re.iloc[0,0]
var_eps1 = r1.scale
icc = var_re1 / (var_re1 + var_eps1)
print(f"ICC = {icc:.3f}")
print(f"-2logL = {-2*r1.llf:.2f}")

# ----------------------------------------------------------------
# 모형 2: 랜덤 절편 + 기울기 모형
# ----------------------------------------------------------------
m2 = smf.mixedlm("hdrs ~ week", df, groups=df["id"],
                  re_formula="~week")
r2 = m2.fit()
print("\n--- 모형 2: 랜덤 절편 + 기울기 ---")
print(r2.summary())
print(f"-2logL = {-2*r2.llf:.2f}")

# ----------------------------------------------------------------
# 우도비 검정
# ----------------------------------------------------------------
lrt_stat = -2 * (r1.llf - r2.llf)
lrt_p = stats.chi2.sf(lrt_stat, df=2)
print(f"\n우도비 검정: χ²={lrt_stat:.2f}, df=2, p={lrt_p:.6f}")

# ----------------------------------------------------------------
# 모형 3: 집단 × 시간 교호작용
# ----------------------------------------------------------------
m3 = smf.mixedlm("hdrs ~ week * dx", df, groups=df["id"],
                  re_formula="~week")
r3 = m3.fit()
print("\n--- 모형 3: 집단 × 시간 교호작용 ---")
print(r3.summary())

lrt_diag = -2 * (r2.llf - r3.llf)
lrt_p_diag = stats.chi2.sf(lrt_diag, df=2)
print(f"\n진단 효과 LRT: χ²={lrt_diag:.2f}, df=2, p={lrt_p_diag:.3f}")

# ----------------------------------------------------------------
# 경험 베이즈 추정 (개인별 기울기)
# ----------------------------------------------------------------
re_effects = r2.random_effects
slopes = np.array([re_effects[i]['week'] for i in range(N)])
pop_slope = r2.params['week']
individual_slopes = pop_slope + slopes

print(f"\n개인 기울기 요약:")
print(f"  평균: {individual_slopes.mean():.2f}")
print(f"  표준편차: {individual_slopes.std():.2f}")
print(f"  범위: [{individual_slopes.min():.2f}, {individual_slopes.max():.2f}]")
print(f"  양수(비호전) 비율: {(individual_slopes > 0).mean():.1%}")

12 결과 요약

12.1 모형 비교표

모형	고정효과	랜덤효과	\(-2\log L\)	\(\Delta(-2\log L)\)	비고
랜덤 절편	절편, week	\(\sigma^2_{v_0}\)	2285.19	—	복합대칭 가정
랜덤 절편+기울기	절편, week	\(\sigma^2_{v_0}\), \(\sigma_{v_0v_1}\), \(\sigma^2_{v_1}\)	2219.04	66.15 (\(p \ll .001\))	구형성 불필요
집단×시간 추가	절편, week, DX, week×DX	동일	2214.94	4.10 (\(p=.13\))	진단 효과 없음

12.2 핵심 메시지

ANOVA/MANOVA의 제약을 하나씩 해소하는 방식으로 MRM이 구축된다:

결측 처리: MRM은 가용한 모든 데이터를 사용 → 46명이 아닌 66명 전원 분석
분산 구조: 랜덤 기울기 추가 → 복합대칭 탈피, 증가하는 분산 포착
개인 추세: 경험 베이즈 → 집단 평균이 아닌 각 피험자의 궤적 추정
공변량: 시불변(집단)·시변(혈장 농도 등) 공변량을 모두 통합

13 관련 주제

선행 지식

후속 주제

다항식 혼합효과 모형 (Ch.5, 미작성) — 비선형 시간 추세
공분산 패턴 모형 (Ch.6, 미작성) — 랜덤 효과 없는 유연한 공분산 구조
자기상관 오차 MRM — 랜덤 효과 + 시계열 오차 결합

관련 개념

REML — 제한 최대가능도 추정
GEE 모형 — 주변 모형, 랜덤 효과 없는 대안