Kwangmin Kim - 단순 선형회귀 모형 — 혼합효과 회귀모형의 출발점

1 왜 출발점이 필요한가

혼합효과 회귀모형(Mixed-effects Regression Model, MRM)은 처음 마주치면 복잡하게 느껴진다. 랜덤 효과, 공분산 구조, 경험 베이즈 추정 등 낯선 개념이 한꺼번에 등장하기 때문이다.

그러나 MRM은 전혀 새로운 모형이 아니다. 익숙한 단순 선형회귀 모형을 종단 데이터에 적용할 때 발생하는 두 가지 근본적 문제를 해결하기 위해 자연스럽게 확장된 모형이다. 따라서 MRM을 이해하는 가장 효율적인 경로는 단순 선형회귀의 가정들을 하나하나 해부하고, 어디가 문제이며, 어떻게 고쳐야 하는지를 추적하는 것이다 (Hedeker & Gibbons, 2006, Ch.4).

2 표기 체계: 이중 첨자의 의미

수식에 들어가기 전에, 종단 데이터에서 사용하는 이중 첨자 표기를 명확히 이해해야 한다.

\[y_{ij}\]

\(i = 1, 2, \ldots, N\): 피험자 번호. \(N\) 은 전체 피험자 수.
\(j = 1, 2, \ldots, n_i\): 피험자 \(i\) 의 관측 시점 번호. \(n_i\) 는 피험자 \(i\) 가 측정된 총 횟수.

왜 \(n_i\) 에 \(i\) 첨자가 붙는가

횡단 데이터에서는 모든 개인이 한 번씩 측정된다. 종단 데이터에서는 피험자마다 측정 횟수가 다를 수 있다 — 어떤 사람은 6번 측정되고, 어떤 사람은 탈락하여 3번만 측정된다. \(n_i\) 는 이 불균형을 수용하는 핵심 표기다.

이것이 MRM이 결측 데이터를 자연스럽게 처리하는 이유다. \(n_i\) 가 피험자마다 달라도 모형이 성립한다.

시간 변수 \(t_{ij}\) 도 마찬가지다. 피험자 \(i\) 가 시점 \(j\) 에서 측정된 시간값 — 예를 들어 “0주, 1주, 2주”처럼 주 단위 값 — 이며, 피험자마다 측정 시점이 달라도 된다. 이것이 ANOVA/MANOVA와의 중요한 차이다. ANOVA/MANOVA에서 시간은 범주형 변수이고 모든 피험자가 동일한 시점에 측정되어야 한다. MRM에서 시간은 연속형 변수이고 피험자별로 다른 시점이 허용된다.

3 단순 선형회귀 모형 (식 4.1)

종단 데이터에 단순 선형회귀를 적용하면 다음과 같다:

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

첨자를 무시하면 이것은 단순히 \(y = \beta_0 + \beta_1 t + \varepsilon\) 이다 — 결과 변수를 시간에 회귀시키는 일반 선형회귀. 이중 첨자는 “누구의 측정값인지(\(i\))”와 “언제 측정된 것인지(\(j\))”를 추적할 뿐이다.

3.1 고정 계수의 의미

모수	의미	직관
\(\beta_0\)	절편	시간이 0일 때의 평균 \(y\) 값
\(\beta_1\)	기울기	시간이 1단위 증가할 때 \(y\) 의 평균 변화량

핵심 주의사항: \(\beta_0\) 와 \(\beta_1\) 에는 \(i\) 첨자가 없다. 즉, 모든 피험자가 동일한 절편과 기울기를 공유한다. 이 모형에서 피험자 \(i=1\) 과 \(i=100\) 은 정확히 같은 기저 수준 \(\beta_0\) 와 같은 변화 속도 \(\beta_1\) 을 가진다는 것이다.

이것이 첫 번째 문제의 씨앗이다. 현실에서 모든 사람이 같은 속도로 변화한다는 가정은 얼마나 현실적인가?

직관: 평행선 가정

\(\beta_0\), \(\beta_1\) 이 모든 사람에게 동일하다면, 각 피험자의 예측 추세선 \(\hat{y}_{ij} = \beta_0 + \beta_1 t_{ij}\) 는 모든 사람에 대해 정확히 같은 직선이다.

실제 데이터에서 관측값이 이 직선에서 벗어나는 이유는 오직 \(\varepsilon_{ij}\) (무작위 오차)뿐이다. 그러나 현실에서 “우울증 환자 A는 치료에 빠르게 반응하고 B는 전혀 반응하지 않는다”는 상황은 오차가 아닌 체계적인 개인 간 이질성이다. 이것을 오차로 처리하는 것은 구조적 신호를 소음으로 오해하는 것이다.

4 iid 오차 가정의 해부

단순 선형회귀의 핵심 가정은 다음과 같다:

\[\varepsilon_{ij} \overset{iid}{\sim} N(0, \sigma^2)\]

“iid”는 세 가지 가정의 압축이다.

4.1 가정 1: 정규성 (Normality)

\[\varepsilon_{ij} \sim N(\cdot, \cdot)\]

오차가 정규 분포를 따른다는 가정이다. 중간 범위 관측값에서 작은 오차, 극단값에서 큰 오차가 나타날 확률을 정규 분포가 표현한다.

직관: 측정 오차는 수많은 독립적인 소규모 원인들(측정 기기 오차, 피험자 상태 변동 등)의 합으로 발생한다. 중심극한정리에 의해 이 합은 정규 분포에 근접한다. 종단 데이터에서 정규성 가정은 대표본에서 비교적 안전하다.

4.2 가정 2: 등분산성 (Homoscedasticity)

\[V(\varepsilon_{ij}) = \sigma^2 \quad \text{(모든 } i, j \text{에 대해 동일)}\]

분산이 피험자(\(i\))와 시점(\(j\))에 무관하게 일정하다. 측정 오차의 크기가 누가 측정되든, 언제 측정되든 항상 같다는 가정이다.

직관: Reisby 데이터에서 Week 0의 표준편차는 4.53인 반면, Week 5에서는 7.22이다. 치료가 진행될수록 개인 간 반응 이질성이 커진다 — 분산이 시간에 따라 증가하는 이분산성이 나타난다. 등분산 가정이 이미 위반된 것이다.

4.3 가정 3: 독립성 (Independence) — 가장 중요

\[\text{Corr}(\varepsilon_{ij}, \varepsilon_{ij'}) = 0 \quad (j \neq j')\] \[\text{Corr}(\varepsilon_{ij}, \varepsilon_{i'j}) = 0 \quad (i \neq i')\]

오차들이 서로 상관이 없다. 특히 같은 피험자(\(i\))의 서로 다른 시점(\(j, j'\)) 측정값 간 오차도 독립이라는 것이다.

이것이 핵심적으로 위반되는 가정이다. 다음 섹션에서 자세히 다룬다.

5 독립 가정이 종단 데이터에서 위반되는 이유

5.1 직관: 같은 사람은 같은 사람이다

피험자 \(i\) 가 Week 0에 우울증 점수가 높다면, Week 1에도 상대적으로 높을 가능성이 크다. 왜냐하면 그 사람의 기저 우울증 경향 — 성격, 스트레스 반응성, 치료 감수성 등의 안정적 특성 — 이 모든 시점에 걸쳐 영향을 미치기 때문이다.

이것은 오차(\(\varepsilon\))의 문제가 아니다. 오차는 “예측 불가능한 무작위 변동”이다. 그러나 한 사람이 일관되게 평균보다 높거나 낮은 값을 가지는 것은 예측 가능한 체계적 패턴이다.

단순 회귀에서는 이 체계적 패턴을 표현할 수단이 없다. \(\beta_0\) 는 모든 사람의 공통 절편이기 때문이다. 그 결과, 피험자 \(i\) 가 평균보다 높은 기저 수준을 가지면, 그 사람의 모든 시점 오차가 같은 방향으로 치우치게 된다.

수리적 메커니즘: 오차 간 양의 상관

피험자 \(i\) 의 실제 기저 수준이 \(\beta_0 + d_i\) (여기서 \(d_i > 0\) 이면 평균보다 높음)라고 하자.

단순 회귀 모형 \(y_{ij} = \beta_0 + \beta_1 t_{ij} + \varepsilon_{ij}\) 에서:

\[\varepsilon_{ij} = y_{ij} - \beta_0 - \beta_1 t_{ij} = (d_i + \text{pure noise}_{ij})\]

\(d_i\) 가 모든 시점에 공통으로 들어가므로:

\[\text{Corr}(\varepsilon_{ij}, \varepsilon_{ij'}) = \frac{Cov(d_i + \text{noise}_{ij},\; d_i + \text{noise}_{ij'})}{V(\varepsilon)} = \frac{V(d_i)}{V(\varepsilon)} > 0\]

같은 피험자의 오차들이 양의 상관을 가진다 — 독립 가정 위반.

Reisby 데이터의 상관 행렬을 보면 이것이 실제로 발생한다. Week 0과 Week 5 사이 상관이 0.18~0.22인데, 이것은 오차가 아닌 피험자 간 이질성에서 오는 상관이다.

5.2 독립 가정 위반의 통계적 결과

독립 가정이 위반된 상태에서 OLS로 추정하면 계수 추정값 \(\hat{\beta}\) 자체는 불편(unbiased)이다. 문제는 표준오차다.

같은 피험자의 반복 측정은 독립 관측만큼의 “새로운 정보”를 제공하지 않는다. 피험자 \(i\) 가 6번 측정되었을 때, OLS는 이를 6개의 독립 관측으로 계산한다. 그러나 실제로 이 6개는 서로 상관된 관측이므로, 정보량은 6개의 독립 관측보다 적다.

OLS의 착각	현실
\(N \times n\) 개의 독립 관측이 있다	실제 독립 정보는 훨씬 적다
표준오차 = \(\sigma / \sqrt{N \times n}\)	실제 표준오차 \(\gg\) OLS 추정값
신뢰구간이 좁다	신뢰구간이 실제보다 좁음 — 과신

결과: 표준오차를 과소 추정 → 검정 통계량이 과대 → 1종 오류(유의하지 않은 효과를 유의하다고 판정) 증가.

얼마나 심각한가

ICC(급내상관계수)가 0.5이고 측정 횟수가 5회라면, 실제 유효 표본 크기는 대략:

\[n_{\text{eff}} = \frac{n}{1 + (n-1) \cdot \text{ICC}} = \frac{5}{1 + 4 \times 0.5} = \frac{5}{3} \approx 1.67\]

OLS는 5개의 관측을 5개 독립 관측으로 계산하지만, 실제 정보량은 2개도 안 된다. ICC가 높을수록 이 왜곡은 커진다.

6 두 번째 문제: 이질성 무시 (Heterogeneity)

독립 가정 위반과는 별개로, 단순 회귀의 \(\beta_1\) (모든 사람에게 동일한 기울기)도 문제다.

6.1 현실의 개인 간 이질성

임상 연구에서 “평균적으로 치료 효과가 있다”는 결론은 가치 있다. 그러나 그 평균 뒤에는 큰 이질성이 숨어 있다:

일부 환자: 치료에 강하게 반응하여 우울증 점수가 급격히 하락
다른 환자: 치료에 반응하지 않거나 오히려 악화
나머지 환자: 중간 수준의 반응

단순 회귀에서 \(\beta_1\) 은 이 모든 것의 평균이다. 개인별 변화 속도 \(b_{1i}\) 는 추정할 수 없다.

Reisby 데이터의 표준편차를 보면 이 이질성이 시간에 따라 증가한다. Week 0에서 sd = 4.53이었다가 Week 5에서 sd = 7.22로 커진다. 이것은 시간이 지날수록 피험자 간 반응 차이가 벌어진다는 신호다 — 모든 사람이 같은 기울기를 가진다는 가정과 정면으로 충돌한다.

직관: 스파게티 플롯이 보여주는 것

Reisby 데이터의 스파게티 플롯(각 피험자의 시계열을 개별 선으로 그린 그래프)을 상상해 보자.

공통 가정: 모든 선이 평행(동일 기울기), 시작점만 다름
실제 관측: 어떤 선은 급격히 하강, 어떤 선은 수평, 어떤 선은 일부 구간에서 상승

단순 회귀는 이 모든 선을 하나의 평균선으로 요약한다. MRM은 각 선을 개별적으로 추정하면서 공통 패턴도 포착한다.

7 모형의 한계 요약

단순 선형회귀 \(y_{ij} = \beta_0 + \beta_1 t_{ij} + \varepsilon_{ij}\) 를 종단 데이터에 적용할 때 발생하는 두 가지 근본적 문제:

문제	원인	결과
독립 가정 위반	같은 피험자의 반복 측정은 상관됨	표준오차 과소 추정, 1종 오류 증가
이질성 무시	\(\beta_0\), \(\beta_1\) 에 \(i\) 첨자 없음	개인별 변화 추세 추정 불가, 집단 평균만 추정

이 두 문제는 같은 원인에서 나온다. 피험자 개인의 특성(\(d_i\), 또는 \(v_{0i}\))이 모형에 포함되지 않기 때문이다. 이 개인 특성을 모형에 명시적으로 포함시키면:

같은 피험자의 관측값 간 상관이 설명됨 → 독립 가정 문제 해소
개인별 절편/기울기 추정 가능 → 이질성 포착

이것이 정확히 혼합효과 회귀모형이 하는 일이다.

8 Reisby 데이터: 단순 회귀의 한계를 직접 확인

Reisby et al. (1977)의 정신과 임상시험 데이터를 사용한다. 66명의 우울증 입원 환자를 대상으로 이미프라민(imipramine) 225mg/day를 4주간 투여하고, Hamilton 우울증 평가 척도(HDRS)를 6시점(Week 0~5)에서 측정했다. HDRS 점수가 높을수록 우울증이 심하다.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 데이터 로드 (Hedeker 홈페이지 reisby.dat)
# 실제 실행 시: pd.read_csv("reisby.dat", ...)
# 여기서는 구조 시연을 위해 요약 통계만 사용

# 교재 Table 4.1: 시점별 HDRS 통계
weeks = [0, 1, 2, 3, 4, 5]
means = [23.44, 21.84, 18.31, 16.42, 13.62, 11.95]
sds   = [4.53,  4.70,  5.49,  6.42,  6.97,  7.22]
ns    = [61,    63,    65,    65,    63,    58]

# 시점별 분산 변화 — 등분산 가정 위반 시각화
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

ax1.errorbar(weeks, means, yerr=sds, fmt='o-', capsize=5)
ax1.set_xlabel("Week")
ax1.set_ylabel("HDRS Mean (±SD)")
ax1.set_title("HDRS Mean and SD over Time")
ax1.grid(alpha=0.3)

ax2.plot(weeks, sds, 'rs-')
ax2.set_xlabel("Week")
ax2.set_ylabel("Standard Deviation")
ax2.set_title("SD increases over time\n(heteroscedasticity)")
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

관찰

평균 감소: Week 0 (23.44) → Week 5 (11.95). 치료 효과가 전반적으로 존재함
표준편차 증가: 4.53 → 7.22. 시간이 지날수록 피험자 간 이질성이 커짐
완전 데이터: 66명 중 46명만 (70%). MANOVA는 46명만 사용하지만 MRM은 66명 모두 사용

import statsmodels.formula.api as smf
import statsmodels.api as sm

# 실제 long-format 데이터가 있다면:
# df = pd.read_csv("reisby_long.csv")

# Step 1: 단순 OLS (독립 가정 위반 상태)
# 모든 관측을 독립으로 취급
model_ols = smf.ols("hdrs ~ week", data=df).fit()
print(model_ols.summary())

# OLS 계수: 모든 피험자에 공통인 β_0, β_1
# 문제: 같은 피험자의 반복 측정을 독립으로 취급 → SE 과소 추정

# Step 2: 실제 표준오차와 OLS 표준오차 비교
# (군집 강건 표준오차 — 피험자 군집을 고려)
model_cluster = smf.ols("hdrs ~ week", data=df).fit(
    cov_type='cluster', cov_kwds={'groups': df['id']}
)
print("\n=== OLS SE vs Cluster-Robust SE ===")
print(f"OLS SE(β_1):     {model_ols.bse['week']:.4f}")
print(f"Cluster SE(β_1): {model_cluster.bse['week']:.4f}")
# 기대: Cluster SE > OLS SE (OLS가 SE를 과소 추정함을 확인)

군집 강건 표준오차와 OLS 표준오차를 비교하면 OLS가 표준오차를 얼마나 과소 추정하는지 직접 확인할 수 있다.

# R: 동일한 비교
library(nlme)
library(dplyr)

# OLS (독립 가정)
fit_ols <- lm(hdrs ~ week, data = reisby_long)

# 군집 강건 SE
library(sandwich)
library(lmtest)
coeftest(fit_ols, vcov = vcovCL(fit_ols, cluster = ~id))

# 결과 해석:
# OLS SE가 실제보다 작으면 → 유의성이 과장됨

9 관측 상관의 패턴 확인

Reisby 데이터의 상관 행렬(Table 4.2)을 보면 복합대칭 가정이 성립하는지 확인할 수 있다.

# 교재 Table 4.2 재현: 시점 간 상관 (pairwise)
# (아래는 교재에서 직접 인용한 값)
corr_matrix = np.array([
    [1.00, 0.49, 0.41, 0.33, 0.23, 0.18],
    [0.49, 1.00, 0.49, 0.41, 0.31, 0.22],
    [0.42, 0.49, 1.00, 0.74, 0.67, 0.46],
    [0.44, 0.51, 0.73, 1.00, 0.82, 0.57],
    [0.30, 0.35, 0.68, 0.78, 1.00, 0.65],
    [0.22, 0.23, 0.53, 0.62, 0.72, 1.00],
])

import seaborn as sns
plt.figure(figsize=(7, 6))
sns.heatmap(corr_matrix, annot=True, fmt=".2f", cmap="Blues",
            xticklabels=[f"Wk{i}" for i in range(6)],
            yticklabels=[f"Wk{i}" for i in range(6)])
plt.title("HDRS Correlations (Pairwise)\n복합대칭이면 모든 비대각 원소가 동일해야 함")
plt.tight_layout()
plt.show()

상관 행렬이 말하는 것

복합대칭(compound symmetry) 가정: 모든 비대각 상관이 동일해야 한다.

실제 값을 보면: - Week 0 ↔︎ Week 1: 0.49 - Week 0 ↔︎ Week 5: 0.18 - Week 4 ↔︎ Week 5: 0.65

멀리 떨어진 시점일수록 상관이 낮아지는 대각 우세(banded) 패턴이 나타난다. 이것은 복합대칭 가정의 위반이며, 단순 랜덤 절편 모형으로도 완전히 설명되지 않는다.

MRM에서 랜덤 절편 + 기울기 모형은 이 증가하는 이질성 패턴을 더 잘 포착한다.

10 MRM으로의 자연스러운 이행

두 문제를 동시에 해결하는 방법은 놀랍도록 자연스럽다. 모형에 피험자별 개인 효과 \(v_{0i}\) 를 추가하면 된다:

\[y_{ij} = \beta_0 + \beta_1 t_{ij} + \underbrace{v_{0i}}_{\text{추가!}} + \varepsilon_{ij}\]

이 하나의 추가가 어떻게 두 문제를 동시에 해결하는지는 다음 포스트(랜덤 절편 MRM)에서 상세히 다룬다.

미리 결론을 말하면:

문제	해결 방식
독립 가정 위반	\(v_{0i}\) 가 피험자 내 상관을 흡수. \(\varepsilon_{ij}\) 는 \(v_{0i}\) 에 조건부로 독립.
이질성 무시	각 피험자의 고유 효과 \(v_{0i}\) 가 추정됨. 기울기도 \(v_{1i}\) 를 추가하면 피험자별로 달라짐.

MRM은 단순 선형회귀 + 개인 랜덤 효과 = 증강된 선형 회귀 모형이다.

11 관련 주제

선행 지식

후속 주제

혼합효과 회귀모형 — 연속 반응변수 개요 (랜덤 절편, 랜덤 기울기, 행렬 표현)

관련 개념