개요

linear mixed model에서 intercept에만 random effects를 주었을 때 추정치들의 standard error가 같아지는 경우에 대해 알아봅시다.

1. 반복측정자료

어떤 임상실험에서 연구자가 특정 기간을 주기로 당뇨병 환자들의 인슐린 농도를 측정한다고 해봅시다. 이러한 경우 환자 각각은 서로 다른 사람이므로, 서로 다른 환자에게서 관찰된 데이터는 서로 독립이라고 봐도 무방합니다.

그러나 동일한 환자에게서 반복적으로 측정된 값들은 서로 독립이라고 보기 어렵습니다. 동일한 환자에게서 측정정되었기 때문에, 각 시점의 데이터들 간에는 특정한 상관관계가 존재할 것이라는 생각을 해볼 수 있겠습니다. 예컨대, 당뇨임에도 치료 효과가 좋은 환자는 시간이 지날수록 인슐린 농도가 점점 높아질 수 있고, 치료에 소극적인 환자의 경우에는 반대의 경향성을 보일 것입니다.

일반적으로 회귀분석에서는 모든 관측지(y, 종속변수)들이 서로 독립이라고 가정합니다. 그러나 반복측정자료의 경우 위와 같은 이유 때문에 회귀분석을 통해 분석하기에 어려움이 존재합니다. 이에 따라 반복측정자료를 분석하는 통계적인 방법론으로는 Mixed Effect Model을 사용하는 방법과, Generalized Estimating Equations을 이용하는 방법 등이 존재합니다. 우선 Mixed Effect Model 중 Linear Mixed Effect Model에 대해 알아보도록 하겠습니다.

2. Fixed & Random Effects

모델에 대해 본격적으로 살펴보기 전에, fixed effects와 random effects에 대해 간단하게 알아보겠습니다.

Fixed Effects

우선 fixed effects의 경우, 우리가 관심 있는 전체 모집단의 평균적인 효과를 추정하는 모수입니다. 예를 들어, 다음과 같은 모델에 대해 생각해보겠습니다.

$y_{i j} = β_{0} + β_{1} * G r o u p_{i} + β_{2} * T i m e_{i j} + β_{3} (G r o u p_{i} * T i m e_{i j}) + ϵ_{i j}$

$y_{i j}$ 는 i번째 환자의 j번째 측정값을 의미하고, $ϵ_{i j} \sim N (0, σ^{2})$ 라고 가정하겠습니다.

위 모델에서

$β_{1}$ (Group effect): 시간에 관계없이 Group간 평균적인 차이를 나타냅니다. 이는 baseline에서 그룹 간 평균 차이를 의미합니다.

$β_{2}$ (Time effect): 어떤 Group에 속하는지와 무관하게 시간이 지남에 따라 전체 평균이 얼마나 변하는지를 나타냅니다. 이 값이 음수라면 평균이 감소, 양수면 평균이 증가하는 경향성을 의미합니다.

$β_{3}$ (Interaction effect): 시간과 Group 간의 상호작용 효과로, 시간에 따라 Group별로 나타나는 차이를 나타냅니다. 즉, 시간에 따라 Group 별로 얼마나 평균이 다르게 나타났는지를 의미합니다.

fixed effects의 경우 위와 같이 estimate 형태로 도출되고, 전체 데이터에 일관되게 적용됩니다.

Random Effects

random effects의 경우 fixed effects와 비슷하면서도 조금 다른 형태를 보입니다. 이번 포스팅에서는 intercept에만 random effect를 주어 여기에 집중해보도록 하겠습니다.

random effects는 이름 그대로 무작위적인 부분 즉, 확률적인 효과를 모델링 합니다. 앞서 언급했듯 fixed effects가 우리가 관심 있는 모집단 전체의 평균적인 효과를 추정하는 모수였다면, random effects는 그 모집단에서 추출된 개별 sample이나 cluster가 갖는 특징이나 변동성을 잡아내는 확률변수(random variable)입니다.

앞서 언급한 당뇨병 환자들을 대상으로 한 임상실험을 다시 생각해봅시다. 모든 환자들이 동일한 치료를 받더라도, 치료 효과가 나타나는 시작점이나(baseline)이나 시간에 따른 변화 정도는 환자 개개인마다 다를 수 있습니다. 이러한 이질성(heterogeneity)을 설명해주는 것이 random effects입니다. 즉, 각 subject가 전체적인 평균( $β_{0}$ )으로부터 얼마나 벗어나 있는지를 나타내는 확률변수로 활용됩니다. 아래 모델을 보겠습니다.

$y_{i j} = (u_{i} + β_{0}) + β_{1} * G r o u p_{i} + β_{2} * T i m e_{i j} + β_{3} (G r o u p_{i} * T i m e_{i j}) + ϵ_{i j}$

여기서 새로 추가된 $u_{i}$ 가 바로 i번째 환자의 random effect입니다. 이는 특정 상수가 아니라 $u_{i} \sim N (0, σ_{u}^{2})$ 을 따르는 확률변수이며, 모든 $i, j$ 에 대해 $ϵ_{i j}$ 과 독립이라고 가정합니다.

위 모델은 결국 다음을 가정합니다.

모든 환자들은 평균적인 시작점을 공유하지만, 각 환자는 자신만의 고유한 시작점(intercept)을 ( $u_{i} + β_{0}$ ) 갖습니다(random effects).
Group, Time, Interaction effects는 모든 환자에게 동일하하게 적용됩니다(fixed effects).

이처럼 모델에 random effects를 포함하게 되면, 단순하게 모든 환자의 평균적인 경향성만을 보는 것에서 한 발 더 나아가, 환자들 사이에 존재하는 변동까지 고려하게 되어 데이터를 훨씬 더 정교하게 설명할 수 있습니다.

결국 random effects가 표현하고자 하는 것은 ‘변동’에 초점을 맞추기에, fixed effects는 기댓값에 해당하는 $E [y_{i j}]$ 를 모델링 한다면, random effects는 분산(공분산)에 해당하는 $V a r (y_{i j})$ 를 모델링 합니다.

실제로 이 모델을 통해 i번째 환자의 첫번째와 두번째 시점의 공분산을 구해보면, $C o v (y_{i 1}, y_{i 2}) = C o v (u_{i} + ϵ_{i 1}, u_{i} + ϵ_{i 2})$ $= C o v (u_{i}, u_{i}) + C o v (u_{i}, ϵ_{i 1}) + C o v (u_{i}, ϵ_{i 2}) + C o v (ϵ_{i 1}, ϵ_{i 2})$ $= σ_{u}^{2}$ 으로 도출되므로, 같은 환자에게서 다른 시간대에 측정된 정보는 상관되어 있음을 확인할 수 있습니다.

3. Linear Mixed Effect Model with Random Intercept

자 이제 intercept에만 random effects가 주어진 상태에서 어떤 경우에 standard error가 같아지는지 알아보겠습니다. 이를 만족하기 위해서는 다음과 같은 엄격한 조건이 모두 성립해야합니다.

Complete Case: 주어진 데이터에서 결측이 단 하나라도 존재해서는 안 됩니다.
Balanced Data: 연구 디자인이 완벽하게 대칭적이고 균형 잡혀있어야 합니다. 반복측정자료에서 균형의 의미는 크게 두 가지를 의미합니다. 첫째는 관측 횟수의 균형으로 모든 환자가 동일한 횟수로 측정되어야하며, 둘째는 관측 시점의 균형으로, 모든 환자가 동일한 시점에 측정되어야 합니다.
Compound Symmtery: 모든 반복측정자료들 사이의 상관관계가 동일해야 합니다. 즉, 모든 측정 시점에서 각 데이터의 분산은 동일해야 하고 $(\forall i, j, V a r (Y_{i j}) = σ^{2} + {σ_{u}}^{2})$ , 각 환자 내에서 임의의 두 시점간 공분산 또한 항상 동일해야 합니다 $(C o v (Y_{i j}, Y_{i k}) = {σ_{u}}^{2}, where j \neq k)$ .

예시 데이터를 확인해보겠습니다. 참고로 예시 데이터의 값들은 모두 임의로 채워넣었기에 결과가 유의하고 않고는 의미가 없음을 미리 말씀드립니다.

library(openxlsx)
data <- openxlsx::read.xlsx("data.xlsx")
dim(data)

[1] 80  4

head(data,10)

   NO          Time Group     Value
1   1        Pre OP 65-79 14.786093
2   1  Post 1 month 65-79 10.477890
3   1 Post 6 months 65-79 10.196241
4   1   Post 1 year 65-79  7.049579
5   2        Pre OP 65-79  9.057036
6   2  Post 1 month 65-79 10.267466
7   2 Post 6 months 65-79  9.657085
8   2   Post 1 year 65-79 10.387103
9   3        Pre OP 65-79  8.508681
10  3  Post 1 month 65-79 10.999663

tail(data,10)

   NO          Time Group     Value
71 18 Post 6 months   ≥80  6.702226
72 18   Post 1 year   ≥80 10.030724
73 19        Pre OP   ≥80  9.320332
74 19  Post 1 month   ≥80  9.939527
75 19 Post 6 months   ≥80 10.699433
76 19   Post 1 year   ≥80 10.021359
77 20        Pre OP   ≥80  8.859701
78 20  Post 1 month   ≥80  9.440627
79 20 Post 6 months   ≥80 10.545564
80 20   Post 1 year   ≥80  8.391252

이 데이터는 총 20명의 환자를 각 4번씩 반복측정한 데이터입니다. NO열은 각 환자의 ID를 의미하고, Time은 총 4개로 수술 직후, 1개월 후, 6개월 후, 1년 후로 나뉘어 있습니다. Group의 경우 65~79세인 경우(NO가 1-10인 환자)와 80세 이상(NO가 11-20인 환자)인 두 가지 경우가 존재합니다. 이 데이터는 결측이 하나도 존재하지 않고 20명의 환자에 대해 각각 모두 같은 시점에 같은 횟수가 측정된 Complete & Balanced data입니다. 우리는 이러한 반복측정자료에서 Time과 Group에 대한 정보로 Value값에 대해 예측하고 싶어하는 상황입니다. Compounds Symmetry만 만족한다면, 앞서 말했던 것처럼 모든 S.E. 값이 동일하게 나올 것입니다. 한 번 확인해보겠습니다.

lme4패키지의 lmer 함수를 사용하여 진행해보면 아래와 같은 결과를 얻을 수 있습니다.

library(lme4)

Loading required package: Matrix

lmer_fit <- lmer(
  Value ~ Group*Time + (1|NO),
  data = data)
summary(lmer_fit)

Linear mixed model fit by REML ['lmerMod']
Formula: Value ~ Group * Time + (1 | NO)
   Data: data

REML criterion at convergence: 310.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2756 -0.4501  0.0719  0.4712  2.6606 

Random effects:
 Groups   Name        Variance Std.Dev.
 NO       (Intercept) 0.04915  0.2217  
 Residual             3.33376  1.8259  
Number of obs: 80, groups:  NO, 20

Fixed effects:
                             Estimate Std. Error t value
(Intercept)                    9.9769     0.5816  17.153
Group65-79                     0.2331     0.8225   0.283
TimePost 1 year               -0.4466     0.8165  -0.547
TimePost 6 months              0.4467     0.8165   0.547
TimePre OP                    -0.9778     0.8165  -1.197
Group65-79:TimePost 1 year     0.2448     1.1548   0.212
Group65-79:TimePost 6 months  -0.1802     1.1548  -0.156
Group65-79:TimePre OP          0.7250     1.1548   0.628

Correlation of Fixed Effects:
            (Intr) G65-79 TmPs1y TmPs6m TmPrOP G65-1y G65-6m
Group65-79  -0.707                                          
TimePost1yr -0.702  0.496                                   
TmPst6mnths -0.702  0.496  0.500                            
TimePre OP  -0.702  0.496  0.500  0.500                     
G65-79:TP1y  0.496 -0.702 -0.707 -0.354 -0.354              
G65-79:TP6m  0.496 -0.702 -0.354 -0.707 -0.354  0.500       
G65-79:TPOP  0.496 -0.702 -0.354 -0.354 -0.707  0.500  0.500

summary의 Random effects를 보면 intercept에 대해서만 적용된 모습을 확인할 수 있습니다. Fixed effects의 Std.Error을 보게 되면 time effect와 interaction effect에 해당하는 부분의 Std. Error가 모두 같은 모습을 확인할 수 있습니다. 어째서 이런 결과가 나타나는 것일까요?

Linear Regression Model

일반적인 선형 회귀 모델을 행렬 형태로 표현하면 다음과 같이 나타낼 수 있습니다.

$Y = X β + ϵ$

여기서 $ϵ \overset{iid}{\sim} N (0, σ^{2} I)$ 이라고 가정하는 것이 일반적입니다. 이후 최소제곱법(LSE)를 사용하여 $\hat{β} = (X^{T} X)^{- 1} X^{T} Y$ 라는 추정량을 얻을 수 있고, $V a r (\hat{β}) = σ^{2} (X^{T} X^{- 1})$ 라고 계산할 수 있습니다.

Linear Mixed Model

지금 우리가 다루고 있는 모델은 여기에 상관관계를 위해 $u$ 라는 term을 추가한 linear mixed model이고, 오차의 독립성 가정을 하지 않기에 선형 회귀 모델과는 조금 다른 구조와 분산-공분산 행렬을 가집니다. 이를 아래와 같이 표현할 수 있습니다.

$Y = X β + Z u + ϵ$

$u \overset{iid}{\sim} N (0, σ_{u}^{2} I), ϵ \overset{iid}{\sim} N (0, σ^{2} I)$ 이고, 이 둘은 독립이라고 가정하겠습니다.

i번째 환자 데이터의 분산-공분산 구조를 살펴보면 다음과 같습니다. $V a r (Y_{i}) = V_{i} = V a r (X_{i} β + Z_{i} u_{i} + ϵ_{i})$ $= V a r (Z_{i} u_{i} + ϵ_{i})$ $= V a r (Z_{i} u_{i}) + V a r (ϵ_{i}) (∵ ϵ ⊥ ⊥ u)$ $= Z_{i} V a r (u_{i}) Z_{i}^{T} + V a r (ϵ_{i})$

총 20명의 환자가 4번씩 반복측정한 상황이고, random effects가 intercept term에만 존재하기에 $u_{i}$ 는 univariate normal distribution을 따르는 확률변수가 되고, $Z_{i} = [\begin{matrix} 1 & 1 & 1 & 1 \end{matrix}]$ 라고 할 수 있으므로 $V_{i} = [\begin{matrix} 1 & 1 & 1 & 1 \end{matrix}] σ_{u}^{2} [\begin{matrix} 1 \\ 1 \\ 1 \\ 1 \end{matrix}] + σ^{2} [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]$

$= [\begin{matrix} σ^{2} + σ_{u}^{2} & σ_{u}^{2} & σ_{u}^{2} & σ_{u}^{2} \\ σ_{u}^{2} & σ^{2} + σ_{u}^{2} & σ_{u}^{2} & σ_{u}^{2} \\ σ_{u}^{2} & σ_{u}^{2} & σ^{2} + σ_{u}^{2} & σ_{u}^{2} \\ σ_{u}^{2} & σ_{u}^{2} & σ_{u}^{2} & σ^{2} + σ_{u}^{2} \end{matrix}]$ 라는 결과를 얻게 됩니다. 이를 통해 random intercept 모델은 compound symmetry 구조를 만족함을 확인할 수 있습니다.

자 그럼 이제 $β$ 의 추정치를 구해보겠습니다. 일반적인 Linear Regression Model에서 OLS(Ordinary Least Square)방법을 이용해 $\hat{β} = (X^{T} X)^{- 1} X^{T} Y$ 를 도출할 수 있다고 했고, 대략적인 과정은 다음과 같습니다.

$Y = X β + ϵ$

$Let Q (β) = (Y - X β)^{T} (Y - X β) = Y^{T} Y - 2 Y^{T} X β + β X^{T} X β$ $Then \arg min_{β} Q (β) = \hat{β} would be the vector that satisfies \frac{\partial Q (β)}{\partial β} = 0$

$since Q (β) is a convex function w.r.t β .$

$We know that X^{T} X is symmetric and by assuming it’s invertible,$

$\frac{\partial Q (β)}{\partial β} = - 2 X^{T} Y + 2 X^{T} X β = 0$ $X^{T} X β = X^{T} Y$

$∴ \hat{β} = (X^{T} X)^{- 1} X^{T} Y$ 자 이제는 이 방법을 살짝 비틀어서 Linear Mixed Model의 경우를 살펴보겠습니다. 일반적인 형태의 모델은 다음과 같았습니다.

$Y = X β + Z u + ϵ$

$Let Z u + ϵ = ϵ^{*}$ $Then Y = X β + ϵ^{*}$ $V a r (Y) = V$ 는 symmertic & positive-semi-definite matrix이므로 Cholesky Decomposition에 의해 다음을 만족하는 가역인 행렬 $Σ$ 가 존재합니다. $V = Σ Σ^{T}$ . 이제 $Σ^{- 1}$ 을 위 모델의 양변에 곱해봅시다. $Σ^{- 1} Y = Σ^{- 1} X β + Σ^{- 1} ϵ^{*}$ 이제 아래와 같이 변수들을 새로 정의하면, $Σ^{- 1} Y = \tilde{Y},$ $Σ^{- 1} X β = \tilde{X} β$

$Σ^{- 1} ϵ^{*} = \tilde{ϵ}$

아래와 같이 OLS의 경우와 아주 유사한 형태의 모델을 새로 정의할 수 있습니다. $\tilde{Y} = \tilde{X} β + \tilde{ϵ}$ 이러한 방법을 OLS를 일반화했다고 하여 Generalized Least Square(GLS)라고 부릅니다. 이제 OLS와 형태가 같기 때문에 $\hat{β}$ 의 추정량은 다음와 같이 도출됩니다.

${\hat{β}}_{G L S} = ({\tilde{X}}^{T} \tilde{X})^{- 1} {\tilde{X}}^{T} \tilde{Y}$

$= ((Σ^{- 1} X)^{T} Σ^{- 1} X)^{- 1} (Σ^{- 1} X)^{T} Σ^{- 1} Y$

$= (X^{T} {Σ^{T}}^{- 1} Σ^{- 1} X)^{- 1} X^{T} {Σ^{T}}^{- 1} Σ^{- 1} Y$

$= (X^{T} (Σ^{T} Σ)^{- 1} X)^{- 1} X^{T} (Σ^{T} Σ)^{- 1} Y$

$= (X^{T} (Σ Σ^{T})^{- 1} X)^{- 1} X^{T} (Σ Σ^{T})^{- 1} Y (∵ Σ^{T} = Σ)$

$= (X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} Y$ 임을 알 수 있습니다. 그렇다면 분산-공분산 행렬은 다음과 같이 유도될 것입니다.

$V a r ({\hat{β}}_{G L S}) = V a r ((X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} Y)$

$= (X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} V a r (Y) ((X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1})^{T}$

$(X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} V {V^{- 1}}^{T} X {(X^{T} V^{- 1} X)^{T}}^{- 1}$

$(X^{T} V^{- 1} X)^{- 1} (X^{T} V^{- 1} X) (X^{T} V^{- 1} X)^{- 1} (∵ V V^{- 1} = I, {V^{- 1}}^{T} = V^{T})$

$= (X^{T} V^{- 1} X)^{- 1}$ 앞에서 보인것 처럼 $V$ 는 Compound Symmertry 구조이고, Complete & Balanced한 데이터를 다루고 있기에 Design matrix인 X의 구조도 변하지 않습니다. 결국 standard error을 결정짓는 것은 방금 전 구한 $V a r ({\hat{β}}_{G L S})$ 의 대각성분과 관련이 있을텐데, 이 값이 같기 때문에 fixed effects의 standard error가 같게 나오는 것입니다.

Why Linear Mixed Model Is Better?

지금까지 잘 읽어오시다가 몇몇 분들은 이런 생각이 드셨을 수도 있습니다.

반복측정자료고 뭐고 그냥 환자 NO를 dummy variable로 만들어서 fitting해도 되지 않나? 굳이 왜 random effects를 쓰지? 이렇게 하면 fixed effects만으로도 표현이 가능한데…

네 맞습니다. 가능합니다. 한 번 그렇게 진행해볼까요?

glm_fit <- glm(
  Value ~ NO + Group + Time + Group*Time,
  data = data)
summary(glm_fit)


Call:
glm(formula = Value ~ NO + Group + Time + Group * Time, data = data)

Coefficients: (1 not defined because of singularities)
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  10.67448    1.04090  10.255 2.78e-14 ***
NO10                          0.43968    1.29108   0.341   0.7348    
NO11                         -0.24069    1.47206  -0.164   0.8707    
NO12                         -0.88679    1.47206  -0.602   0.5494    
NO13                          0.48099    1.47206   0.327   0.7451    
NO14                         -0.71105    1.47206  -0.483   0.6310    
NO15                          0.28560    1.47206   0.194   0.8469    
NO16                         -2.18372    1.47206  -1.483   0.1438    
NO17                          0.75654    1.47206   0.514   0.6094    
NO18                         -2.92089    1.47206  -1.984   0.0523 .  
NO19                         -0.43488    1.47206  -0.295   0.7688    
NO2                          -0.78528    1.29108  -0.608   0.5456    
NO20                         -1.12076    1.47206  -0.761   0.4498    
NO3                          -0.50924    1.29108  -0.394   0.6948    
NO4                          -1.20306    1.29108  -0.932   0.3556    
NO5                          -0.20593    1.29108  -0.160   0.8739    
NO6                          -1.05182    1.29108  -0.815   0.4188    
NO7                          -0.05504    1.29108  -0.043   0.9662    
NO8                          -1.48829    1.29108  -1.153   0.2541    
NO9                           0.21424    1.29108   0.166   0.8688    
Group65-79                         NA         NA      NA       NA    
TimePost 1 year              -0.44664    0.81655  -0.547   0.5866    
TimePost 6 months             0.44668    0.81655   0.547   0.5866    
TimePre OP                   -0.97777    0.81655  -1.197   0.2364    
Group65-79:TimePost 1 year    0.24482    1.15478   0.212   0.8329    
Group65-79:TimePost 6 months -0.18020    1.15478  -0.156   0.8766    
Group65-79:TimePre OP         0.72501    1.15478   0.628   0.5328    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 3.333764)

    Null deviance: 260.10  on 79  degrees of freedom
Residual deviance: 180.02  on 54  degrees of freedom
AIC: 345.91

Number of Fisher Scoring iterations: 2

Group 변수에서 NA가 나왔습니다. 아마 다중공산성 문제 등이 발생해서 design matirx가 singular 해진 것 같습니다. 여기서부터 문제가 느껴지시나요? NO의 level이 너무 많아지게되어서 모델 구조가 필요 이상으로 복잡해지고, 원하지도 않은 수많은 계수들을 추정하게 됩니다. 이렇게 되면 정작 우리가 중요하게 생각하고 있는 변수들의 효과에 대한 해석도 어려워집니다.

Group 변수에서 NA가 발생한 것 자체가 추정의 불안정성을 보여주는 예시라고 볼 수 있겠습니다. 이와 별개로 만약 반복측정된 횟수가 정말 적은 경우라면, NA가 뜨지 않더라고 해당 변수에 해당하는 표준 오차가 굉장히 커질 것이기에 값 자체에 대한 신뢰성 또한 떨어집니다.

이것이 우리가 lmer을 쓰는 이유 중 하나입니다. fixed effects만으로 처리하기에는 어려운 부분을 random effects로 말끔하게 처리하여 문제를 해결해준다는 점에서 큰 장점을 가집니다.

4. GEE VS lmer

GEE

이번에는 GEE(Generaluzed Estimating Equation)를 사용해 똑같은 데이터를 분석해보겠습니다. 코드는 아래와 같습니다. 참고로 corsrt = exchangeable이라고 설정해주어 한 subject내에서 측정된 모든 관측지의 상관관계는 동일하다고 가정했습니다. GEE는 심지어 반복측정환자들끼리 잘 정리되어 있는 데이터가 아니라면 제대로 돌아가지 않습니다. 예시로 사용한 데이터는 반복측정된 환자 순서대로 잘 정리되어 있음을 밝힙니다.

library(geepack)
gee_mod <- geeglm(
  Value ~ Group*Time,
  id        = NO,
  data      = data,
  corstr    = "exchangeable"
)
summary(gee_mod)


Call:
geeglm(formula = Value ~ Group * Time, data = data, id = NO, 
    corstr = "exchangeable")

 Coefficients:
                             Estimate Std.err    Wald Pr(>|W|)    
(Intercept)                    9.9769  0.8306 144.290   <2e-16 ***
Group65-79                     0.2331  0.9386   0.062    0.804    
TimePost 1 year               -0.4466  0.9007   0.246    0.620    
TimePost 6 months              0.4467  0.4787   0.871    0.351    
TimePre OP                    -0.9778  0.8964   1.190    0.275    
Group65-79:TimePost 1 year     0.2448  1.2324   0.039    0.843    
Group65-79:TimePost 6 months  -0.1802  0.6812   0.070    0.791    
Group65-79:TimePre OP          0.7250  1.3258   0.299    0.584    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)    3.045  0.7025
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha  0.01453 0.07475
Number of clusters:   20  Maximum cluster size: 4

lmer을 사용했을 때와 다른 점을 느끼셨나요? 두 모델 모두 estimate들의 값은 같게 나오지만, standard error의 값이 다르게 나온다는 점입니다. 두 방법 모두 반복측정자료를 분석하는데 사용되는데, 대체 어떤 차이가 있는 것일까요?

저는 To GEE or Not to GEE: Comparing Population Average and Mixed Models for Estimating the Associations Between Neighborhood Risk Factors and Health라는 논문에서 어느 정도 해답을 찾을 수 있었습니다.

lmer과 같은 Mixed Model을 neighbor-specific model로 설명합니다. 이름에서 알 수 있듯이, 특정한 지역(neighbor or cluster) 혹은 객체(subject)가 각자 자신만의 고유한 반응 패턴을 가질 수 있다고 가정합니다. 결국 앞서 계속 살펴봤던 random effects model과 같은 말입니다. 각 지역(혹은 개인)마다 다른 intercept나 slope를 가질 수 있다고 가정하여 개별 단위 수준에서 추정을 진행합니다. 각 환자마다 당뇨병 치료의 영향이 다를 수 있음을 인정하는 것입니다.

반면, GEE와 같은 모델들은 Population-Average Model로 설명합니다. 이 경우, 각 지역이나 객체의 차이보다는 전체 모집단에서 나타나는 평균적인 관계에 더 주목합니다. 당뇨병 치료에서 어떤 요인이 한 단위 증가할 때, 어떤 환자인지는 중요하지 않고, 평균적으로 얼마나 인슐린 농도가 변하는지를 탐구하고자 하는 것입니다.

논문에 따르면, linear model에서 상당히 흥미로운 부분을 찾을 수 있습니다.

A simple random intercept linear model implies equal variances for all observations and equal covariances of all possible paired observations within the statistical unit (neighborhood) … This will yield the same estimates as the exchangeable working correlation model in GEE.

random intercept만을 포함하는 lmer 모델은 compound symmertry를 만족하기에 모든 관측치의 분산이 같고, 한 subject내의 모든 관측치는 동일한 상관관계를 공유합니다. 나아가 이는 GEE에서 상관구조를 ’exchangeable’로 설정했을 때의 가정과 정확히 일치합니다. 이러한 이유에서 이 둘은 동일한 회귀 계수 추정치를 내놓았던 것입니다.

그러나 통계적인 추론 과정에서 사용되는 standard error는 달라집니다. lmer의 경우 앞서 $u \sim N (0, σ_{u}^{2})$ 라고 가정한 것처럼, random effects의 분포를 가정합니다. 당연히 이 가정이 맞다는 전제에서 MLE(Maximum Likelihood Estimation)나 REML(Restrictied Maximum Likelihood)을 이용해 standard error를 계산합니다. 모델의 구조적 가정이 그대로 반영되기 때문에, 모든 시점에서 동일한 표준 오차를 내놓을 수밖에 없는 것입나다.

GEE의 경우 lmer과 달리 random effects의 확률적인 분포를 가정하지 않습니다. 관측된 데이터를 기반으로 robust한 standard error을 계산하는데요, 이를 sandwich estimator라고 부릅니다. 이 방식을 사용하게 되면 분포 가정에 있어 비교적 자유롭고, 설령 설정한 상관구조가 조금 잘못되었더라도 비교적 정확한 표준 오차를 제공합니다.

BLUE vs BLUP

BLUE(Best Linear Unbiased Estimator): fixed effects에 대한 추정량입니다. 우리가 알고자 하는 모집단 전체의 평균적인 경향성을 나타내는데, 예를 들어 평균적으로 치료 수준을 한 단위 높였을 때 인슐린 농도가 어느 정도 변하는가?에 대한 답입니다. lmer과 GEE 모두 fixed effects가 존재하기에, 두 모델에서 모두 구할 수 있습니다.

BLUP(Best Linear Unbiased Predictor): random effects에 대한 예측치입니다. BLUE와 다르게 추정이 아니라 예측이라는 점에 주목해야 합니다. BLUE가 모집단 전체의 평균적인 경향성에 주목했다면, BLUP는 모집단 평균에서 벗어난 각 개인의 고유한 편차를 나타냅니다. 예컨대, 3번째 환자는 평균보다 얼마나 더 인슐린 농도가 늘었을까?에 대한 답입니다. random effects에 대한 부분이기에, population-average에 해당하는 GEE에서는 구할 수 없고, lmer에서만 구할 수 있습니다. 한 번 살펴 보겠습니다.

증명의 기초는 Henderson이 1959년에 발명한 Mixed Model Equations(MME)를 따라갑니다. 증명이 너무 길고 이해가 어려울 수도 있어 복잡한 계산 과정은 과감히 생략하도록 하겠습니다. 행렬형태의 식으로 표현하기 위해, 앞선 노테이션과 조금 다르게, $u$ 의 분산-공분산 행렬을 $G$ 로, $ϵ$ 의 분산-공분산 행렬은 $R$ 이라고 하겠습니다.

lmer은 likelihood를 기반으로 작동하는 만큼, 이번에는 likelihood로 전개를 해보겠습니다.(사실 LS를 쓰나 Likelihood를 쓰나 결국 똑같은 형태가 나옵니다.) n은 전체 sample size를, r은 random effects의 차원을 의미합니다. 우선 다음과 같이 conditional 형태로 likelihood를 적어봅시다.

$f (Y, u) = f (Y ∣ u) \cdot f (u)$ $= \frac{1}{(2 π)^{\frac{n}{2}} det R^{1 / 2}} \exp (- \frac{1}{2} (Y - X β - Z u)^{T} V^{- 1} \frac{1}{2} (Y - X β - Z u))$

$\times \frac{1}{(2 π)^{\frac{r}{2}} det G^{1 / 2}} \exp (- \frac{1}{2} u^{T} G u)$

이제 여기에 로그를 씌우게 되면 $β$ 와 $u$ 모두에 대해 2차식으로 볼 수 있습니다. likelihood function의 concave함을 이용해 각각 $β$ 와 $u$ 모두로 편미분을 해주고 각각 hat을 씌워주면, 아래와 같이 나타낼 수 있습니다.

$\begin{aligned} X^{T} R^{- 1} X \hat{β} + X^{T} R^{- 1} Z \hat{u} & = X^{T} R^{- 1} Y \\ Z^{T} R^{- 1} X \hat{β} + (Z^{T} R^{- 1} Z + G^{- 1}) \hat{u} & = Z^{T} R^{- 1} Y \end{aligned}$ 이제 이를 행렬 형태로 나타내면, 아래와 같이 쓸 수 있겠습니다.

$[\begin{matrix} X^{T} R^{- 1} X & X^{T} R^{- 1} Z \\ Z^{T} R^{- 1} X & Z^{T} R^{- 1} Z + G^{- 1} \end{matrix}] [\begin{matrix} \hat{β} \\ \hat{u} \end{matrix}] = [\begin{matrix} X^{T} R^{- 1} Y \\ Z^{T} R^{- 1} Y \end{matrix}]$ 위 행렬식이 MME의 기본적인 형태입니다. 이를 정말 열심히 열심히 정리하면 다음과 같은 BLUE와 BLUP를 얻을 수 있습니다. $B L U E = \hat{β} = (X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} Y : GLS estimator$ $B L U P = \hat{u} = G Z^{T} V^{- 1} (Y - X \hat{β})$

여기서 한 가지 의문이 들 수 있는데요, 일반적으로 Gauss-Markov Theorem을 만족할 때 OLS를 통해 구한 estimator가 unbiased하고 가장 작은 분산을 갖기에 BLUE라고 부릅니다. 반복측정자료는 같은 subject에서 측정된 데이터 간 상관관계가 존재하기 때문에 일반적으로 Gauss-Markov Theorem을 만족하지 않습니다. 그렇다면 과연 lmer을 통해 구한 ${\hat{β}}_{G L S}$ 는 BLUE가 맞을까요? 놀랍게도 맞습니다. 간단하게만 알아보면, Least Square을 일반화해서 GLS를 진행한 것처럼, Gauss-Markov Theorem을 일반화한 Aitken’s Theorem이 존재합니다! 이 정리를 통해 GLS estimator가 BLUE임을 보일 수 있습니다. 간단한 증명은 다음과 같습니다.

우선 unbiasedness부터 보면, $E [{\hat{β}}_{G L S}] = E [(X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} (X β + ϵ^{*})]$

$= E [(X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} X β] + E [(X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} ϵ^{*}]$

$= E [β] + (X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} E [ϵ^{*}] = β$ 이렇게 쉽게 unbiasedness를 보일 수 있습니다. 나아가 분산을 보겠습니다. $b$ 를 또다른 linear unbiased estimator로 이면서 다음과 같은 형태를 만족한다고 해봅시다. $b = [(X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} + A] Y$ unbiasedness를 만족해야 하기 때문에, 계산을 해보면 $A X = 0$ 이 되어야 합니다. 그렇다면 분산을 구해보면, $V a r (b) = [(X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} + A] V [(X^{T} V^{- 1} X)^{- 1} X^{T} V^{- 1} + A]^{- 1}$ $= (X^{T} V^{- 1} X)^{- 1} + A V A^{T} + (X^{T} V^{- 1} X)^{- 1} X^{T} A^{T} + A X (X^{T} V^{- 1} X)^{- 1}$ $= (X^{T} V^{- 1} X)^{- 1} + A V A^{T} (∵ A X = 0)$ $\geq (X^{T} V^{- 1} X)^{- 1} (∵ V is psd)$ $= V a r ({\hat{β}}_{G L S})$

따라서 ${\hat{β}}_{G L S}$ 가 BLUE가 됨을 보일 수 있습니다. BLUP 또한 BLUE와 비슷한 방법으로 또다른 linear unbiased parameter을 들고 와서 증명이 가능합니다. BLUP의 형태를 다시 한 번 보겠습니다.

$B L U P = \hat{u} = G Z^{T} V^{- 1} (Y - X \hat{β})$ 이제 이 식을 우리가 계속 다뤘던 intercept에만 random effects가 존재하는 모델에 맞추어 바꿔보도록 하겠습니다. $G$ 또한 원래의 노테이션으로 바꿔서 i번째 환자의 입장에서 정리해보면 다음과 같이 표현할 수 있겠습니다.

${\hat{u}}_{i} = (\frac{n_{i} σ_{u}^{2}}{n_{i} σ_{u}^{2} + σ^{2}}) \times \frac{1}{n_{i}} \sum_{j = 1}^{n_{i}} (y_{i j} - x_{i j}^{T} \hat{β})$ 이렇게 간단한게 스칼라 형태로 정리 되는데요, 이 식이 가지는 의미는 무엇일까요?

우선 뒷부분부터 보면 우리가 평소에 회귀분석에서 볼 수 있는 잔차의 평균 형태입니다. i번째 환자가 평균적으로 예측에서 얼마나 벗어나있는지에 대해 알려주는 정보로 해석할 수 있습니다. 또한 일반적인 선형 회귀 모델과 다르게 우리 모델은 $Z u$ 라는 random effects를 표현하는 term이 존재했으므로, 이 부분을 fixed effects의 효과를 제거하고 random effects와 $ϵ$ 에 해당하는 error만 남긴 부분으로 볼 수도 있습니다.

첫번째 부분이 BLUP의 핵심이라고 볼 수 있습니다. 이는 가중치의 형태로 수축(Shrinkage) 효과를 보여주는데요, 직관적으로 $n_{i}$ 의 값이 커지거나 $σ_{u}^{2}$ 커질수록 이 값은 1에 가까워집니다. 각 환자의 데이터가 많거나 환자 내에서 측정한 데이터끼리의 상관관계가 강할수록 가중치가 1에 가까워져 각 환자의 평균 잔차를 더 신뢰하는 형태가 됩니다. 반대로, 0에 가까울수록 각 환자의 데이터보다는 전체 평균에 가까워지도록 예측치를 수축시킵니다.

결국 BLUP는 각 환자의 반복측정된 데이터의 수와 그 데이터끼리 상관된 정도에 따라 어떻게 예측을 진행할지를 결정해준다고 볼 수 있습니다. 데이터가 부족한 상황에서는 예측이 불확실한 노이즈에 과도하게 영향을 받는 것을 막고, 전체 데이터의 정보에 따라 안정적이고 보수적인 예측을 가능하게 해준다는 점에서 의미가 있습니다.

lmer은 모집단의 전반적인 경향성뿐만 아니라 그 안에서 개개인이 얼마나 다양한 패턴을 보여주는지 또한 제시할 수 있다는 점에서 GEE와는 다른 차별점을 가진다고 할 수 있겠습니다.

마치며

결국 lmer에서 intercept에만 random effects가 존재할 때 추정치들의 표준 오차가 같아지는 현상은 무언가 잘못된 것이 아니라, 모델 구조상 자연스러운 결과입니다. 반복측정자료 연구에서 주된 관심사가 개체간 변동성을 고려하는 즉, subject(neighbor)-specific하다면, 이는 받아들일 수 있는 결과입니다.

다만, 전체 집단의 평균적인 변화에 관심이 있다면 GEE와 같은 population-average 모델이 더 좋은 방법이 될 수 있습니다. random effects의 분포를 가정하지 않으면서 robust한 추정량을 제공하는 GEE가 대안이 될 수도 있습니다.

결국 가장 중요한 것은 내가 지금 무엇을 하고 있는지에 대해 정확히 아는 것입니다. 이 질문에 대해 정확히 대답할 수 있어야 이에 따라 적합한 모델과 그 결과에 대한 올바른 해석을 할 수 있을 것입니다.

Reuse

CC BY-NC 4.0

Citation

BibTeX citation:

@online{jee2025,
  author = {Jee, Mingu},
  title = {반복측정자료 {모델링:} Linear Mixed Effects Model},
  date = {2025-08-12},
  url = {https://blog.zarathu.com/posts/2025-08-08-lmer and gee/},
  langid = {en}
}

For attribution, please cite this work as:

Jee, Mingu. 2025. “반복측정자료 모델링: Linear Mixed Effects Model.” August 12, 2025. https://blog.zarathu.com/posts/2025-08-08-lmer and gee/.