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

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

주 출처

보조 출처

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

13.1.1 개요

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

13.1.2 모형식과 가정

  • 각 개체의 반응값을 \({\bf y}_{i}=\left(y_{i 1}, \cdots, y_{i m_{i}}\right)\)라 하자. (\(i=1, \ldots, N\))
  • Recall: linear mixed effect model (LMM) \[ y_{ij}= {\bf x}_{ij}^T \boldsymbol{\beta} + {\bf z}_{ij}^T {\bf u}_{i} + \epsilon_{ij}, \quad j=1,\ldots, m_i, \quad i=1, \cdots, N \]
  • 여기서 고정효과 관련 공변량 \({\bf x}_{ij}\), 고정효과 계수 \(\boldsymbol{\beta}\), 랜덤효과 관련 공변량 \({\bf z}_{ij}\)은 non-random으로 가정된다. 랜덤효과 \({\bf u}_{i}\)는 확률변수로 가정된다. 만약 \(E(y_{ij}|{\bf u}_i) = \mu_{ij}\)로 쓰면, LMM은 다음과 같이 다시 쓸 수 있었다. \[ \mu_{ij}= {\bf x}_{ij}^T \boldsymbol{\beta} + {\bf z}_{ij}^T {\bf u}_{i}, \quad j=1,\ldots, m_i, \quad i=1, \cdots, N \]
  • GLMM의 모형식: known link function \(g: \mathbb{R} \rightarrow \mathbb{R}\)에 대하여, \[ g(\mu_{ij}) = g(E(y_{ij}|{\bf u}_i)) = {\bf x}_{ij}^T \boldsymbol{\beta} + {\bf z}_{ij}^T {\bf u}_{i}, \quad j=1,\ldots, m_i, \quad i=1, \cdots, N, \]
  • GLMM의 모형식 (벡터 표현): \[ g(E({\bf y}_i|{\bf u}_i) ) = {\bf X}_i \boldsymbol{\beta} + {\bf Z}_i {\bf u}_{i}, \quad i=1, \cdots, N \]
  • (위 식에서 \(g(\cdot)\)는 elementwise하게 적용된다고 가정하였다)
  • 확률 가정 \[ \begin{gather} E({\bf u}_i) = {\bf 0}, \qquad {\rm Var}({\bf u}_i) = {\bf R}_i, \nonumber \\ {\bf u}_i ~~ {\rm iid} \nonumber \\ y_{i1}, \ldots, y_{im_i} \text{~are independent given~} {\bf u}_i, \nonumber \end{gather} \]

13.1.3

  • 이항 경시적 자료(binary categorical data): \[ \begin{gather} y_{ij} \mid {\bf u}_i \sim B(1, \mu_{ij}); \nonumber \\ (\text{i.e.},~ P(y_{ij}=1 \mid {\bf u}_i) =: \mu_{ij},~ P(y_{ij}=0 \mid {\bf u}_i) =: 1-\mu_{ij}) \nonumber \\ \text{where }~~{\rm logit}(\mu_{ij}) := \log \left( \frac{\mu_{ij}}{1 - \mu_{ij}} \right) = {\bf x}_{ij}^T \boldsymbol{\beta} + {\bf z}_{ij}^T {\bf u}_{i}, \quad j=1,\ldots, m_i, \quad i=1, \cdots, N. \nonumber \end{gather} \]
  • 계수형 경시적 자료(count data): \[ \begin{gather} y_{ij} \mid {\bf u}_i \sim {\rm Poisson}(\mu_{ij}) \nonumber \\ \log(\mu_{ij}) = {\bf x}_{ij}^T \boldsymbol{\beta} + {\bf z}_{ij}^T {\bf u}_{i}, \quad j=1,\ldots, m_i, \quad i=1, \cdots, N. \nonumber \end{gather} \]

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\left(Y_{i j}=1 \mid U_{i}\right) \right\} = \log \frac{P\left(Y_{i j}=1 \mid U_{i}\right)}{P\left(Y_{i j}=0 \mid U_{i}\right)}=\beta_{0}+\beta_{1} X_{i j}+U_{i} \]

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