Linear mixed effect model

김당찬·2023년 6월 27일
0

Fixed and Random effect

ii번째 그룹에 대한 jj번째 관측값 y11,,ynny_{11},\ldots,y_{nn} 들이 주어질 때, 이들의 평균을 모델링하는 다음과 같은 모형을 생각해보자.

E[yij]=μ+αi+ϵij\mathrm{E}[y_{ij}] = \mu + \alpha_{i}+\epsilon_{ij}

이때 각 parameter μ,α\mu,\alpha 는 평균에 영향을 미치는 모수이고 μ\mu는 global mean에, α\alpha는 각 그룹의 평균에 영향을 미친다. ϵij\epsilon_{ij} 는 개별 관측값에 영향을 미치는 오차항이다. 여기서 두 가지 경우를 살펴보도록 하자.

1. α\alpha is a fixed effect

α\alpha가 고정효과라는 것은, 각 그룹에 대한 평균의 영향이 그룹 내에서는 모두 일정하다는 것을 의미한다. 예를 들어 n개의 환자 그룹에 서로 다른 n가지 치료방법을 적용한다 했을 때, 그룹 내에서는 일정한 효과가 있다고 가정하고 그룹 간 효과는 유의미하게 차이가 있다고 설정한 상황을 생각해 볼 수 있다. 이 경우 n개의 parameter αi\alpha_i를 설정하게 되며 이를 고정효과 모형이라고 한다.

2. α\alpha is a random effect

α\alpha가 임의효과(random effect)라는 것은, 각 그룹 간의 차이가 일정하다고 설정하기 보단, 일종의 변동(variation)으로 간주하는 것을 의미한다. 앞선 예시에서와 같이 n개의 환자 그룹이 존재하는데, 이번에는 이들을 각각 서울에 있는 임의의 n개의 (동급의) 일반병원에서 치료하는 상황을 고려해보자. 이 경우는 고정효과와는 다르게 각 병원의 차이가 실험적으로 유의미하지 않으며, 임의로 선택한 상황이므로 개별 그룹 간 효과 역시 임의인 상황이다. 이런 경우 αi\alpha_i 를 random variable로 보게 되며, 이를 임의효과 모형이라고 한다. 즉 다음과 같은 상황을 생각해볼 수 있다.

yij=μ+αi+ϵijαiiidN(0,σa2)ϵijiidN(0,σ2)αiϵij\begin{aligned} y_{ij}&= \mu+\alpha_i+\epsilon_{ij}\\ \alpha_{i}&\overset{iid}{\sim} N(0,\sigma_{a}^{2})\\ \epsilon_{ij}&\overset{iid}{\sim} N(0,\sigma^{2})\\ \alpha_i&\perp\epsilon_{ij} \end{aligned}

Mixed effect model

Mixed effect model이란 앞서 다룬 고정효과와 임의효과가 모두 포함되어있는 모형을 말한다. 다음과 같이 일반적인 모형을 살펴보자.

Y=Xβ+Zγ+eY=X\beta+Z\gamma+e

여기서 X,ZX,Z는 알려진 design matrix이고 β\beta는 고정효과 벡터를, γ\gamma는 임의효과 벡터를 각각 나타낸다. 이때 임의효과 벡터의 경우 확률벡터의 일종이므로, 공분산행렬과 평균벡터를 지정해주어야 하는데, E(γ)=0,cov(γ)=:D,cov(γ,e)=0,cov=(e)=R\mathrm{E}(\gamma)=0, cov(\gamma)=:D,cov(\gamma,e) = 0, cov=(e)=R 로 두도록 하자. 그러면 관측값과 임의효과의 결합분포는 다음과 같이 주어진다.

(Yγ)((Xβ0),(ZDZT+RZDDZTD))\begin{pmatrix}Y\\ \gamma \end{pmatrix} \sim \begin{pmatrix} \begin{pmatrix}X\beta \\0 \end{pmatrix}, \begin{pmatrix}ZDZ^T+R&ZD\\ DZ^{T}& D \end{pmatrix} \end{pmatrix}

그런데 실제로 각 모수(행렬)들을 추정하는 과정에서 임의효과의 공분산행렬 DD가 복잡한 경우, 즉 그룹간 임의효과의 상관관계가 존재하는 경우 해당 성분의 추정이 어려워지는 문제가 발생하기 때문에, Variance component model을 활용하여 공분산행렬들에 구조를 주어 추정가능하게끔 변환하는 과정을 거친다. 여기서는 cov(γ)=D=diag{σi2I}cov(\gamma)=D=diag\{\sigma_{i}^{2}I\} 형태를 주어 cov(Y)=V=lσl2ZlZlT+Rcov(Y)=V=\sum_{l}\sigma_{l}^{2}Z_{l}Z_{l}^{T}+R 로 변환하고, 이를 이용해 최대가능도 추정량을 구할 수 있다.

Mixed model equation

XβX\beta의 best linear unbiased estimator(BLUE)와 ZγZ\gamma의 best linear unbiased predictor(BLUP)은 다음과 같은 손실함수를 최소화하는 β,γ\beta,\gamma를 찾으면 구할 수 있는데, 이를 mixed model equation이라고 한다.

minβ,γ    (YXβZγ)TR1(YXβZγ)+γTD1γ\min_{\beta,\gamma}\;\; (Y-X\beta-Z\gamma)^{T}R^{-1}(Y-X\beta-Z\gamma)+\gamma^{T}D^{-1}\gamma

이는 Y,γY,\gamma의 joint likelihood와 동일한데, 일반적으로 임의효과인 γ\gamma는 관측이 불가능하므로 실제 가능도를 상정하기 어렵지만, 여기서는 joint distribution이 존재한다고 가정하여 다음과 같은 분해를 이용한 것이다.

f(Y,γβ)=f(Yγ,β)f(γβ)f(Y,\gamma|\beta)=f(Y|\gamma,\beta)f(\gamma|\beta)

Estimation example

처음 다루었던 예시 모형을 살펴보도록 하자.

yij=μ+αi+ϵijy_{ij}=\mu +\alpha_{i}+\epsilon_{ij}

이 경우에 α\alpha가 고정효과인지, 임의효과인지에 따라 추정 형태가 어떻게 달라지는지 살펴보자. 먼저 α\alpha가 고정효과인 경우, 다음과 같이 최대가능도 추정량을 구할 수 있다.

α^i=yˉi.yˉ..\hat\alpha_{i}=\bar{y}_{i.}-\bar{y}_{..}

여기서 yˉ..,yˉi.\bar y_{..}, \bar y_i. 는 각각 전체 관측치의 평균과 i번째 그룹 내 관측치의 평균을 의미한다. 반면, α\alpha를 임의효과로 가정하는 경우는 다음과 같이 주어진다(mixed model equation을 풀면 된다).

α~i=σ^a2σ^a2+σ^2/ni(yˉi.yˉ..)\tilde\alpha_{i}=\frac{\hat\sigma_{a}^{2}}{\hat\sigma_{a}^{2}+\hat\sigma^{2}/n_{i}}(\bar y_{i.}-\bar y_{..})

이를 고정효과의 추정 결과와 비교해보면, additional term σ^a2σ^a2+σ^2/ni\frac{\hat\sigma_{a}^{2}}{\hat\sigma_{a}^{2}+\hat\sigma^{2}/n_{i}} 가 1보다 작으므로 추정량의 수축(shrinkage)이 발생했다고 볼 수 있다. 다만, 그룹 간 분산 σa2\sigma_{a}^{2}이 커질수록 그 수축효과가 작아지는 것 역시 확인가능한데, 이는 그룹간 분산이 클수록 임의효과가 random effect와 같이 작용한다는 것을 의미한다. 좀 더 생각해보면, 그룹간 분산을 바탕으로 그룹 간 평균의 모델링을 임의효과로 설정할지, 고정효과로 설정할지 판단할 수 있다는 것을 의미한다.

profile
블로그 이사했습니다 https://ddangchani.github.io

0개의 댓글