Linear Mixed Model (3): 추정과 모델 선택

ML vs REML, LRT, AIC/BIC, 수렴 문제

LMM의 파라미터 추정 방법인 ML과 REML의 차이, 모델 간 비교를 위한 Likelihood Ratio Test(LRT), 정보 기준(AIC/BIC), 그리고 실무에서 자주 만나는 수렴 문제와 해결법을 다룬다.

Statistics
Longitudinal Data
저자

Kwangmin Kim

공개

2026년 03월 07일

1 Linear Mixed Model (3): 추정과 모델 선택

1.1 ML vs REML

LMM의 파라미터 \((\beta, \sigma^2_u, \sigma^2_e)\)를 추정하는 두 가지 방법이 있다.

1.1.1 Maximum Likelihood (ML)

고정 효과와 분산 성분을 동시에 최대화한다.

\[\hat{\theta}_{ML} = \arg\max_\theta \log L(\beta, \sigma^2_u, \sigma^2_e | Y)\]

문제: 고정 효과 \(\beta\)의 추정 불확실성을 무시하므로 분산 성분을 과소추정한다.

직관적 예시:
- 데이터 10개, 고정 효과 파라미터 5개 추정
- 실제로 분산 추정에 활용 가능한 자유도 = 10 - 5 = 5
- ML은 이 손실을 반영하지 않음 → 분산 과소추정

1.1.2 Restricted Maximum Likelihood (REML)

고정 효과를 먼저 제거한 잔차에서 분산 성분을 추정한다.

\[\hat{\theta}_{REML} = \arg\max_\theta \log L_{REML}(\sigma^2_u, \sigma^2_e | Y, \hat{\beta})\]

  • 고정 효과 추정의 불확실성을 자유도로 반영
  • 분산 성분 추정에 더 정확 (unbiased)
  • 기본값: lme4::lmer(), statsmodels.mixedlm() 모두 기본이 REML

1.1.3 언제 ML을 써야 하는가

목적 방법 이유
분산 성분 추정 REML 불편 추정
랜덤 효과 구조 비교 (동일 fixed effects) REML REML로 랜덤 구조 LRT 가능
고정 효과 구조 비교 (다른 fixed effects) ML REML은 fixed effects 다르면 비교 불가
최종 모델 해석 REML 분산 추정이 더 정확

핵심 규칙:

Fixed Effects 구조를 비교할 때 → ML로 재적합 후 LRT
Random Effects 구조를 비교할 때 → REML 사용
최종 해석 → REML

1.1.4 수치 예시

import statsmodels.formula.api as smf

# REML (기본값)
model_reml = smf.mixedlm(
    "satisfaction ~ personalized + week",
    data=df, groups=df["user_id"]
).fit(reml=True)   # 기본값

# ML
model_ml = smf.mixedlm(
    "satisfaction ~ personalized + week",
    data=df, groups=df["user_id"]
).fit(reml=False)

print(f"REML 분산 추정: σ²_u = {model_reml.cov_re.iloc[0,0]:.3f}")
print(f"ML 분산 추정:   σ²_u = {model_ml.cov_re.iloc[0,0]:.3f}")
# REML이 항상 같거나 더 큰 분산을 추정 (과소추정 교정)
출력 예시:
REML 분산 추정: σ²_u = 0.421
ML 분산 추정:   σ²_u = 0.398   ← 약간 과소추정

1.2 모델 선택: LRT (Likelihood Ratio Test)

두 중첩된(nested) 모델을 비교할 때 사용한다.

1.2.1 수식

\[\Lambda = -2(\log L_0 - \log L_1) \sim \chi^2_{df_1 - df_0}\]

  • \(L_0\): 단순 모델의 likelihood
  • \(L_1\): 복잡 모델의 likelihood
  • 자유도 = 추가된 파라미터 수

1.2.2 Fixed Effects 비교 (ML 사용)

# ML로 두 모델 적합
m1 = smf.mixedlm("satisfaction ~ week", data=df,
                  groups=df["user_id"]).fit(reml=False)

m2 = smf.mixedlm("satisfaction ~ week + personalized", data=df,
                  groups=df["user_id"]).fit(reml=False)

# LRT
from scipy.stats import chi2

lr_stat = 2 * (m2.llf - m1.llf)
df_diff = m2.df_modelwc - m1.df_modelwc   # 파라미터 차이
p_value = 1 - chi2.cdf(lr_stat, df=df_diff)

print(f"LRT: χ²({df_diff}) = {lr_stat:.2f}, p = {p_value:.4f}")
# LRT: χ²(1) = 48.23, p < 0.001
# → personalized 추가가 모델을 유의하게 개선함

1.2.3 Random Effects 비교 (REML 사용)

library(lme4)
library(lmtest)

# REML로 두 모델 적합
m_ri <- lmer(satisfaction ~ personalized + week + (1 | user_id),
             data = df, REML = TRUE)

m_rs <- lmer(satisfaction ~ personalized + week + (1 + week | user_id),
             data = df, REML = TRUE)

# LRT (REML 기반)
anova(m_ri, m_rs, refit = FALSE)   # refit=FALSE: REML 유지
결과:
        Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
m_ri     5 12234 12267  -6112    12224
m_rs     7 12198 12244  -6092    12184   40.2      2  1.9e-09 ***

→ Random Slope 모델이 유의하게 나음 (χ²(2) = 40.2, p < 0.001)

주의: LRT로 Random Effects 비교 시, \(H_0\)가 경계에 있어 (\(\sigma^2 \geq 0\)) p-value가 보수적이다. 실제 유의성은 더 강할 수 있다.

1.3 정보 기준: AIC와 BIC

비중첩 모델이나 여러 모델을 동시에 비교할 때 사용한다.

1.3.1 수식

\[\text{AIC} = -2\log L + 2k\] \[\text{BIC} = -2\log L + k \log(n)\]

  • \(k\): 추정된 파라미터 수
  • \(n\): 관측치 수
  • 낮을수록 좋음

1.3.2 AIC vs BIC 비교

AIC BIC
패널티 \(2k\) \(k \log n\)
큰 n에서 파라미터 패널티 고정 파라미터 패널티 증가
선호 경향 복잡한 모델 (더 flexible) 단순한 모델 (더 conservative)
목적 예측 성능 진짜 모델 찾기
추천 예측이 목적 설명이 목적

1.3.3 실무 예시

import pandas as pd

models = {
    "Null (intercept only)": smf.mixedlm("satisfaction ~ 1", data=df, groups=df["user_id"]).fit(reml=False),
    "Main effects": smf.mixedlm("satisfaction ~ personalized + week", data=df, groups=df["user_id"]).fit(reml=False),
    "With interaction": smf.mixedlm("satisfaction ~ personalized * week", data=df, groups=df["user_id"]).fit(reml=False),
    "With segment": smf.mixedlm("satisfaction ~ personalized + week + segment", data=df, groups=df["user_id"]).fit(reml=False),
}

comparison = pd.DataFrame({
    name: {"AIC": m.aic, "BIC": m.bic, "logLik": m.llf, "params": m.df_modelwc}
    for name, m in models.items()
}).T.sort_values("AIC")

print(comparison)
                      AIC      BIC    logLik  params
With interaction   12185.2  12231.4  -6085.6     7
With segment       12190.1  12243.3  -6087.1     8
Main effects       12198.3  12237.5  -6092.2     6
Null               12890.7  12916.8  -6441.4     3

→ AIC 기준: With interaction이 최선
→ Main effects 대비 Δ AIC = 13.1 (> 10: 강한 지지)

1.3.4 ΔAIC 해석 기준 (Burnham & Anderson)

\(\Delta\)AIC 해석
0 ~ 2 두 모델 비슷
2 ~ 4 약간의 지지 차이
4 ~ 10 상당한 지지 차이
> 10 강한 지지 차이

1.4 모델 선택 전략: Bottom-up vs Top-down

1.4.1 Bottom-up (권장)

단순 모델에서 시작해 복잡도를 점진적으로 높인다.

Step 1: Null model (랜덤 절편만)
        → ICC 확인

Step 2: Fixed Effects 추가 (ML로 LRT)
        → 이론적으로 중요한 변수부터

Step 3: Random Slope 필요 여부 확인 (REML로 LRT)
        → 이론적 근거 있을 때만 추가

Step 4: 최종 모델 REML로 재적합
        → 보고용 파라미터 추출

1.4.2 Top-down

Full model에서 시작해 불필요한 항을 제거한다.

  • 가설이 명확할 때 사용
  • Overfitting 위험 있음

1.4.3 실무 권장 코드 (Python)

import statsmodels.formula.api as smf
from scipy.stats import chi2

def lrt(model_simple, model_complex):
    """Likelihood Ratio Test"""
    stat = 2 * (model_complex.llf - model_simple.llf)
    df = model_complex.df_modelwc - model_simple.df_modelwc
    p = 1 - chi2.cdf(stat, df)
    return {"LRT_stat": round(stat, 3), "df": df, "p_value": round(p, 4)}

# Step 1: Null
m0 = smf.mixedlm("satisfaction ~ 1", data=df, groups=df["user_id"]).fit(reml=False)

# Step 2: Fixed effects
m1 = smf.mixedlm("satisfaction ~ week", data=df, groups=df["user_id"]).fit(reml=False)
m2 = smf.mixedlm("satisfaction ~ week + personalized", data=df, groups=df["user_id"]).fit(reml=False)
m3 = smf.mixedlm("satisfaction ~ week + personalized + segment", data=df, groups=df["user_id"]).fit(reml=False)

print("week 추가:", lrt(m0, m1))
print("personalized 추가:", lrt(m1, m2))
print("segment 추가:", lrt(m2, m3))
출력:
week 추가:          {'LRT_stat': 12.4, 'df': 1, 'p_value': 0.0004}
personalized 추가:  {'LRT_stat': 48.2, 'df': 1, 'p_value': 0.0000}  ← 강한 효과
segment 추가:       {'LRT_stat': 31.5, 'df': 2, 'p_value': 0.0000}  ← 세그먼트도 유의

1.5 수렴 문제 (Convergence Issues)

Random Slope 모델에서 자주 발생하는 실무 문제.

1.5.1 수렴 실패 증상

ConvergenceWarning: Singular matrix
또는
Model failed to converge with max|grad| = 0.003 (tol = 0.002)

1.5.2 원인과 해결책

원인 증상 해결책
분산 추정값이 0에 근접 isSingular() = TRUE 해당 랜덤 효과 제거
파라미터가 너무 많음 느린 수렴 모델 단순화
스케일 차이 큼 수렴 불안정 변수 표준화
그룹 수 부족 추정 불안정 더 많은 데이터 필요

1.5.3 해결 방법

방법 1: 변수 표준화

from sklearn.preprocessing import StandardScaler

df['week_std'] = (df['week'] - df['week'].mean()) / df['week'].std()
df['satisfaction_std'] = (df['satisfaction'] - df['satisfaction'].mean()) / df['satisfaction'].std()

model = smf.mixedlm(
    "satisfaction_std ~ personalized + week_std",
    data=df, groups=df["user_id"], re_formula="~week_std"
).fit()

방법 2: 공분산 제거 (Diagonal)

절편-기울기 공분산을 0으로 고정:

# R: || 사용
model_diag <- lmer(
  satisfaction ~ personalized + week + (1 | user_id) + (0 + week | user_id),
  data = df
)

방법 3: 랜덤 기울기 제거

# 수렴 안 되면 Random Intercept Only로 후퇴
model_ri <- lmer(
  satisfaction ~ personalized + week + (1 | user_id),
  data = df
)

방법 4: 옵티마이저 변경

model <- lmer(
  satisfaction ~ personalized + week + (1 + week | user_id),
  data = df,
  control = lmerControl(optimizer = "bobyqa",
                        optCtrl = list(maxfun = 2e5))
)
# Python: 옵티마이저 변경
model = smf.mixedlm(...).fit(method=['lbfgs', 'nm', 'cg'])

1.6 샘플 크기 고려사항

Mixed Model에서 “샘플 크기”는 단순하지 않다.

1.6.1 두 가지 샘플 크기

관측치 수 (\(N\)) 그룹 수 (\(J\))
Fixed Effects 추정 \(N\) 중요 보조적
Random Effects 추정 보조적 \(J\)가 중요
일반적 권장 충분히 많을수록 최소 10~30개

1.6.2 구체적 권장 기준

목적 그룹 수 그룹당 관측치
랜덤 절편 추정 ≥ 10 ≥ 5
랜덤 기울기 추정 ≥ 20 ≥ 10
복잡한 공분산 ≥ 30 ≥ 15
안정적 추정 ≥ 50 ≥ 20

실무 예시:

AI Agent 개인화 실험:
- 사용자 수 (그룹 수): 500명
- 사용자당 세션 수: 평균 8회
- 총 관측치: 4,000건

→ Random Intercept: ✅ (500명 > 10)
→ Random Slope: ✅ (500명 > 20, 세션 8회 > 10)
→ 안정적 추정 가능

1.7 고정 효과 추정값 해석

1.7.1 t-value vs z-value

소프트웨어 통계량 p-value 제공
lme4 (R) t-value ❌ (자유도 불명확)
nlme (R) t-value ✅ (근사 자유도)
statsmodels (Python) z-value ✅

lme4에서 p-value 구하기:

library(lmerTest)   # lme4 위에 덮어쓰기

model <- lmer(satisfaction ~ personalized + week + (1 | user_id), data = df)
summary(model)   # 이제 p-value 포함

Satterthwaite 자유도 근사 (lmerTest 기본):

  • 고정 효과의 자유도를 Satterthwaite 방법으로 근사
  • 샘플이 충분히 크면 z-test와 결과 거의 같음

1.7.2 95% 신뢰구간

# Wald CI (근사)
confint(model, method = "Wald")

# Profile CI (정확하지만 느림)
confint(model, method = "profile")

# Bootstrap CI (가장 정확, 매우 느림)
confint(model, method = "boot", nsim = 1000)
# Python: statsmodels는 Wald CI 제공
model.conf_int()

1.8 전체 분석 워크플로우 요약

# 완전한 분석 파이프라인
import statsmodels.formula.api as smf
import pandas as pd

# 1. ICC 확인 (REML)
null = smf.mixedlm("Y ~ 1", data=df, groups=df["id"]).fit(reml=True)
icc = null.cov_re.iloc[0,0] / (null.cov_re.iloc[0,0] + null.scale)
print(f"ICC = {icc:.3f}")  # 0.05 이상이면 Mixed Model 필요

# 2. Fixed Effects 구조 결정 (ML)
m1 = smf.mixedlm("Y ~ X", data=df, groups=df["id"]).fit(reml=False)
m2 = smf.mixedlm("Y ~ X + Z", data=df, groups=df["id"]).fit(reml=False)
# LRT로 Z 추가 여부 결정

# 3. Random Effects 구조 결정 (REML)
m_ri = smf.mixedlm("Y ~ X + Z", data=df, groups=df["id"]).fit(reml=True)
m_rs = smf.mixedlm("Y ~ X + Z", data=df, groups=df["id"],
                    re_formula="~X").fit(reml=True)
# LRT로 Random Slope 필요 여부 결정

# 4. 최종 모델 (REML)
final = smf.mixedlm("Y ~ X + Z", data=df, groups=df["id"],
                     re_formula="~X").fit(reml=True)

# 5. 결과 보고
print(final.summary())
print("\n분산 성분:")
print(f"  σ²_u0 = {final.cov_re.iloc[0,0]:.3f} (랜덤 절편)")
print(f"  σ²_e  = {final.scale:.3f} (잔차)")

다음: [04-mixed-model-practice.qmd] — AI Agent 개인화 실험 end-to-end 실무 예시

Subscribe

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