ML 추정과 REML — 분산 성분 추정의 두 방법

Maxwell Ch.15.2 ML vs ANOVA Approach · REML

Mixed-effects model 의 추정 방법 — Maximum Likelihood (ML) 와 Restricted Maximum Likelihood (REML) 의 차이, 적용 상황, ANOVA 기반 추정과의 관계, LRT (Likelihood Ratio Test) 의 수행 방법, REML 의 작은 표본 우위를 단계적으로 정리한다.

Experimentation
DOE
저자

Kwangmin Kim

공개

2026년 05월 08일

1 정의

정의: ML vs REML
방법 Likelihood 분산 성분 추정 Fixed effect 추정
ML \(\log L(\boldsymbol{\beta}, \boldsymbol{\theta} | \mathbf{Y})\) 약간 편향 (분산 underestimate) 비편향
REML Fixed effect 의 자유도 손실을 보정 비편향 (별도 계산)
ANOVA-based 모형 비교의 EMS 균등 데이터에서 정확. 비균등에서 편향 비편향

2 ML 추정

전체 likelihood 를 fixed effect 와 variance component 에 대해 최대화.

\[ L(\boldsymbol{\beta}, \boldsymbol{\theta}) = \frac{1}{(2\pi)^{n/2} |\mathbf{V}|^{1/2}} \exp\left(-\frac{1}{2}(\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})\right) \]

여기서 \(\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^T + \mathbf{R}\) 은 variance components 의 함수.

문제: \(\boldsymbol{\beta}\) 추정의 자유도 손실이 분산 추정에 반영되지 않아 분산 성분이 약간 과소 추정.

2.1 ML 의 편향

\(n\) 명 데이터에서 fixed effect 모수 수 \(p\) 가 있으면, 잔차 자유도 \(n - p\).

ML 추정량: \(\hat\sigma^2_{\text{ML}} = SS_E / n\) (편향). 비편향 추정량: \(\hat\sigma^2 = SS_E / (n - p)\).

ML 이 약 \((n-p)/n\) 배 underestimate. 작은 표본에서 큰 편향.

3 REML

분산 성분만의 likelihood (fixed effect 의 자유도 손실 보정):

\[ \log L_{\text{REML}}(\boldsymbol{\theta}) = -\frac{1}{2}\log|\mathbf{V}| - \frac{1}{2}\log|\mathbf{X}^T \mathbf{V}^{-1} \mathbf{X}| - \frac{1}{2}(\mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}})^T \mathbf{V}^{-1} (\mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}}) \]

추가 항 \(-\frac{1}{2}\log|\mathbf{X}^T \mathbf{V}^{-1} \mathbf{X}|\) 가 자유도 손실을 보정.

분산 성분이 비편향. 작은 표본에서 우수.

직관: REML 의 보정

\(n\) 명 데이터에서 잔차 분산 \(\sigma^2\) 의 ML 추정량은 \(SS_E / n\), 비편향 추정량은 \(SS_E / (n - p)\) (\(p\) = fixed effect 모수 수). REML 은 이 자유도 보정을 mixed model 의 모든 분산 성분에 일반화한다.

→ ANOVA 의 \(MS_E = SS_E / df_E\) 와 동일한 정신. ANOVA 가 EMS 식으로 자유도 보정 한 것을 likelihood 일반화.

4 사용 시점

상황 권장
분산 성분 추정 REML
Fixed effect 의 LRT 비교 ML (REML 로는 LRT 불가능)
공분산 구조 비교 (LRT) REML
최종 결과 보고 REML
Wald 검정 (단일 fixed effect) 둘 다 OK

5 ANOVA-Based Estimation

균등 데이터에서 ANOVA 의 EMS 식으로 분산 성분 추정 (G-MAX10-1 참조).

방법 균등 비균등 결측 음수 추정
ANOVA 정확 편향 어려움 가능
ML 정확 정확 처리 가능 가능 (rare)
REML 정확 정확 처리 가능 가능 (rare)

REML 이 일반 권장. ANOVA 는 educational interest.

6 가설 데이터 — 비교

\(n = 20\), \(T = 5\), random intercept 모형. True \(\sigma^2_0 = 36\), \(\sigma^2 = 9\).

추정 결과 (가상):

방법 \(\hat\sigma^2_0\) \(\hat\sigma^2\)
True 36 9
ANOVA 35 9.2
ML 33 9
REML 35.5 9

ML 이 약간 underestimate, REML 이 unbiased.

7 LRT (Likelihood Ratio Test)

정의: LRT

두 nested 모형 비교:

\[ LRT = -2(\log L_{\text{reduced}} - \log L_{\text{full}}) \]

자유도 = 모수 수 차이. \(\chi^2\) 분포.

7.1 Fixed Effect 비교

ML 추정 필수. REML 의 likelihood 는 fixed effect 의 함수가 아니므로 LRT 무의미.

# Python statsmodels (ML refit)
md1 = mixedlm("Y ~ time", data, groups="subject").fit(reml=False)
md2 = mixedlm("Y ~ time + group", data, groups="subject").fit(reml=False)
LRT = 2 * (md2.llf - md1.llf)

7.2 Variance Component 비교

REML OK. 단 \(H_0: \sigma^2 = 0\) (boundary) 검정은 표준 \(\chi^2\) 분포가 아닌 mixture distribution 사용.

# Variance component LRT
md1 = mixedlm("Y ~ time", data, groups="subject").fit()  # random intercept
md2 = mixedlm("Y ~ time", data, groups="subject", re_formula="~time").fit()
                                                                     # + slope
LRT = 2 * (md2.llf - md1.llf)
# p-value: 0.5 * P(chi^2_1 > LRT) + 0.5 * P(chi^2_2 > LRT)
# (boundary mixture)
함정: Boundary Test 의 mixture distribution

\(\sigma^2 = 0\) 검정은 모수 공간의 boundary 에서 — chi-square 가 부적합.

올바른 분포: \(0.5 \chi^2_0 + 0.5 \chi^2_1\) (자유도 1 추가 시).

표준 LRT (just \(\chi^2_1\)) 사용 시 \(p\) 값 conservative (실제보다 큼). variance component 가 0 보다 살짝 큰 데이터를 detect 하기 어려움.

R RLRsim 패키지가 시뮬레이션 기반 정확한 검정.

8 가설 데이터 — LRT 적용

random intercept vs random intercept + slope 비교.

8.1 ML 추정

md_int (random intercept):     log L = -120, AIC = 248
md_slope (intercept + slope):  log L = -115, AIC = 240

LRT = \(2 \times (-115 - (-120)) = 10\). 자유도 2 (slope variance + intercept-slope cov).

\(p\)-value (mixture): 약 0.005. 유의 → random slope 가 의미.

8.2 Wald 검정 vs LRT

Fixed effect 검정에서: - Wald: \(\hat\beta / SE(\hat\beta) \to z\). 빠름. - LRT: 두 모형 likelihood 비교. 정확.

작은 표본에서는 LRT 가 더 정확. 큰 표본에서는 두 결과 비슷.

9 REML 의 limitation

함정: REML 의 LRT 한계

REML 의 likelihood 는 fixed effect 의 함수가 아님. 따라서:

  • Fixed effect 가 다른 두 모형의 REML LRT → 무의미.
  • Variance component 가 다른 두 모형의 REML LRT → OK.

해결: - Fixed effect 비교 → ML 로 refit 후 LRT. - 최종 보고 → REML 결과.

R lme4: update(model, REML=FALSE) 로 ML refit. Python: mixedlm(..., reml=False).

10 Python 코드

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import mixedlm

np.random.seed(2026)
n = 30
T = 5
sigma_0 = 6
sigma_e = 3

records = []
for subj in range(n):
    u0 = np.random.normal(0, sigma_0)
    for t in range(T):
        y = 50 + 3 * t + u0 + np.random.normal(0, sigma_e)
        records.append({"subject": subj, "time": t, "Y": y})
data = pd.DataFrame(records)

# REML (default)
md_reml = mixedlm("Y ~ time", data=data, groups=data["subject"]).fit()
print("=== REML ===")
print(f"sigma_0² = {md_reml.cov_re.iloc[0,0]:.3f} (true: {sigma_0**2})")
print(f"sigma²   = {md_reml.scale:.3f} (true: {sigma_e**2})")

# ML
md_ml = mixedlm("Y ~ time", data=data, groups=data["subject"]).fit(reml=False)
print("\n=== ML ===")
print(f"sigma_0² = {md_ml.cov_re.iloc[0,0]:.3f}")
print(f"sigma²   = {md_ml.scale:.3f}")
print(f"log-likelihood = {md_ml.llf:.3f}")

# LRT for fixed effect: time vs no time
md_no_time = mixedlm("Y ~ 1", data=data, groups=data["subject"]).fit(reml=False)
LRT = 2 * (md_ml.llf - md_no_time.llf)
from scipy import stats as sps
p_lrt = 1 - sps.chi2.cdf(LRT, 1)
print(f"\nLRT for time effect: chi² = {LRT:.2f}, df=1, p = {p_lrt:.4f}")

# Variance component LRT (random slope)
md_slope_ml = mixedlm("Y ~ time", data=data, groups=data["subject"],
                       re_formula="~time").fit(reml=False)
LRT_slope = 2 * (md_slope_ml.llf - md_ml.llf)
# Mixture chi-square (0.5 chi^2_1 + 0.5 chi^2_2)
p_slope_naive = 1 - sps.chi2.cdf(LRT_slope, 2)  # naive
p_slope_mix = 0.5 * (1 - sps.chi2.cdf(LRT_slope, 1)) + 0.5 * (1 - sps.chi2.cdf(LRT_slope, 2))
print(f"\nLRT for random slope: chi² = {LRT_slope:.2f}")
print(f"  naive p (df=2): {p_slope_naive:.4f}")
print(f"  mixture p:      {p_slope_mix:.4f}")

11 가정과 한계

  • ML 의 작은 표본 편향.
  • REML 의 LRT 한계 (fixed effect 비교 불가).
  • Variance component 의 boundary 검정: 표준 \(\chi^2\) 부적합.
  • 수렴 문제: 작은 표본·복잡한 random structure.

12 모형 비교 절차 — 종합

Step 1: Random structure 선택
  - REML 추정.
  - LRT 또는 AIC/BIC 로 random component 비교.
  - $\sigma^2 = 0$ boundary 검정 주의.

Step 2: Fixed effect 선택
  - ML 추정 (LRT 가능하도록).
  - LRT 또는 AIC.

Step 3: 최종 모형
  - REML 로 refit.
  - Wald 검정 또는 신뢰구간 보고.

13 ML 매핑

ML 의 ML vs REML

ML 의 cross-validation: - 각 fold 의 정확도 → response. - 모델 (fixed) × fold (random) → mixed model.

REML 추정으로 fold 의 random variance 를 정확하게 추정.

작은 fold 수 (예: 5-fold) 에서 REML 이 ML 보다 정확. 모델 차이 검정 시 ML refit + LRT.

14 본 시리즈

G-MAX15-0  Multilevel 개관
G-MAX15-1  Y = Xβ + Zu + ε
G-MAX15-2  ML vs REML  ← 현재 글
G-MAX15-3  Growth Curve + Covariance Structures
G-MAX15-4  Model Comparison + Time-varying Covariates

15 관련 주제

선행 지식

후속 주제

다른 카테고리 연결

16 더 읽을 거리

  • Patterson, H. D., Thompson, R. (1971). “Recovery of inter-block information when block sizes are unequal.” Biometrika 58(3): 545-554 — REML 원조.
  • Harville, D. A. (1977). “Maximum likelihood approaches to variance component estimation and to related problems.” JASA 72(358): 320-338.
  • Searle, S. R., Casella, G., McCulloch, C. E. (2006). “Variance Components.” Wiley.
  • Verbeke, G., Molenberghs, G. (2003). “The use of score tests for inference on variance components.” Biometrics 59(2): 254-262 — boundary test.
  • Pinheiro, J. C., Bates, D. M. (2000). “Mixed-Effects Models in S and S-PLUS.” Springer.

Subscribe

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