Linear Mixed Effect Models With Longitudinal Data

1 Giới Thiệu

1.1 Dữ liệu: Bệnh tự kỉ Autism ở trẻ em

Dữ liệu thực hành cho bài viết này là dataset về bệnh tự kỉ ở trẻ em 2-13 tuổi, được download từ site này: “http://www-personal.umich.edu/~bwest/autism.csv”

Sau khi download về máy tính, các bạn đọc file dữ liệu autism.csv vào R để phân tích.

library(tidyverse)
library(lme4)
library(nlme)
library(lmerTest)
# Đọc dữ liệu
autism=read.csv("F:/R/Data/autism.csv")
# Bảng cấu trúc dữ liệu
head(autism)

##   age vsae sicdegp childid
## 1   2    6       3       1
## 2   3    7       3       1
## 3   5   18       3       1
## 4   9   25       3       1
## 5  13   27       3       1
## 6   2   17       3       3

dim(autism)

## [1] 612   4

Dataset này có 612 dòng và 4 biến số:

• age: tuổi

• childid: mã số bệnh nhi

• sicdegp chẩn đoán bệnh và

• vsae: điểm số đánh giá phát triển hành vi của trẻ vào các thời điểm: tháng thứ 2, thứ 3, thứ 5, thứ 9 và tháng thứ 13 cho mỗi bệnh nhi.

Trong đó chẩn đoán bệnh sicdegp có 3 nhóm:

• Nhóm 1: chẩn đoán autism điển hình
• Nhóm 2: rối loạn phát triển phổ biến (PDR)
• Nhóm 3: không phân loại (nhóm đối chứng)

Ở lần khám đầu tiên (tháng thứ 2), có 116 trẻ thuộc nhóm có bệnh: Nhóm 1 + Nhóm 2 và 40 bệnh nhi thuộc nhóm 3

Vậy đây là dataset có biến số liên tục vsae được đo lặp lại ở 5 thời điểm, gọi là longitudinal data.

1.2 Mục tiêu phân tích

Chúng ta muốn biết:

• Biến đổi của vsae theo thời gian như thế nào
• Điểm số vsae thay đổi khác nhau như thế nào giữa các nhóm trẻ được chẩn đoán các nhóm bệnh khác nhau

2 Mô tả về mẫu nghiên cứu

Phân tích tóm tắt về các biến số:

attach(autism)
sicdegp.f <- factor(sicdegp)
age.f <- factor(age)
# Add the new variables to the data frame object.
autism.updated <- data.frame(autism, sicdegp.f, age.f)
# Cấu trúc bảng dữ liệu
head(autism.updated)

##   age vsae sicdegp childid sicdegp.f age.f
## 1   2    6       3       1         3     2
## 2   3    7       3       1         3     3
## 3   5   18       3       1         3     5
## 4   9   25       3       1         3     9
## 5  13   27       3       1         3    13
## 6   2   17       3       3         3     2

# Summarise data
summary(age.f)

##   2   3   5   9  13 
## 156 150  91 120  95

table(childid)

## childid
##   1   2   3   4   6   8   9  10  12  13  14  15  16  17  18  19  21  22 
##   5   5   5   4   4   3   4   4   3   3   5   5   4   5   2   4   4   3 
##  24  27  28  30  31  32  33  35  36  37  38  39  40  41  42  43  44  45 
##   4   5   4   4   4   3   5   2   3   2   5   4   5   3   4   5   3   4 
##  46  47  48  49  50  51  55  57  58  59  60  61  62  63  64  65  66  67 
##   5   3   4   5   4   3   5   4   4   5   3   2   5   4   4   4   5   4 
##  70  71  76  77  78  80  81  82  85  86  87  88  91  92  95  96  97  99 
##   3   3   3   5   3   4   4   4   2   4   5   2   4   4   5   5   5   4 
## 100 101 104 105 106 107 108 109 110 111 112 113 114 115 116 117 119 120 
##   4   4   4   3   4   4   4   5   4   5   4   5   1   4   4   4   3   3 
## 122 123 124 126 128 129 131 133 134 135 136 139 141 142 143 145 146 147 
##   5   4   3   4   3   2   5   3   4   4   4   5   5   4   4   5   4   4 
## 148 150 151 152 153 154 155 156 158 159 161 165 167 168 169 170 171 172 
##   5   5   5   4   5   3   4   3   4   3   4   3   2   5   5   2   5   4 
## 173 174 175 178 179 180 181 182 183 184 185 186 187 190 191 193 194 195 
##   2   4   5   3   4   4   4   4   5   4   4   4   4   5   5   2   2   3 
## 196 197 198 200 201 202 204 205 207 208 209 210 211 212 
##   4   4   4   2   4   4   4   4   4   2   4   3   1   5

# Số BN tham gia tại mỗi thời điểm nghiên cứu, từ tháng 2 - tháng thứ 13
table(sicdegp.f, age.f)

##          age.f
## sicdegp.f  2  3  5  9 13
##         1 50 48 29 37 28
##         2 66 64 36 48 41
##         3 40 38 26 35 26

attach(autism.updated)
# Tóm tắt điểm số vsae:
summary(vsae)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1.00   10.00   14.00   26.41   27.00  198.00       2

# vsae trung bình ở các thời điểm đo
round(tapply(vsae, age.f, mean, na.rm=TRUE), 2)

##     2     3     5     9    13 
##  9.09 15.26 21.48 39.55 60.60

# Min vsae ở các thời điểm đo
tapply(vsae, age.f, min, na.rm=TRUE)

##  2  3  5  9 13 
##  1  4  4  3  7

# Max vsae ở các thời điểm đo
tapply(vsae, age.f, max, na.rm=TRUE)

##   2   3   5   9  13 
##  20  63  77 171 198

Có 158 bệnh nhi tham gia vào nghiên cứu này. Có 2 bệnh nhi không bắt đầu nghiên cứu tại tháng thứ 2 mà vào nghiên cứu ở những lần khám sau. Lần khám đầu tiên, 156 bệnh nhi được phân thành 3 nhóm như sau:

• Nhóm 1 (autism): 50

• Nhóm 2 (rối loạn phát triển phổ biến): 66

• Nhóm 3 (nhóm không phân loại, nhóm chứng): 40

Số lần khám đánh giá trên mỗi bệnh nhi khác nhau, nhiều nhất 5 lần, ít nhất là 1 lần. Có 2 giá trị vsae bị missing (có tham gia khám nhưng không ghi nhận giá trị của vsae).

2.1 Biến thiên của vsae trong mẫu nghiên cứu

Mô tả khuynh hướng phát triển chung của điểm số vsae theo thời gian trên biểu đồ:

p1 = ggplot(autism.updated, aes(x=age, y=vsae, group=childid)) +
  geom_line(aes(color=childid))+
  geom_point()+
  scale_x_continuous("Age", breaks=seq(0, 14, 1))+
  ggtitle("VSAE By Indivisuals")

p2 = ggplot(autism.updated, aes(x=age, y=vsae)) +
  geom_smooth(method = 'loess')+
  geom_point()+
  scale_x_continuous("Age", breaks=seq(0, 14, 1))+
  ggtitle("Overall VSAE Trend")
  
gridExtra::grid.arrange(p1,p2,ncol=2)

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing missing values (geom_point).

Điểm số vsae tăng dần theo thời gian với biến thiên khác nhau khá lớn và ngày càng lớn hơn tại các thời điểm về sau.

2.2 So sánh biến thiên của vsae trong ba nhóm bệnh nhi

Cũng với dạng biểu đồ trên nhưng được đánh giá và so sánh ba nhóm bệnh nhi:

attach(autism.updated)

## The following objects are masked _by_ .GlobalEnv:
## 
##     age.f, sicdegp.f

## The following objects are masked from autism.updated (pos = 3):
## 
##     age, age.f, childid, sicdegp, sicdegp.f, vsae

## The following objects are masked from autism:
## 
##     age, childid, sicdegp, vsae

ggplot(autism.updated, aes(x=age, y=vsae, group=childid)) +
  geom_line(aes(color=childid))+
  geom_point(aes(color=childid))+
  facet_grid(. ~ sicdegp.f)

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_point).

ggplot(autism.updated, aes(x=age, y=vsae)) +
  geom_smooth(method = 'loess')+
  geom_point()+
  facet_grid(. ~ sicdegp.f)

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing missing values (geom_point).

Nhóm autism điển hình (nhóm 1), vsae tăng chậm nhất, so với nhóm 2 và nhóm chứng (3).

Nhìn vào biểu đồ chúng ta thấy tương quan giữa biến outcome (vsae) và biến tiên đoán như age, trong các nhóm bệnh sicdegp không thể là hồi qui tuyến tính, mà có thể có dạng bậc hai (bình phương).

Dữ liệu đo được ở những thời điểm khác nhau sẽ có tương quan với nhau, vì thế mô hình ANOVA không phù hợp để phân tích.

Mô hình Linear Mixed Effect sẽ được chọn lựa để phân tích trong trường hợp này.

3 Xây dựng mô hình Linear Mixed Effect Models

Nhắc lại cấu trúc mô hình fixed effect model và mixed effect model

3.1 Vì sao gọi là RANDOM EFFECT và FIXED EFFECT ?

Trong trường hợp này, thời điểm đo lường chỉ số vsae trong nghiên cứu là các tháng 2, 3, 5, 9, 13 được chọn ngẫu nhiên ra từ quá trình phát triển của bệnh nhi, không thể đại diện cho toàn quá trình. Vì thế, mục tiêu nghiên cứu là muốn sử dụng dữ liệu từ những thời điểm hạn chế này để suy luận ra toàn bộ quá trình phát triển bệnh lí của trẻ em. Biến số age trong dataset này là nguồn random variance của điểm số vsae, vì thế age được gọi là random variable.

Trong khi đó, biến số nhóm bệnh sicdegp gồm 3 nhóm, là đại diện cho toàn bộ các phân loại chẩn đoán trong bệnh lí này. Vì thế nhóm bệnh (sicdegp) là fixed variable.

Mô hình LMM được xây dựng vừa có random variable (age) và fixed variable (sicdegp) gọi là mô hình hỗn hợp, mixed effect model.

3.2 Chuẩn bị dữ liệu

Nhóm sicdegp = 3 là nhóm chứng. Chúng ta phải chuyển label thành sicdegp = 0 để hàm lmer() nhận ra nhóm chứng này vì mặc định nhóm có thứ tự nhỏ nhất là nhóm chứng (0).

Vì age nhỏ nhất là age = 2, để mô hình chạy tốt như mặc định ta phải tạo ra age2 = age -2, như thế age2 sẽ có giá trị nhỏ nhất là age2 = 0, lúc này intercept sẽ có ý nghĩa.

# Recode the SICDEGP factor for model fitting.
age2 <- age - 2
sicdegp2 <- sicdegp
sicdegp2[sicdegp == 3] <- 0
sicdegp2[sicdegp == 2] <- 2
sicdegp2[sicdegp == 1] <- 1
sicdegp2.f <- factor(sicdegp2)
table(sicdegp2.f )

## sicdegp2.f
##   0   1   2 
## 165 192 255

3.3 Xây dựng mô hình lmer

Dự đoán là mô hình này có biến thiên bậc hai, chúng ta đưa age2^2 vào mô hình để kiểm chứng.

Chú ý: với (age2 + I(age2^2) | childid), chúng ta muốn mô hình đánh giá hai biến số age2 và age2^2 có liên quan có ý nghĩa với “childid”, tức là các nhóm bệnh nhi hay không.

Nói cách khác, chỉ số vsae thay đổi theo thời gian nhưng với phương cách khác nhau giữa các nhóm bệnh nhi. Hai nhóm bệnh sẽ được so sánh với nhóm chứng.

model1.lmr <- lmer(vsae ~ age2 + I(age2^2) + sicdegp2.f +
                            age2*sicdegp2.f + I(age2^2)*sicdegp2.f + 
                            (age2 + I(age2^2) | childid),
                          REML = T, data = autism.updated)

summary(model1.lmr)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## vsae ~ age2 + I(age2^2) + sicdegp2.f + age2 * sicdegp2.f + I(age2^2) *  
##     sicdegp2.f + (age2 + I(age2^2) | childid)
##    Data: autism.updated
## 
## REML criterion at convergence: 4608
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3015 -0.3905 -0.0392  0.2837  6.8667 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  childid  (Intercept)  1.4332  1.1972              
##           age2        14.2647  3.7769    0.06      
##           I(age2^2)    0.1562  0.3952    0.89 -0.40
##  Residual             37.6414  6.1353              
## Number of obs: 610, groups:  childid, 158
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            13.78643    0.82251 294.32000  16.762  < 2e-16 ***
## age2                    5.59928    0.78666 136.00000   7.118 5.73e-11 ***
## I(age2^2)               0.20417    0.08537 110.30000   2.392  0.01847 *  
## sicdegp2.f1            -5.43483    1.11108 296.81000  -4.891 1.64e-06 ***
## sicdegp2.f2            -4.03533    1.04597 296.13000  -3.858  0.00014 ***
## age2:sicdegp2.f1       -3.28366    1.08196 141.93000  -3.035  0.00286 ** 
## age2:sicdegp2.f2       -2.75204    1.01795 139.42000  -2.704  0.00771 ** 
## I(age2^2):sicdegp2.f1  -0.13881    0.11833 113.76000  -1.173  0.24322    
## I(age2^2):sicdegp2.f2  -0.13396    0.11023 111.85000  -1.215  0.22682    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) age2   I(g2^2) scd2.1 scd2.2 a2:2.1 a2:2.2 I(2^2):2.1
## age2        -0.400                                                      
## I(age2^2)    0.447 -0.623                                               
## sicdegp2.f1 -0.740  0.296 -0.331                                        
## sicdegp2.f2 -0.786  0.314 -0.351   0.582                                
## ag2:scdg2.1  0.291 -0.727  0.453  -0.407 -0.229                         
## ag2:scdg2.2  0.309 -0.773  0.482  -0.229 -0.407  0.562                  
## I(g2^2):2.1 -0.322  0.450 -0.721   0.447  0.254 -0.629 -0.347           
## I(g2^2):2.2 -0.346  0.483 -0.774   0.256  0.447 -0.351 -0.627  0.559

AIC(model1.lmr)

## [1] 4639.978

BIC(model1.lmr)

## [1] 4710.594

Giải thích kết quả:

Intercept = 13.78643, p < 0.001 có ý nghĩa thống kê,

age2 = 5.59928 suy ra age = 7.59928: mỗi năm, điểm số vsae tăng trung bình 7.59928 điểm

I(age2^2) = 0.20417: mỗi năm, điểm số vsae tăng trung bình 0.20417 điểm cho age2^2

sicdegp2.f1 = -5.43483: Nhóm bệnh nhi autism điển hình sẽ kém hơn 5.43483 điểm só với nhóm chứng

sicdegp2.f2 = -4.03533: Nhóm bệnh nhi rối loạn phát triển phổ biến sẽ kém hơn 4.03533 điểm só với nhóm chứng

age2 có liên quan có ý nghĩa thống kê với nhóm bệnh. Trong đó nhóm autism giảm điểm nhiều hơn nhóm rối loạn phát triển phổ biến khi được so sánh với nhóm chứng, tương ứng là -3.28366 và -2.75204.

3.4 Mô hình lmer() chỉ dành cho random intercepts, không có random slops

Khác nhau ở cú pháp: (1 | childid)

model2.lmr <- lmer(vsae ~ age2 + I(age2^2) + sicdegp2.f +
                            age2*sicdegp2.f + I(age2^2)*sicdegp2.f + 
                            (1 | childid),
                          REML = T, data = autism.updated)

summary(model2.lmr)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## vsae ~ age2 + I(age2^2) + sicdegp2.f + age2 * sicdegp2.f + I(age2^2) *  
##     sicdegp2.f + (1 | childid)
##    Data: autism.updated
## 
## REML criterion at convergence: 5448
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5464 -0.4111 -0.0128  0.2655  5.5564 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  childid  (Intercept) 158.6    12.59   
##  Residual             348.0    18.65   
## Number of obs: 610, groups:  childid, 158
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            12.86222    3.11329 343.20000   4.131 4.53e-05 ***
## age2                    7.60779    1.39133 459.30000   5.468 7.49e-08 ***
## I(age2^2)              -0.05433    0.12853 464.10000  -0.423  0.67270    
## sicdegp2.f1            -4.52171    4.19387 342.10000  -1.078  0.28172    
## sicdegp2.f2            -2.45744    3.94985 342.90000  -0.622  0.53425    
## age2:sicdegp2.f1       -4.92155    1.92243 462.20000  -2.560  0.01078 *  
## age2:sicdegp2.f2       -5.92929    1.80760 463.00000  -3.280  0.00112 ** 
## I(age2^2):sicdegp2.f1   0.04776    0.17808 466.70000   0.268  0.78866    
## I(age2^2):sicdegp2.f2   0.28964    0.16672 466.60000   1.737  0.08299 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) age2   I(g2^2) scd2.1 scd2.2 a2:2.1 a2:2.2 I(2^2):2.1
## age2        -0.512                                                      
## I(age2^2)    0.402 -0.963                                               
## sicdegp2.f1 -0.742  0.380 -0.299                                        
## sicdegp2.f2 -0.788  0.403 -0.317   0.585                                
## ag2:scdg2.1  0.370 -0.724  0.697  -0.507 -0.292                         
## ag2:scdg2.2  0.394 -0.770  0.741  -0.292 -0.506  0.557                  
## I(g2^2):2.1 -0.290  0.695 -0.722   0.399  0.229 -0.963 -0.535           
## I(g2^2):2.2 -0.310  0.742 -0.771   0.230  0.400 -0.537 -0.964  0.556

AIC(model2.lmr)

## [1] 5470.017

BIC(model2.lmr)

## [1] 5518.565

3.5 So sánh hai mô hình lmer:

anova(model1.lmr,model2.lmr)

## Data: autism.updated
## Models:
## model2.lmr: vsae ~ age2 + I(age2^2) + sicdegp2.f + age2 * sicdegp2.f + I(age2^2) * 
## model2.lmr:     sicdegp2.f + (1 | childid)
## model1.lmr: vsae ~ age2 + I(age2^2) + sicdegp2.f + age2 * sicdegp2.f + I(age2^2) * 
## model1.lmr:     sicdegp2.f + (age2 + I(age2^2) | childid)
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## model2.lmr 11 5472.1 5520.6 -2725.0   5450.1                             
## model1.lmr 16 4635.1 4705.7 -2301.6   4603.1 846.96      5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p < 0.001, như vậy, việc đưa random slops vào mô hình sẽ tăng tính hiệu quả của mô hình, bằng chứng là các chỉ số AIC, BIC, Loglik đều giảm.

4 So sánh fixed model vs. mixed models:

Xây dựng mô hình glm() chỉ có fixed effect:

# Loại NA ra khỏi dữ liệu để chạy glm()
autism.updated <- subset(data.frame(autism, sicdegp2.f, age2),!is.na(vsae))
# Mo hình glm                           
model2.glm <- gls(vsae ~ age2 +  I(age2^2) + sicdegp2.f, data = autism.updated)
summary(model2.glm)

## Generalized least squares fit by REML
##   Model: vsae ~ age2 + I(age2^2) + sicdegp2.f 
##   Data: autism.updated 
##        AIC      BIC   logLik
##   5581.419 5607.851 -2784.71
## 
## Coefficients:
##                  Value Std.Error   t-value p-value
## (Intercept)  24.184714 2.2165469 10.910987  0.0000
## age2          3.499495 0.9140148  3.828707  0.0001
## I(age2^2)     0.093775 0.0841288  1.114665  0.2654
## sicdegp2.f1 -22.240835 2.4939930 -8.917762  0.0000
## sicdegp2.f2 -16.397184 2.3417777 -7.002024  0.0000
## 
##  Correlation: 
##             (Intr) age2   I(2^2) scd2.1
## age2        -0.477                     
## I(age2^2)    0.380 -0.965              
## sicdegp2.f1 -0.618  0.021 -0.015       
## sicdegp2.f2 -0.658  0.031 -0.028  0.570
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.13544492 -0.47685476 -0.03361707  0.24422291  5.99198611 
## 
## Residual standard error: 23.42649 
## Degrees of freedom: 610 total; 605 residual

So sánh hai mô hình:

• lmer model 1: AIC = 4639.978 ; BIC = 4710.594; residual: 37.6414

• glm model: AIC = 5581.419 ; BIC: 5607.851; residual: 605

Qua đó cho thấy linear mixed effect model chính xác hơn trong ước tính biến số vsae trong trường hợp này.

Kết luận: Mô hình linear Mixed Effect kết hợp fixed và random effects phù hợp với dữ liệu longitudinal để đánh giá biến thiên của biến số liên tục (vsae) trong mối tương quan với các biến số categorical và continuous khác.

Cám ơn các bạn đã theo dỏi

---