13  일반화 선형 혼합 효과 모형 (generalized linear mixed effect models, GLMM)

Acknowledgement: 본 절의 구성에는 다음 교재를 발췌하였다.

주 출처

• 김양진 (2020). R과 SAS를 활용한 경시적 자료분석. 자유아카데미. - 6장 일반화 선형혼합효과모형

보조 출처

• 이근백 (2022). 경시적 자료분석 - R 활용. 자유아카데미. - 5장 연속형 경시적 자료분석을 위한 선형혼합모형

13.1 일반화 선형 혼합 모형 이론

13.1.1 개요

• 이전 단원의 선형혼합모형(Linear mixed effect model, LMM)의 outcome을 binary(0/1), count형 등으로 확장하려 한다.
• $Y_{ij}$가 연속형인 경우의 선형모형에서는 오차항 ${ϵ}_{ij}$의 공분산구조(예를 들어 exchangeable, AR(1))를 통해 반복치들 간의 연관관계를 모형화할 수 있었다. 그러나 로지스틱 회귀모형을 비롯한 일반화 선형모형(GLM)의 모형식은 오차항을 포함하지 않는다. 개체별로 공유되는 효과는 random effect로 표현해야 하겠다.

13.1.2 모형식과 가정

• 각 개체의 반응값을 ${\mathbf{y}}_i=(y_{i1},\cdots ,y_{i{m}_i})$라 하자. ($i=1,\dots ,N$)
• Recall: linear mixed effect model (LMM) $$y_{ij}={\mathbf{x}}_{ij}^T\mathit{\beta }+{\mathbf{z}}_{ij}^T{\mathbf{u}}_i+{ϵ}_{ij},j=1,\dots ,m_i,i=1,\cdots ,N$$
• 여기서 고정효과 관련 공변량 ${\mathbf{x}}_{ij}$, 고정효과 계수 $\mathit{\beta }$, 랜덤효과 관련 공변량 ${\mathbf{z}}_{ij}$은 non-random으로 가정된다. 랜덤효과 ${\mathbf{u}}_i$는 확률변수로 가정된다. 만약 $E(y_{ij}|{\mathbf{u}}_i)={\mu }_{ij}$로 쓰면, LMM은 다음과 같이 다시 쓸 수 있었다. $${\mu }_{ij}={\mathbf{x}}_{ij}^T\mathit{\beta }+{\mathbf{z}}_{ij}^T{\mathbf{u}}_i,j=1,\dots ,m_i,i=1,\cdots ,N$$
• GLMM의 모형식: known link function $g:\mathbb{R}\to \mathbb{R}$에 대하여, $$g({\mu }_{ij})=g(E(y_{ij}|{\mathbf{u}}_i))={\mathbf{x}}_{ij}^T\mathit{\beta }+{\mathbf{z}}_{ij}^T{\mathbf{u}}_i,j=1,\dots ,m_i,i=1,\cdots ,N,$$
• GLMM의 모형식 (벡터 표현): $$g(E({\mathbf{y}}_i|{\mathbf{u}}_i))={\mathbf{X}}_i\mathit{\beta }+{\mathbf{Z}}_i{\mathbf{u}}_i,i=1,\cdots ,N$$
• (위 식에서 $g(\cdot )$는 elementwise하게 적용된다고 가정하였다)
• 확률 가정 $$\begin{array}{c}E({\mathbf{u}}_i)=\mathbf{0},\mathrm{V}\mathrm{a}\mathrm{r}({\mathbf{u}}_i)={\mathbf{R}}_i,\\ {\mathbf{u}}_i\mathrm{i}\mathrm{i}\mathrm{d}\\ y_{i1},\dots ,y_{i{m}_i}\text{ are independent given }{\mathbf{u}}_i,\end{array}$$

13.1.3 예

• 이항 경시적 자료(binary categorical data): $$\begin{array}{c}{y}_{ij}\mid {\mathbf{u}}_i\sim B(1,{\mu }_{ij});\\ (\text{i.e.},P(y_{ij}=1\mid {\mathbf{u}}_i)=:{\mu }_{ij},P(y_{ij}=0\mid {\mathbf{u}}_i)=:1-{\mu }_{ij})\\ \text{where }\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}({\mu }_{ij}):=\log \left(\frac{{\mu }_{ij}}{1-{\mu }_{ij}}\right)={\mathbf{x}}_{ij}^T\mathit{\beta }+{\mathbf{z}}_{ij}^T{\mathbf{u}}_i,j=1,\dots ,m_i,i=1,\cdots ,N.\end{array}$$
• 계수형 경시적 자료(count data): $$\begin{array}{c}{y}_{ij}\mid {\mathbf{u}}_i\sim \mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{n}({\mu }_{ij})\\ \log ({\mu }_{ij})={\mathbf{x}}_{ij}^T\mathit{\beta }+{\mathbf{z}}_{ij}^T{\mathbf{u}}_i,j=1,\dots ,m_i,i=1,\cdots ,N.\end{array}$$

13.1.4 분석 예제: ICHS (호흡기 질환과 비타민 결핍의 관계)

• 먼저 자료를 불러오자.

ichs = read.table("data/ichs.dat", header=T)

dim(ichs)

[1] 1500    6

head(ichs)

  ID RESPONSE TIME GENDER VITA AGE
1  1        0    0      1    0   4
2  1        0    3      1    0   4
3  1        0    6      1    0   4
4  1        1    9      1    0   4
5  1        0   12      1    0   4
6  1        0   15      1    0   4

table(ichs[ichs$TIME==0, ]$VITA)

  0   1 
159  91 

table(ichs$RESPONSE, ichs$TIME)

   
      0   3   6   9  12  15
  0 183 182 172 176 174 169
  1  67  68  78  74  76  81

• 호흡기 질환과 비타민 결핍과의 관계를 조사하기 위한 자료이다.
• 랜덤 절편(random intercept)을 포함한 GLMM을 고려해 보자.

$$\text{logit}\left\{P(Y_{ij}=1\mid U_i)\right\}=\log \frac{P(Y_{ij}=1\mid U_i)}{P(Y_{ij}=0\mid U_i)}={\beta }_0+{\beta }_1{X}_{ij}+U_i$$

• 여기서 $U_i\sim N(0,{\sigma }_u^2)$: 개인별 기저 효과를 나타내는 임의효과 항이다.
• 만약 $X_{ij}$을 구성하는 변수가 한개이고 실험군/대조군 여부$(1/0)$를 알려 준다 하자 (즉, $X_{ij}=1$은 $i$ 번째 개인이 $j$ 번째 관측시점에서 비타민 결핍임을 의미). 그 때의 회귀계수 ${\beta }_1$은 다음과 같이 해석할 수 있다.
• 먼저 $X_{ij}=0$일 때를 보면, $$\exp ({\beta }_0+U_i)=\frac{P(Y_{ij}=1\mid X_{ij}=0,U_i)}{P(Y_{ij}=0\mid X_{ij}=0,U_i)},$$
• 즉 $\exp ({\beta }_0+U_i)$는 $i$ 번째 어린이가 비타민 결핍이 아닐 때, 호흡기 질환의 비감염 확률에 대한 감염 확률의 비(오즈, odds)를 의미한다.
• 한편, $X_{ij}=1$일 때를 보면, $$\exp ({\beta }_0+{\beta }_1+U_i)=\frac{P(Y_{ij}=1\mid X_{ij}=1,U_i)}{P(Y_{ij}=0\mid X_{ij}=1,U_i)},$$
• 는 $i$ 번째 어린이가 비타민 결핍일 때, 호흡기 질환의 비감염 확률에 대한 감염 확률의 비 (오즈)를 보여 준다.
• 두 오즈 (6.2)와 (6.3)의 비(오즈비)는 $$\exp \left({\beta }_1\right)=\frac{\exp ({\beta }_0+{\beta }_1+U_i)}{\exp ({\beta }_0+U_i)}$$
• 즉, 이는 $i$ 번째 어린이의 비타민 결핍 여부에 따른 호흡기 질환에 대한 오즈비를 의미한다.

13.1.4.1 랜덤 절편 모형

• R에서 GLMM을 적용하기 위해 lme4 패키지에 있는 glmer 함수를 적용한다.
• 아래 분석 결과에 의하면, 두 가지 유의한 변수(time, gender)가 호흡기 질환과 관련이 있다.
  • 시간이 지날수록 호흡기 질환 발생 가능성이 증가한다.
  • 여자 어린이의 발생 가능성이 더 낮다.

library(lme4)

Loading required package: Matrix

fit2 = glmer(RESPONSE ~ TIME + VITA + GENDER + AGE + (1|ID), data = ichs, family = binomial)
summary(fit2)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: RESPONSE ~ TIME + VITA + GENDER + AGE + (1 | ID)
   Data: ichs

     AIC      BIC   logLik deviance df.resid 
  1356.7   1388.5   -672.3   1344.7     1494 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1543 -0.3835 -0.1748  0.3054  2.7126 

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 7.422    2.724   
Number of obs: 1500, groups:  ID, 250

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -1.17278    0.53057  -2.210  0.02708 * 
TIME         0.03365    0.01587   2.120  0.03400 * 
VITA         0.63254    0.42912   1.474  0.14047   
GENDER      -1.15364    0.41948  -2.750  0.00596 **
AGE         -0.14792    0.10809  -1.368  0.17116   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr) TIME   VITA   GENDER
TIME   -0.235                     
VITA   -0.306  0.010              
GENDER -0.356 -0.017  0.010       
AGE    -0.711 -0.008 -0.033 -0.024

13.2 랜덤 절편 + 랜덤 기울기 모형

• 아래 분석 결과에 의하면, 랜덤 절편만을 포함한 모형이 더 작은 AIC 를 갖는다. 즉, 연구자는 랜덤 절편 모형을 더 선호할 것이다.

fit22 = glmer(RESPONSE ~ TIME + VITA + GENDER + AGE + (TIME|ID), data = ichs,
             family = binomial)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00295204 (tol = 0.002, component 1)

summary(fit22)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: RESPONSE ~ TIME + VITA + GENDER + AGE + (TIME | ID)
   Data: ichs

     AIC      BIC   logLik deviance df.resid 
  1360.4   1402.9   -672.2   1344.4     1492 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0868 -0.3699 -0.1718  0.3088  2.5999 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 ID     (Intercept) 8.46239  2.90902       
        TIME        0.00125  0.03535  -0.66
Number of obs: 1500, groups:  ID, 250

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -1.22598    0.54799  -2.237  0.02527 * 
TIME         0.04391    0.02566   1.711  0.08712 . 
VITA         0.61214    0.43245   1.415  0.15692   
GENDER      -1.14445    0.42230  -2.710  0.00673 **
AGE         -0.15823    0.10950  -1.445  0.14847   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr) TIME   VITA   GENDER
TIME   -0.327                     
VITA   -0.272 -0.089              
GENDER -0.357  0.064 -0.005       
AGE    -0.645 -0.129 -0.014 -0.035
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00295204 (tol = 0.002, component 1)