Kwangmin Kim - 랜덤 절편 MRM — 모형 구조, 불완전 데이터, 복합대칭, 추정과 검정

1 출발점 요약: 왜 랜덤 절편이 필요한가

이전 포스트에서 단순 선형회귀 \(y_{ij} = \beta_0 + \beta_1 t_{ij} + \varepsilon_{ij}\) 가 종단 데이터에서 두 가지 이유로 부적합함을 확인했다.

독립 가정 위반: 같은 피험자의 반복 측정은 상관됨 → OLS 표준오차 과소 추정
이질성 무시: \(\beta_0, \beta_1\) 에 \(i\) 첨자 없음 → 모든 사람이 동일한 절편·기울기

가장 단순한 해결책은 모형에 피험자별 개인 효과를 하나 추가하는 것이다.

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

이것이 랜덤 절편 MRM이다 (식 4.2, Hedeker & Gibbons, 2006, Ch.4).

2 모형 구조: 각 항목의 역할

2.1 새로 추가된 항 \(v_{0i}\) 란 무엇인가

\(v_{0i}\) 는 피험자 \(i\) 의 개인 수준 효과(individual-specific effect)다. 이것이 모형에 추가된 의미를 단계별로 이해한다.

직관: 피험자 \(i\) 의 관측값들이 집단 평균 추세 \(\beta_0 + \beta_1 t\) 에서 얼마나 벗어나 있는지를 나타내는 상수 편차다. 어떤 환자는 집단 평균보다 항상 우울증 점수가 높고, 어떤 환자는 항상 낮다. 이 개인별 수준 차이가 \(v_{0i}\) 다.

\(v_{0i}\) 의 두 가지 해석

편차로서의 해석: \(v_{0i}\) 는 피험자 \(i\) 가 집단 평균 절편 \(\beta_0\) 에서 얼마나 벗어나 있는지를 나타낸다. \(v_{0i} > 0\) 이면 평균보다 높은 기저 수준, \(v_{0i} < 0\) 이면 낮은 기저 수준.

확률 변수로서의 해석: 연구 참여자들은 더 큰 모집단(population)에서 무작위 표본이다. 따라서 \(v_{0i}\) 도 이 모집단에서 추출된 무작위 값이다 → 랜덤 효과(random effect). 고정 효과(fixed effect)가 아닌 이유다.

가장 많이 사용되는 분포 가정: \[v_{0i} \overset{iid}{\sim} N(0, \sigma^2_v)\]

평균 0: 집단적으로 모집단 평균 \(\beta_0\) 에서 편향 없이 분포
분산 \(\sigma^2_v\): 개인 간 기저 수준의 이질성 크기

2.2 분포 가정 정리

\[v_{0i} \overset{iid}{\sim} N(0, \sigma^2_v) \quad \text{— 피험자 간 독립, 동일 분포}\]

\[\varepsilon_{ij} \mid v_{0i} \overset{iid}{\sim} N(0, \sigma^2) \quad \text{— } v_{0i} \text{에 조건부로 독립}\]

\[v_{0i} \perp \varepsilon_{ij} \quad \text{— 랜덤 효과와 잔차는 독립}\]

\(\varepsilon_{ij}\) 의 조건부 독립이 핵심이다. 단순 회귀에서 \(\varepsilon_{ij}\) 들의 상관은 숨어있는 \(v_{0i}\) 때문에 발생했다. \(v_{0i}\) 를 모형에 명시적으로 포함시키면, \(v_{0i}\) 가 주어진 조건에서 잔차들은 진짜로 독립이 된다.

3 위계적 표현: 수준-1과 수준-2

같은 모형을 두 층위로 분해하면 구조가 더 명확해진다.

3.1 수준-1 모형 (피험자 내, Within-Subject) — 식 4.3

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

\(b_{0i}\): 피험자 \(i\) 의 개인 절편 (피험자마다 다름)
\(b_{1i}\): 피험자 \(i\) 의 개인 기울기 (지금은 모두 동일)
\(\varepsilon_{ij}\): 이 개인 추세 주변의 무작위 변동

직관: “피험자 \(i\) 의 관측값은 그 사람만의 직선(\(b_{0i}\), \(b_{1i}\)) 주변에 분포한다.”

3.2 수준-2 모형 (피험자 간, Between-Subject) — 식 4.4

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

절편: 집단 평균 \(\beta_0\) + 개인 편차 \(v_{0i}\)
기울기: 모든 사람에게 동일한 집단 기울기 \(\beta_1\)

직관: “각 사람의 출발점(\(b_{0i}\))은 다르지만, 변화 속도(\(b_{1i}\))는 모두 같다.”

3.3 합치면 단일 방정식 (식 4.2)

수준-1에 수준-2를 대입하면:

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

위계적 표현의 가치는 해석의 편의에 있다. “수준-2 공변량(피험자 수준 변수)이 절편 \(b_{0i}\) 와 기울기 \(b_{1i}\) 를 어떻게 설명하는가”라는 질문을 자연스럽게 확장할 수 있다. 이것을 “기울기를 결과변수로”(slopes-as-outcomes) 모형이라고도 부른다.

시각적 이해: 평행선 다발

랜덤 절편 모형에서 각 피험자의 추세선은:

\[\hat{y}_{ij} = (\beta_0 + v_{0i}) + \beta_1 t_{ij}\]

기울기는 모두 동일: \(\beta_1\) (모든 선이 평행)
절편만 다름: \(\beta_0 + v_{0i}\) (선들이 평행하게 위아래로 분산)

\(\sigma^2_v\) 가 크면 선들이 넓게 퍼지고, \(\sigma^2_v \approx 0\) 이면 모든 선이 거의 겹친다.

집단 평균 추세 \(\beta_0 + \beta_1 t\) 와 각 개인 추세의 차이는 시간에 무관하게 상수 \(v_{0i}\) 로 유지된다 — 이것이 “랜덤 절편” 모형이 함의하는 평행성이다.

4 불완전 데이터 구조 (Incomplete Data Across Time)

이것이 MRM의 가장 중요한 실용적 장점이며, ANOVA/MANOVA와 근본적으로 다른 지점이다.

4.1 불완전 데이터의 세 가지 유형

종단 연구에서 “완전 데이터”는 오히려 예외다. 세 가지 불완전 상황이 흔히 발생한다.

유형	설명	예시
불균등 관측 횟수	피험자마다 \(n_i\) 가 다름	어떤 환자는 6번, 어떤 환자는 3번 측정
단속적 결측 (Intermittent)	중간 시점에서 결측, 이후 재참여	Week 2 빠지고 Week 3~5 있음
탈락 (Attrition)	특정 시점 이후 완전 탈락	Week 3 이후 데이터 없음

모형 수식에서 이것이 어떻게 표현되는지 다시 보자:

\[y_{ij}, \quad j = 1, \ldots, n_i\]

\(n\) 에 \(i\) 첨자가 붙어 있다. 피험자마다 관측 횟수가 다를 수 있다는 것이 모형에 내재되어 있다. 또한 시간 변수 \(t_{ij}\) 에도 \(i\) 첨자가 있으므로 피험자마다 측정 시점이 달라도 된다.

ANOVA/MANOVA와의 대비:

방법	결측 처리
반복측정 ANOVA	구형성 가정 + 완전 데이터 필요 (결측 피험자 제외)
MANOVA	완전 데이터만 사용 (목록 삭제, listwise deletion)
MRM	모든 피험자의 가용 데이터를 모두 사용

Reisby 데이터에서 66명 중 완전 데이터가 있는 피험자는 46명(70%)이다. MANOVA는 46명만 분석하지만 MRM은 66명 전체를 사용한다.

4.2 결측 데이터 메커니즘 (Missing Data Mechanism)

MRM이 결측 데이터를 “자동으로” 처리한다는 것은 어떤 조건에서 성립하는가? 결측 메커니즘을 이해해야 한다.

Rubin(1976)의 분류체계:

결측 메커니즘 3가지

MCAR (Missing Completely At Random, 완전 무작위 결측)

결측이 어떤 변수와도 무관하다. 동전 던지기처럼 순전히 무작위.

\[P(\text{결측} \mid y_{\text{관측}}, y_{\text{결측}}, X) = P(\text{결측})\]

MAR (Missing At Random, 무작위 결측)

결측이 관측된 데이터에 의존할 수 있으나, 결측된 값 자체에는 의존하지 않는다.

\[P(\text{결측} \mid y_{\text{관측}}, y_{\text{결측}}, X) = P(\text{결측} \mid y_{\text{관측}}, X)\]

예: 우울증이 심한 사람(높은 초기 HDRS)이 중도 탈락할 가능성이 높다면, 그 결측은 관측된 초기 값에 의존하므로 MAR.

MNAR (Missing Not At Random, 비무작위 결측)

결측이 결측된 값 자체에 의존한다.

\[P(\text{결측} \mid y_{\text{결측}}) \neq \text{const}\]

예: 현재의 우울증 점수가 높아서 측정을 거부하는 경우. 측정되지 않은 그 값이 결측 여부를 결정한다.

MRM의 결측 처리 능력:

Laird(1988)는 MRM을 ML로 추정하면 MAR 가정 아래에서 유효한 통계적 추론을 제공함을 보였다.

\[\text{(MRM + ML 추정)} + \text{MAR 가정} \Rightarrow \text{편향 없는 추정}\]

이것이 왜 성립하는가? ML 추정에서 가능도 함수는 각 피험자가 실제로 관측한 데이터만을 사용하여 구성된다. 관측된 데이터의 가능도를 최대화할 때, MAR 메커니즘은 무시 가능(ignorable)하다 — 결측 발생 과정을 별도로 모형화하지 않아도 올바른 추론이 가능하다.

왜 MAR이 MNAR보다 안전한가

MNAR 상황에서는 MRM도 편향된다. 예를 들어:

치료 효과가 없는(개선 없는) 환자들이 탈락할 때
탈락 이유가 측정되지 않은 현재 상태(악화)에 있을 때

이 경우 남은 표본은 치료에 반응한 사람들로 편향된다. MRM은 이 선택 편향을 자동으로 교정하지 못한다.

실용적 원칙: 결측 이유가 관측된 공변량(진단군, 기저선 점수 등)과 관측된 결과값으로 설명 가능하면 MAR으로 간주할 수 있다. 탈락 이유가 “지금 상태가 너무 나빠서” 처럼 결측 시점의 측정되지 않은 값에 직접 연결되면 MNAR을 의심해야 한다.

4.3 탈락 vs 단속적 결측

탈락(attrition)과 단속적 결측(intermittent missing data)은 발생 메커니즘이 다를 수 있다.

유형	특성	MAR 가정 타당성
탈락	Week \(k\) 이후 영구 결측	MAR 가정이 더 강하게 요구됨 — 탈락 이유가 이후 잠재적 결과와 관련될 수 있음
단속적 결측	특정 시점 결측 후 재참여	MAR 가정이 더 합리적 — 일시적 사정(이동, 방문 일정) 때문인 경우가 많음

MRM은 두 유형 모두 원칙적으로 처리 가능하나, 탈락에서 MAR 가정의 검토가 더 중요하다.

4.4 코드: 불완전 데이터 구조 확인

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 불완전 데이터 시뮬레이션 예시
np.random.seed(42)
N = 10          # 피험자 수
max_T = 6       # 최대 시점 수

# 각 피험자의 관측 횟수를 다르게 설정
n_obs = np.random.randint(3, 7, size=N)  # 피험자마다 3~6회 측정

data = []
for i in range(N):
    # 피험자 i의 관측 시점 (랜덤 선택)
    times = sorted(np.random.choice(range(max_T), size=n_obs[i], replace=False))
    for j, t in enumerate(times):
        data.append({'id': i, 'week': t, 'obs_num': j + 1})

df = pd.DataFrame(data)
print("불완전 데이터 구조:")
print(df.pivot_table(index='id', columns='week',
                     values='obs_num', fill_value=np.nan))

# 각 피험자의 관측 횟수
print(f"\n피험자별 관측 횟수: {dict(df.groupby('id').size())}")

# Reisby 데이터의 불완전 구조 확인
# 교재 Table 4.1 기반
weeks = [0, 1, 2, 3, 4, 5]
ns    = [61, 63, 65, 65, 63, 58]

print("Reisby 데이터 시점별 피험자 수")
print(f"전체 피험자: 66명")
print(f"완전 데이터:  46명 ({46/66*100:.0f}%)")
print(f"결측 피험자:  {66-46}명 ({(66-46)/66*100:.0f}%)")
print()
for w, n in zip(weeks, ns):
    missing = 66 - n
    print(f"Week {w}: n={n:2d}, 결측={missing:2d}명 ({missing/66*100:.0f}%)")

# R: lme4로 불완전 데이터 자동 처리 확인
library(lme4)
library(dplyr)

# 완전 데이터와 불완전 데이터의 피험자 수
n_complete <- sum(tapply(reisby_long$hdrs, reisby_long$id,
                         function(x) !any(is.na(x))))
n_total    <- length(unique(reisby_long$id))

cat("완전 데이터 피험자:", n_complete, "\n")
cat("전체 피험자:", n_total, "\n")

# lme4는 가용 데이터를 모두 사용 (별도 처리 불필요)
fit_ri <- lmer(hdrs ~ week + (1 | id), data = reisby_long, REML = FALSE)
# 모형 적합 시 n_total명의 데이터가 모두 사용됨

5 복합대칭 공분산 구조와 ICC

5.1 분산-공분산 구조 유도 (식 4.5)

랜덤 절편 모형 \(y_{ij} = \beta_0 + \beta_1 t_{ij} + v_{0i} + \varepsilon_{ij}\) 에서 관측값의 분산과 공분산을 계산한다.

분산 (같은 피험자, 같은 시점):

\[V(y_{ij}) = V(v_{0i} + \varepsilon_{ij}) = \sigma^2_v + \sigma^2\]

\(v_{0i}\) 와 \(\varepsilon_{ij}\) 가 독립이므로 분산이 더해진다
모든 피험자, 모든 시점에서 동일하다 (시점 불변, 피험자 불변)

공분산 (같은 피험자, 다른 시점):

\[C(y_{ij}, y_{ij'}) = C(v_{0i} + \varepsilon_{ij},\; v_{0i} + \varepsilon_{ij'})\] \[= V(v_{0i}) + C(\varepsilon_{ij}, \varepsilon_{ij'}) = \sigma^2_v + 0 = \sigma^2_v\]

\(v_{0i}\) 가 두 시점에 공통으로 들어가므로 공분산이 \(\sigma^2_v\) 로 고정된다
모든 시점 쌍에서 동일 (시차에 무관)

이 두 특성을 합치면 복합대칭(Compound Symmetry, CS) 구조가 된다.

\[V_i = \begin{bmatrix} \sigma^2_v + \sigma^2 & \sigma^2_v & \cdots & \sigma^2_v \\ \sigma^2_v & \sigma^2_v + \sigma^2 & \cdots & \sigma^2_v \\ \vdots & & \ddots & \vdots \\ \sigma^2_v & \sigma^2_v & \cdots & \sigma^2_v + \sigma^2 \end{bmatrix}\]

복합대칭의 강한 가정

복합대칭은 두 가지를 동시에 가정한다:

등분산: 모든 시점의 분산이 \(\sigma^2_v + \sigma^2\) 로 동일
등공분산: 모든 시점 쌍의 공분산이 \(\sigma^2_v\) 로 동일 (→ 시차와 무관하게 같은 상관)

Reisby 데이터에서 분산은 시간에 따라 증가(4.53→7.22)하고, 상관은 시차가 클수록 작아진다. 두 가정 모두 위반된다. 따라서 랜덤 절편 모형은 이 데이터의 출발점이지, 최종 모형이 아니다.

5.2 급내상관계수 (ICC)

분산과 공분산을 상관으로 표현하면:

\[\text{ICC} = \rho = \frac{C(y_{ij}, y_{ij'})}{V(y_{ij})} = \frac{\sigma^2_v}{\sigma^2_v + \sigma^2}\]

해석: ICC는 총 분산 중 피험자 수준에 귀속되는 비율이다.

ICC 값	의미
ICC \(\approx 0\)	개인 간 이질성이 거의 없음. 반복 측정이 독립에 가까움 → OLS도 무방
ICC \(= 0.46\) (Reisby)	총 변동의 46%가 피험자 간 차이. 반복 측정 간 실질적 상관 → MRM 필수
ICC \(\approx 1\)	거의 모든 변동이 피험자 수준. 한 번만 측정해도 충분

ICC = 0.46 의 실제적 의미 (Reisby 데이터)

랜덤 절편 모형에서 \(\hat{\sigma}^2_v = 16.15\), \(\hat{\sigma}^2 = 19.04\) 이므로:

\[\text{ICC} = \frac{16.15}{16.15 + 19.04} = 0.46\]

같은 피험자의 임의의 두 시점 HDRS 점수 간 예측 상관이 0.46이다. 이것이 의미하는 것은:

피험자 간 기저 우울증 수준의 차이가 총 변동의 46%를 설명한다
랜덤 절편을 무시하는 OLS는 이 46%를 오차로 처리한다
OLS 유효 표본 크기: \(n_{\text{eff}} = 6 / (1 + 5 \times 0.46) \approx 1.7\) → 6번 측정이 실제로는 약 1.7번의 정보량

이것이 OLS 표준오차가 MRM보다 작게 나오는 이유다.

6 분산 분해: OLS vs 랜덤 절편 MRM

Reisby 데이터에 단순 OLS와 랜덤 절편 MRM을 각각 적합하면 결과가 어떻게 다른가.

6.1 OLS 결과 (집단 내 군집 무시)

\[\hat{\beta}_0 = 23.60 \;(SE = 0.55), \quad \hat{\beta}_1 = -2.41 \;(SE = 0.18), \quad \hat{\sigma}^2 = 35.40\]

OLS는 단 하나의 분산 \(\hat{\sigma}^2 = 35.40\) 에 모든 변동을 집어넣는다.

6.2 랜덤 절편 MRM 결과 (Table 4.3, Hedeker & Gibbons, 2006)

\[\hat{\beta}_0 = 23.55 \;(SE = 0.64), \quad \hat{\beta}_1 = -2.38 \;(SE = 0.14)\] \[\hat{\sigma}^2_v = 16.15 \;(SE = 3.41), \quad \hat{\sigma}^2 = 19.04 \;(SE = 1.53)\] \[-2\log L = 2285.19\]

MRM은 OLS의 단일 오차 분산 \(35.40\) 을 두 성분으로 분해한다:

\[\underbrace{35.40}_{\text{OLS 오차 분산}} \approx \underbrace{16.15}_{\text{피험자 간 분산 } \hat{\sigma}^2_v} + \underbrace{19.04}_{\text{피험자 내 분산 } \hat{\sigma}^2}\]

“한 통계학자의 오차가 다른 통계학자의 연구 주제다”

Hedeker & Gibbons(2006, p.56)는 이 분해를 이렇게 표현한다:

“What the ordinary regression model lumped together into error variance (35.40), the random-intercept model separates into within-subjects and between-subjects variances (19.04 and 16.15, respectively). This illustrates a golden rule of statistics: One statistician’s error term is another’s career!”

OLS에서는 피험자 간 이질성 \(\sigma^2_v = 16.15\) 가 오차로 처리된다. MRM에서는 이것이 명시적으로 추정되어 해석 가능한 모수가 된다.

6.3 두 분석의 계수 비교

모수	OLS	랜덤 절편 MRM
\(\hat{\beta}_0\)	23.60	23.55
\(\hat{\beta}_1\)	-2.41	-2.38
\(SE(\hat{\beta}_0)\)	0.55	0.64
\(SE(\hat{\beta}_1)\)	0.18	0.14
\(\hat{\sigma}^2\)	35.40	19.04 (피험자 내)
피험자 간 분산	—	16.15

계수 추정값은 매우 비슷하다. 점추정의 방향과 크기는 OLS와 MRM 모두 비슷하게 나온다. 차이는 표준오차다. \(SE(\hat{\beta}_0)\) 는 MRM에서 0.64로 OLS보다 크고, \(SE(\hat{\beta}_1)\) 은 0.14로 OLS보다 작다. 이 방향이 직관과 반대처럼 보일 수 있다.

왜 기울기 SE는 MRM에서 더 작은가

\(\beta_0\) (절편)는 피험자 간 비교에 의존 → 군집을 무시하면 SE 과소 추정 \(\beta_1\) (기울기)는 피험자 내 시간 변화에 의존 → OLS가 피험자 내 오차 \(\sigma^2\) 와 피험자 간 오차 \(\sigma^2_v\) 를 합산한 큰 분산을 사용하므로 기울기 SE가 오히려 과대 추정될 수 있다.

MRM은 기울기를 피험자 내 변화로부터 추정하므로, 피험자 간 이질성(\(\sigma^2_v\))이 기울기 추정의 불확실성에 기여하지 않는다. 따라서 기울기 SE가 OLS보다 작게 나온다.

7 추정: ML 추정과 가능도 함수

7.1 주변 가능도 (Marginal Likelihood)

MRM에서 \(v_{0i}\) 는 관측되지 않는 잠재 변수다. 추정을 위해 \(v_{0i}\) 를 적분하여 소거한다.

\(y_i = (y_{i1}, \ldots, y_{in_i})'\) 의 주변(marginal) 분포:

\[y_i \sim N(X_i \beta, \; V_i), \quad V_i = \sigma^2_v \mathbf{1}_{n_i} \mathbf{1}_{n_i}' + \sigma^2 I_{n_i}\]

여기서 \(X_i = [1, t_{i1}; 1, t_{i2}; \ldots; 1, t_{in_i}]\) 이고 \(\mathbf{1}_{n_i}\) 는 \(n_i \times 1\) 벡터.

주변 로그 가능도:

\[\ell(\beta, \sigma^2_v, \sigma^2) = -\frac{1}{2} \sum_{i=1}^N \left[\log |V_i| + (y_i - X_i\beta)' V_i^{-1} (y_i - X_i\beta)\right]\]

핵심: \(V_i\) 의 차원이 피험자마다 \(n_i \times n_i\) 로 다를 수 있다. 이것이 불완전 데이터를 자연스럽게 처리하는 수학적 메커니즘이다. 관측 횟수가 적은 피험자는 작은 \(V_i\) 를, 많은 피험자는 큰 \(V_i\) 를 사용하여 가능도가 계산된다.

8 가설 검정

8.1 고정효과 \(\beta\): Wald 검정

\[Z = \frac{\hat{\beta}_k}{SE(\hat{\beta}_k)} \sim N(0,1)\]

고정효과에 대해서는 Wald 검정이 적합하다. 단, 표본이 작을 때는 \(t\) 분포를 사용하는 소프트웨어도 있다(자유도 계산은 소프트웨어마다 다름).

Reisby 결과에서:

\(\hat{\beta}_1 = -2.38\), \(SE = 0.14\) → \(Z = -17.0\), \(p < 0.0001\)
“환자들의 주당 HDRS 감소율이 0과 유의하게 다르다” = 치료 기간에 걸쳐 유의한 개선

8.2 분산 성분 \(\sigma^2_v\): Wald 검정 부적합

분산은 0 이상의 값만 가질 수 있다는 경계값 제약 때문에 Wald 검정을 사용하지 않는다.

\[\sigma^2_v \geq 0 \quad \text{(경계값 0에서 정규 근사 성립 안 함)}\]

경계값 문제 (Boundary Problem)

\(H_0: \sigma^2_v = 0\) 검정을 Wald 검정으로 하면 \(\hat{\sigma}^2_v / SE(\hat{\sigma}^2_v)\) 를 표준 정규분포와 비교한다.

그러나 귀무가설이 경계값(boundary) \(\sigma^2_v = 0\) 에서 성립하므로, 검정 통계량의 표집 분포가 정규 분포를 따르지 않는다. 정규 분포를 가정하면 \(p\)-값이 보수적으로(크게) 추정된다 — 실제보다 랜덤 효과를 덜 발견하게 된다.

분산 성분에는 우도비 검정(LRT)을 사용한다 — 단, \(p\)-값을 2로 나누는 보정을 적용한다.

\[\text{LRT 통계량} = -2\log L_{\text{제한}} - (-2\log L_{\text{확장}}) \sim \chi^2_{df}\]

경계값 문제로 인해 얻은 \(\chi^2\) \(p\)-값을 2로 나눈다:

\[p_{\text{보정}} = p_{\chi^2} / 2\]

적용: 랜덤 절편 모형(\(\sigma^2_v\) 포함) vs 단순 OLS(\(\sigma^2_v = 0\))를 LRT로 비교. 자유도 = 1 (추가 모수 \(\sigma^2_v\) 1개).

8.3 중첩 모형 비교 (Nested Model Comparison)

랜덤 절편 모형에서 랜덤 절편+기울기 모형으로 확장할 때도 LRT를 사용한다:

랜덤 절편: \(-2\log L = 2285.19\) (모수: \(\beta_0, \beta_1, \sigma^2_v, \sigma^2\) — 4개)
랜덤 절편+기울기: \(-2\log L = 2219.04\) (추가 모수: \(\sigma^2_{v_1}, \sigma_{v_0 v_1}\) — 6개)

\[\text{LRT} = 2285.19 - 2219.04 = 66.15, \quad df = 2, \quad p < 0.0001\]

랜덤 기울기를 추가함으로써 모형이 유의하게 개선된다.

ML vs REML: 모형 비교 시 주의

고정효과 구조가 다른 두 모형 비교: ML 사용
분산 성분 구조가 다른 두 모형 비교 (같은 고정효과): REML 사용

REML 추정값으로 구성된 \(-2\log L\) 은 고정효과 구조가 다른 두 모형 비교에 사용할 수 없다. REML 가능도는 고정효과를 소거한 잔차의 가능도이므로, 고정효과 구조가 다르면 비교 대상이 달라진다.

9 Reisby 데이터 분석: 코드와 해석

9.1 Python 구현 (statsmodels)

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

# --- 데이터 로드 (long format) ---
# df 컬럼: id (피험자), week (0~5), hdrs (HDRS 점수), dx (진단: 0=비내인성, 1=내인성)

# --- Model 1: OLS (독립 가정 위반) ---
m_ols = smf.ols("hdrs ~ week", data=df).fit()

# --- Model 2: 랜덤 절편 MRM ---
m_ri = smf.mixedlm("hdrs ~ week",
                   data=df,
                   groups=df["id"]).fit(reml=False)  # ML 추정

print("=== OLS 결과 ===")
print(f"β_0 = {m_ols.params['Intercept']:.3f} (SE = {m_ols.bse['Intercept']:.3f})")
print(f"β_1 = {m_ols.params['week']:.3f}   (SE = {m_ols.bse['week']:.3f})")
print(f"σ²  = {m_ols.mse_resid:.3f}")

print("\n=== 랜덤 절편 MRM 결과 ===")
print(f"β_0 = {m_ri.params['Intercept']:.3f} (SE = {m_ri.bse['Intercept']:.3f})")
print(f"β_1 = {m_ri.params['week']:.3f}   (SE = {m_ri.bse['week']:.3f})")
print(f"σ²_v = {m_ri.cov_re.values[0][0]:.3f}")
print(f"σ²   = {m_ri.scale:.3f}")
print(f"ICC  = {m_ri.cov_re.values[0][0] / (m_ri.cov_re.values[0][0] + m_ri.scale):.3f}")
print(f"-2logL = {-2 * m_ri.llf:.2f}")

# --- LRT: 랜덤 절편 필요성 검정 ---
from scipy import stats

# OLS의 -2logL (독립 오차만 있는 경우)
ll_ols = m_ols.llf
ll_ri  = m_ri.llf

lrt_stat = -2 * (ll_ols - ll_ri)
# 경계값 보정: p-값을 2로 나눔
p_raw    = stats.chi2.sf(lrt_stat, df=1)
p_adj    = p_raw / 2

print(f"\n=== LRT: σ²_v = 0 검정 ===")
print(f"LRT 통계량 = {lrt_stat:.3f}")
print(f"χ² p-값   = {p_raw:.6f}")
print(f"보정 p-값 = {p_adj:.6f}  (경계값 보정)")

9.2 R 구현 (lme4)

library(lme4)
library(lmerTest)

# --- Model 1: OLS ---
m_ols <- lm(hdrs ~ week, data = reisby_long)

# --- Model 2: 랜덤 절편 MRM (ML 추정) ---
m_ri  <- lmer(hdrs ~ week + (1 | id),
              data = reisby_long,
              REML = FALSE)

# 결과 요약
summary(m_ri)

# ICC 계산
vc <- as.data.frame(VarCorr(m_ri))
sigma2_v <- vc[vc$grp == "id",    "vcov"]
sigma2_e <- vc[vc$grp == "Residual", "vcov"]
icc <- sigma2_v / (sigma2_v + sigma2_e)
cat("ICC =", round(icc, 3), "\n")

# --- LRT: 랜덤 절편 필요성 검정 ---
# OLS를 lme4 형식으로 (랜덤 효과 없음)
m_null <- lm(hdrs ~ week, data = reisby_long)

# 직접 LRT
lrt_stat <- -2 * (logLik(m_null) - logLik(m_ri))
# 경계값 보정
p_adj <- pchisq(lrt_stat, df = 1, lower.tail = FALSE) / 2
cat("LRT =", round(lrt_stat, 3), ", 보정 p-값 =", round(p_adj, 6), "\n")

9.3 집단 평균 추세 시각화

weeks_grid = np.linspace(0, 5, 100)
beta0, beta1 = m_ri.params['Intercept'], m_ri.params['week']

# 집단 평균 추세
mean_trend = beta0 + beta1 * weeks_grid

# 개인별 랜덤 절편 추출
re = m_ri.random_effects  # {id: 절편 편차}

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# 왼쪽: 개인별 추세 (스파게티 플롯)
ax = axes[0]
for pid, effect in re.items():
    ind_intercept = beta0 + effect[0]
    ax.plot(weeks_grid, ind_intercept + beta1 * weeks_grid,
            alpha=0.3, color='steelblue', linewidth=0.8)
ax.plot(weeks_grid, mean_trend, 'r-', linewidth=2.5, label='집단 평균 추세')
ax.set_xlabel("Week")
ax.set_ylabel("Estimated HDRS")
ax.set_title("랜덤 절편 MRM:\n개인 추세 (평행선 다발)")
ax.legend()

# 오른쪽: 관측값 vs 예측값
ax = axes[1]
preds = m_ri.fittedvalues
ax.scatter(preds, df['hdrs'], alpha=0.3, s=15)
ax.plot([5, 28], [5, 28], 'r--', label='y = x')
ax.set_xlabel("Fitted values")
ax.set_ylabel("Observed HDRS")
ax.set_title("관측값 vs 예측값")
ax.legend()

plt.tight_layout()
plt.show()

10 모형의 한계와 다음 단계

랜덤 절편 모형은 MRM의 출발점이지 최종 모형이 아니다.

한계	증거	해결책
복합대칭 가정 위반	상관이 시차에 따라 감소 (0.49 → 0.18)	랜덤 기울기 추가 또는 공분산 패턴 모형
분산 증가 포착 불가	sd 4.53 → 7.22로 증가	랜덤 기울기 추가 → 시간에 따른 분산 증가 모형화
개인별 변화율 동일 가정	일부는 급격 호전, 일부는 무반응	랜덤 기울기 \(v_{1i}\) 추가

이 한계들은 다음 포스트 — 랜덤 절편+기울기 MRM에서 다룬다.

11 관련 주제

선행 지식

후속 주제

관련 개념