Sử dụng package Gamlss

Bài 4: Thành phần bất định (Random effects)

1 Giới thiệu về Mixed model và Random Effects

Chúng ta đã đi qua 3 chặng đường trên hành trình khám phá package GAMLSS, các bạn đã nhận ra gamlss là một công cụ dựng mô hình hồi quy độc đáo, không giống bất cứ package nào khác trong R. Cho đến thời điểm này, chúng ta chỉ mới xét đến thành phần xác định (fixed effect) trong mô hình Gamlss. Trong chặng thứ 4 này, Nhi sẽ giới thiệu với các bạn về thành phần bất định hay ngẫu nhiên (random terms hay random effect) trong một mô hình gọi là hỗn hợp (mixed model).

Loại mô hình này còn được biết dưới nhiều tên gọi khác nhau trong sách vở, như multilevel model, hierarchical model, mixed model, nested model, random effect model, subject specific model, variance component model, random intercept và random slope model, repeated measure ANOVA… (nghe có vẻ đáng sợ phải không các bạn ?). Sự “phong phú” về thuật ngữ này đã khoác lên chủ đề này một màu sắc bí hiểm. Mixed model còn là thử thách với những người không chuyên vì bạn không thể học sâu về nó mà không sử dụng code lập trình (SAS hoặc R).

Điều tồi tệ hơn nữa đó là trong Y văn, tác giả rất hiếm khi mô tả đích danh Mixed model, và gần như không bao giờ trình bày nội dung mô hình mà chỉ nói một cách chung chung là ANOVA hoặc linear model (dù mixed model và random effect là chìa khóa cho hầu hết thử nghiệm lâm sàng, longitudinal study, split-plot và repeated measure design). Ngay cả khi họ dùng mixed model, thì cũng chỉ trình bày fixed effect mà giấu đi phần random trong mô hình. Do đó, không có nhiều sinh viên Y khoa và bác sĩ biết dùng Mixed model và random effect.

Một người học hồi quy tuyến tính chưa thể đạt đến cảnh giới tối cao nếu chưa thực hành qua Mixed model, và như Nhi đã nói, đây là phương pháp quan trọng, nguồn gốc của những phân tích ANOVA trong thử nghiệm lâm sàng.

Tóm lại, GLMM rất quan trọng và cần thiết nếu bạn muốn dùng hồi quy để giải quyết các nghiên cứu lâm sàng, thí nghiệm. GLMM phức tạp, nhưng Nhi sẽ cố gắng trình bày về Mixed model một cách đơn giản nhất có thể, và do chúng ta đang sử dụng R, khoảng cách đến thành công không còn xa bao nhiêu.

Mixed model là một chủ đề rất rộng, có thể trình bày về nó trong cả quyển sách. Nhi giới thiệu một vài tựa sách mà bạn có thể tìm đọc như Mixed Models, Theory and Application with R của Eugene Demidenko, Generalized, linear and Mixed model của Charles E. McCulloch, Applied Mixed models in Medicine của Helen Brown và Robin Prescott, Introduction to Mixed Modelling – Beyond Regression and Analysis of Variance, Analysis of Generalized Linear Mixed Models, Mixed effects models and extensions in Ecology with R của Alain F.Ziir và cộng sự, vân vân… Nhi không có khả năng đi quá sâu vào chi tiết vào chủ đề này chỉ trong 1 bài viết ngắn. Mục tiêu của bài thực hành này do đó chỉ nhắm đến việc giới thiệu một số thuật ngữ, khái niệm chính, và cú pháp của một số hàm trong Gamlss. Để đơn giản và nhất quán, Nhi sẽ sử dụng Generalized linear mixed model (GLMM) như tên gọi của phương pháp.

Một cách cơ bản, mô hình hồi quy tuyến tính cho phép chúng ta ước lượng (và kiểm soát) một phần sự biến thiên của biến kết quả, ta gọi đây là phần xác định hay hiệu ứng xác định (fixed effect). Tuy nhiên, như ta biết, không có mô hình nào hoàn hảo: ngay cả thành phần xác định của một mô hình phù hợp nhất vẫn không cho phép giải thích toàn bộ sự biến thiên của biến kết quả, phần sai số còn dư thường được cho là hệ quả của tính ngẫu nhiên (trong khi chọn mẫu chẳng hạn). Tuy nhiên các bạn có thể xem thành tố ngẫu nhiên này như 1 biến số (yếu tố) và đưa nó vào trong mô hình (dưới dạng xác suất điều kiện), như một cố gắng để giải thích nhiều hơn nữa phần sai số còn dư. Như vậy, hiệu ứng ngẫu nhiên (random effect) là những biến số mang tính bất định, ngẫu nhiên được gắn thêm vào mô hình nguyên thủy ở nhiều vị trí khác nhau (thí dụ random intercept, random slope của 1 predictor), cho phép chúng ta ước lượng được sự biến thiên ngẫu nhiên bên ngoài những quy luật mà mô hình fixed effect đã mô tả.

Một thí dụ đơn giản: thử nghiệm lâm sàng về hiệu quả điều trị của loại thuốc X trên n bệnh nhân, fixed effect ở đây là phân nhóm điều trị (1 nhóm Placebo vs nhóm điều trị), còn random effect có thể xem như là cơ địa của bệnh nhân.

Một trong những dấu hiệu nhận dạng một bài toán cần tới mixed model và random effect, đó là tồn tại sự phân nhóm (cụm) trong cấu trúc dữ liệu. Sự phân cụm này có thể là do chủ ý của thiết kế thí nghiệm hoặc hoàn toàn ngẫu nhiên. Việc nhận dạng cấu trúc phân nhóm này nhiều lúc là do suy luận và giả định của chúng ta. Trong thí dụ nêu trên, ta có thể giả định mỗi bệnh nhân không hẳn là 1 cá thể, nhưng là 1 mẫu đại diện cho 1 nhóm trong quần thể thực tế.

Một số loại random effect đặc biệt

1)Nested effect

Một yếu tố (F) lồng trong yếu tố khác (G), hiệu ứng của G bị giới hạn bên trong một cấp bậc nhất định của G, thí dụ G1 nhưng không hiện diện ở cấp bậc khác như G2, G3…. Thí dụ minh họa, một nghiên cứu khảo sát sự thay đổi thị lực dưới tác động của một điều trị phẫu thuật nhãn khoa. Biến số kết quả được đánh giá riêng cho mắt bên trái và mắt bên phải (yếu tố Side) cho mỗi bệnh nhân, trước và sau khi phẫu thuật, như vậy hiệu ứng ngẫu nhiên của yếu tố “Side” gồm 2 levels: Trái/phải được lồng trong yếu tố “Id của Bệnh nhân”, do đó mô hình được viết dưới dạng:

Outcome ~ Treatment + (1|Id) + (1|Id:Side )

Crossed effect:

Đây là hiệu ứng ngẫu nhiên của các yếu tố giao nhau, hay không có cơ sở để lập luận rằng chúng lồng ghép với nhau. Thí dụ để khảo sát mức độ tương hợp giữa 3 bác sĩ chẩn đoán hình ảnh, mỗi người trong số họ sẽ lần lượt đọc 3 hình ảnh Ct scan của cùng một bệnh nhân được chụp trong 3 lần khác nhau bằng 3 thiết bị khác nhau. Như vậy yếu tố ngẫu nhiên sẽ là crossed effect = (1|Doctor) + (1|Replication) chứ không phải là nested effect = (1|Doctor:Replication) vì mỗi level của replication là như nhau cho cả 3 levels của Doctor.

Multilevel model hay Hierarchical model:

Các thuật ngữ này phổ biến trong khoa học xã hội và kinh tế hơn là Y học, thí dụ thường gặp nhất trong sách vở là nghiên cứu về các chỉ số giáo dục với đơn vị là học sinh với yếu tố nẫu nhiên được phân cấp theo thành phố, trường, lớp … tuy nhiên khái niệm về cấu trúc đa cấp trong mô hình cũng có thể gặp trong nghiên cứu Y khoa và dịch tễ học. Thí dụ cộng đồng dân cư có thể phân cấp thành Thành phố, Quận, phường, tổ, v..v ,random effect có dạng cấu trúc đa cấp cũng tương tự như hiệu ứng lồng ghép (nested) và có thể được viết như:

Điểm thi đại học ~ Giới tính + (1|Thành phố) + (1|Thành phố:Trường) + (1|Thành phố:Trường:Lớp)

Repeated measure và Longitudinal study

Đây là thiết kế nghiên cứu rất phổ biến, khi một kết quả (outcome) được khảo sát nhiều lần (tại nhiều thời điểm khác nhau) trên cùng một đối tượng. Thí dụ: Một thử nghiệm thuốc trong đó 1 biomarker hay đại lượng lâm sàng được khảo sát tại thời điểm t0, sau đó bệnh nhân được phân nhóm điều trị (placebo/thuốc X) và ta lặp lại phép đo mỗi 3 tháng trong vòng 1 năm, như vậy ta có 5 thời điểm Baseline,3,6,9,12. Mô hình với random slope do đó sẽ được viết là: Outcome ~ Baseline*Treatment + Time+ (Treatment |Patient Id)

2 Thí dụ minh họa: Dataset Epileptic

Epileptic là 1 dataset kinh điển cho mô hình random effect, nó đã được sử dụng bở nhiều tác giả như Thall ,Vail [1990], Breslow và Clayton[1993], Lee và Nelder [1996, 2000]. Đây là một nghiên cứu thử nghiệm lâm sàng trong đó có 2 nhóm bệnh nhân (n=59) mắc bệnh động kinh. Tần suất cơn động kinh được đếm trong mỗi 2 tuần và lặp lại 5 lần. Sau 2 tuần đầu tiên, việc điều trị bắt đầu với Placebo và loại thuốc điều trị cần khảo sát. Sau đó số cơn động kinh được ghi nhận sau mỗi 2 tuần trong 4 lần tái khám.

library(tidyverse)

library(gamlss)

my_theme <- function(base_size =8, base_family = "sans"){
  theme_bw(base_size = base_size, base_family = base_family) +
    theme(
      panel.grid.major = element_line(color = "gray"),
      panel.grid.minor = element_blank(),
      panel.background = element_rect(fill = NA),
      strip.background = element_rect(fill = "#510166", color = NA, size =0.5),
      strip.text = element_text(face = "bold", size = 8, color = "white"),
      legend.position = "bottom",
      legend.justification = "center",
      legend.background = element_blank(),
      legend.margin = margin(0.5,0.5,0.5,0.5)
    )
}

3 Thăm dò số liệu

Như thường lệ, Nhi sẽ bắt đầu thăm dò số liệu bằng biểu đồ

Ghi chú: kết quả ở thời điểm ban đầu trước điều trị được chuyển sang thang logarithm (lbase) Tuổi của bệnh nhân cũng được chuyển sang thang logarithm

epil=read.csv("https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/MASS/epil.csv")

epil$base=as.factor(epil$base)

epil%>%ggplot(aes(x=y,fill=trt))+
  geom_density(alpha=.8)+
  geom_histogram(aes(y=..density..,fill=trt),colour="black",binwidth = 3,alpha=0.7,show.legend = F)+
  my_theme()+scale_fill_manual(values=c("#ff9c07","#d60a29"))+
  facet_wrap(~period,ncol=2)+
  ggtitle("Density + Histogram")

epil%>%ggplot(aes(x=trt,y=y,group=factor(lbase)))+
  geom_smooth(aes(color=lbase),size=1.2,alpha=0.2,se=F,method="lm")+
  my_theme()+
  scale_colour_gradient(low="#ffa528",high="#db1154")+
  facet_wrap(~period,scales="free")

epil%>%ggplot(aes(x=trt,y=y,group=factor(lage)))+
  geom_smooth(aes(color=lage),size=1,alpha=0.5,se=F,method="lm")+
  my_theme()+
  scale_colour_gradient(low="#ffa528",high="#db1154")+
  facet_wrap(~period,scales="free")

epil%>%ggplot(aes(x=period,y=y,group=1))+
  geom_path(aes(color=trt),alpha=0.8,size=1)+
  geom_point(aes(color=trt,shape=trt),alpha=0.5,size=2)+
  my_theme(5)+
  scale_color_manual(values=c("#ff9c07","#d60a29"))+
  facet_wrap(~subject,scales="free",ncol=7)

epil%>%ggplot(aes(x=period,y=y,group=1))+
  geom_path(aes(color=trt),alpha=0.4,size=1)+
  geom_point(aes(color=trt,shape=trt),alpha=0.5,size=2)+
  my_theme(5)+
  scale_color_manual(values=c("#ff9c07","#d60a29"))

epil%>%ggplot(aes(x=period,y=y,group=trt,color=trt))+
  geom_smooth(alpha=0.5,se=F,method="glm")+
  my_theme()+
  scale_color_manual(values=c("#ff9c07","#d60a29"))

epil%>%ggplot(aes(x=period,y=log(y),group=factor(lbase),color=lbase))+
  geom_smooth(alpha=0.5,se=F,method="glm")+
  my_theme()+
  scale_colour_gradient2(low="#f7620c",mid="#d60652",high="#62017a")

epil%>%ggplot(aes(x=period,y=log(y),group=lage,color=lage))+
  geom_smooth(alpha=0.5,se=F,method="glm")+
  scale_colour_gradient2(high="#aa03ba",low="#fc6e02",mid="#d60652")+
  my_theme()

epil%>%ggplot(aes(x=period,y=log(y),group=lbase,color=lbase))+
  geom_smooth(alpha=0.5,se=F,method="glm")+
  my_theme()+facet_wrap(~trt,ncol=2)+
  scale_colour_gradient2(high="#aa03ba",low="#fc6e02",mid="#d60652")

epil%>%ggplot(aes(x=lage,y=log(y),group=trt))+
  geom_smooth(aes(color=trt,fill=trt),alpha=0.3,method="lm")+
  theme_bw()+
  scale_color_manual(values=c("#ff9c07","#d60a29"))+
  scale_fill_manual(values=c("#f7b80c","#d60a29"))

epil%>%ggplot(aes(x=lage,y=log(y),group=trt,color=trt))+
  geom_point(alpha=0.5)+
  geom_smooth(alpha=0.5,se=F,method="glm")+
  my_theme()+facet_wrap(~period,scales="free")+
  scale_color_manual(values=c("#ff9c07","#d60a29"))

epil%>%ggplot(aes(x=lbase,y=log(y),group=trt,color=trt))+
  geom_point(alpha=0.5)+
  geom_smooth(alpha=0.5,se=F,method="glm")+
  my_theme()+facet_wrap(~period,scales="free")+
  scale_color_manual(values=c("#ff9c07","#d60a29"))

Các biểu đồ gợi ý rằng :

Biến kết quả là một số đếm, với dấu hiệu overdispersion nhẹ, nên có thể họ phân phối phù hợp cho nó sẽ là negative binomial

Giá trị baseline có tác động lên giá trị biến kết quả trong quá trình điều trị tại cả 4 thời điểm, nhưng chưa rõ liệu có tương tác giữa nó và hiệu ứng điều trị hay không

Nhìn chung, nhóm điều trị giảm rõ nét tần suất cơn động kinh so với nhóm Placebo, nhưng hiệu ứng này không hằng định ở mọi bệnh nhân

Tuổi bệnh nhân không có ảnh hưởng đáng kể đến kết quả điều trị

Nhi chuyển dạng biến số period thành một contrast variable với 4 cấp bậc giá trị:

epil$visit=(2*epil$period-5)/10

factor(epil$visit)%>%levels(.)

## [1] "-0.3" "-0.1" "0.1"  "0.3"

4 Hàm gamlssNP()

Hàm gamlssNP() có công dụng là mô hình hóa marginal likelihood, bản chất của nó là ước tính tổng quát (hiệu ứng ngẫu nhiên được phân tích trong bản thân algorithm của Gamlss).

Tùy lựa chọn của người dùng, hàm gamlssNP cho phép áp dụng 2 phương pháp:

1. Phép cầu phương Gaussian (Quadrature) nhằm ước lượng marginal likelihhod bằng cách thay thế tích phân trong khoảng phân phối chuẩn của random effect bằng phép cộng diện tích của nhiều tứ giác (con số do người dùng quy ước).

Công dụng của phương pháp này là dựng mô hình với Random intercept (1 yếu tố ngẫu nhiên liên tục và có phân phối chuẩn) cho tham số vị trí trung tâm Mu.

2. phương pháp phi tham số (non parametric) cho yếu tố ngẫu nhiên rời rạc ; với công dụng là mô hình hóa random intercept hoặc cả random intercept và random slope cho từ 1 đến 4 tham số của 1 họ phân phối (Mu, Sigma, Nu và Tau).

Cú pháp hàm gamlssNP

Hàm gamlssNP được hỗ trợ bởi module gamlss.mx có công dụng:

Mô hình hóa random intercept liên tục, phân phối chuẩn, cho Mu Hay: Dựng mô hình phi tham số với random intercept và random slope cho cả Mu, Sigma, Nu và Tau của 1 họ phân phối tùy chọn.

Tùy chỉnh random = có nội dung là 1 công thức để biểu thị cho random effect, thí dụ random intercept ở cấp độ phân nhóm: random = ~ 1|F

Random slope ở cấp độ phân nhóm: random = ~ X|F

Nếu random effect ở cấp độ quan sát: không cần khai báo cho tùy chỉnh random Tùy chỉnh Mixture : phươn pháp ước lượng:

Phép cầu phương Gaussian: gq, dùng cho random intercept phân phối chuẩn Np= nonparametric random effect models K= số điểm trong phép ước lượng np hay gp Khi muốn mở rộng random effect cho các tham số khác như Sigma, Nu và Tau, ta phải khai báo công thức cho các thành tố này, phần random effect được biểu thị bằng tên 1 matrix gọi là MASS, thí dụ: sigma.fo ~ X+ MASS cho biết MASS nằm ở vị trí Intercept, còn sigma.fo ~ X* MASS cho biết MASS nằm cả ở intercept và slope.

library(gamlss.mx)

epil$subject=as.factor(epil$subject)


m0<-gamlssNP(y~1, 
               random=~1|subject,
               K=15,
               mixture="gq",
               data=epil, 
               family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..14 ..15 ..16 ..
## 17 ..18 ..19 ..
## EM algorithm met convergence criteria at iteration   19

## Global deviance trend plotted.

m00<-gamlssNP(y~trt, 
              random=~1|subject,
              K=15,
              mixture="gq",
              data=epil, 
              family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..14 ..15 ..16 ..
## 17 ..18 ..19 ..20 ..21 ..22 ..23 ..24 ..25 ..26 ..27 ..28 ..29 ..
## EM algorithm met convergence criteria at iteration   29

## Global deviance trend plotted.

m10<-gamlssNP(y~trt+visit, 
             random=~1|subject,
             K=15,
             mixture="gq",
             data=epil, 
             family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..14 ..15 ..16 ..
## 17 ..18 ..19 ..20 ..21 ..22 ..23 ..24 ..25 ..26 ..27 ..28 ..29 ..
## EM algorithm met convergence criteria at iteration   29

## Global deviance trend plotted.

m11<-gamlssNP(y~trt+lbase+visit, 
             random=~1|subject,
             K=15,
             mixture="gq",
             data=epil, 
             family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..14 ..
## EM algorithm met convergence criteria at iteration   14

## Global deviance trend plotted.

m12<-gamlssNP(y~trt*lbase+visit, 
             random=~1|subject,
             K=15,
             mixture="gq",
             data=epil, 
             family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..
## EM algorithm met convergence criteria at iteration   13

## Global deviance trend plotted.

m13<-gamlssNP(y~trt*lbase+visit+lage, 
              random=~1|subject,
              K=15,
              mixture="gq",
              data=epil, 
              family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..
## EM algorithm met convergence criteria at iteration   12

## Global deviance trend plotted.

GAIC(m0,m00,m10,m11,m12,m13,k=log(nrow(epil)))

##     df      AIC
## m11  6 1284.977
## m12  7 1289.034
## m13  8 1292.597
## m0   3 1329.775
## m00  4 1331.816
## m10  5 1334.508

LR.test(m11,m12)

##  Likelihood Ratio Test for nested GAMLSS models. 
##  (No check whether the models are nested is performed). 
##  
##        Null model: deviance= 1252.194 with  6 deg. of freedom 
##  Altenative model: deviance= 1250.788 with  7 deg. of freedom 
##  
##  LRT = 1.406041 with 1 deg. of freedom and p-value= 0.2357148

LR.test(m10,m11)

##  Likelihood Ratio Test for nested GAMLSS models. 
##  (No check whether the models are nested is performed). 
##  
##        Null model: deviance= 1307.189 with  5 deg. of freedom 
##  Altenative model: deviance= 1252.194 with  6 deg. of freedom 
##  
##  LRT = 54.995 with 1 deg. of freedom and p-value= 1.207923e-13

summary(m11)

## Warning in summary.gamlss(m11): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  "NBI Mixture with NO" 
## 
## Call:  
## gamlssNP(formula = y ~ trt + lbase + visit, random = ~1 | subject,  
##     family = NBI, data = epil, K = 15, mixture = "gq") 
## 
## Fitting method: RS()

## Warning in summary.gamlss(m11): observations with zero weight not used for
## calculating dispersion

## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.81816    0.05329  34.119  < 2e-16 ***
## trtprogabide -0.32936    0.07393  -4.455 8.64e-06 ***
## lbase         1.01523    0.04878  20.813  < 2e-16 ***
## visit        -0.27334    0.16426  -1.664   0.0962 .  
## z             0.49906    0.03785  13.184  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.0228     0.1753  -11.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3540 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  230 
##                       at cycle:  1 
##  
## Global Deviance:     1252.194 
##             AIC:     1264.194 
##             SBC:     1284.977 
## ******************************************************************

#method np

m20<-gamlssNP(y~trt+lbase+visit, 
              random=~1|subject,
              K=15,
              mixture="np",
              data=epil, 
              family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..14 ..15 ..16 ..
## 17 ..18 ..19 ..20 ..21 ..22 ..23 ..24 ..25 ..26 ..27 ..28 ..29 ..30 ..31 ..32 ..33 ..
## 34 ..35 ..36 ..37 ..38 ..39 ..40 ..41 ..
## EM algorithm met convergence criteria at iteration   41

## Global deviance trend plotted.

## EM Trajectories plotted.

m21<-gamlssNP(y~trt*lbase+visit, 
              random=~1|subject,
              K=15,
              mixture="np",
              data=epil, 
              family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..14 ..15 ..16 ..
## 17 ..18 ..19 ..20 ..21 ..22 ..23 ..24 ..25 ..26 ..27 ..28 ..29 ..30 ..31 ..32 ..33 ..
## 34 ..35 ..36 ..37 ..38 ..39 ..40 ..41 ..42 ..43 ..44 ..45 ..46 ..47 ..48 ..49 ..50 ..
## 51 ..52 ..53 ..54 ..
## EM algorithm met convergence criteria at iteration   54

## Global deviance trend plotted.

## EM Trajectories plotted.

m22<-gamlssNP(y~trt*lbase+visit, 
             random=~trt|subject,
             K=20,
             mixture="np",
             data=epil, 
             family=NBI)

## 1 ..2 ..3 ..4 ..5 ..6 ..7 ..8 ..9 ..10 ..11 ..12 ..13 ..14 ..15 ..16 ..
## 17 ..18 ..19 ..20 ..21 ..22 ..23 ..24 ..25 ..26 ..27 ..28 ..29 ..30 ..31 ..32 ..33 ..
## 34 ..35 ..36 ..37 ..38 ..39 ..40 ..41 ..42 ..43 ..44 ..45 ..46 ..47 ..48 ..49 ..50 ..
## 51 ..52 ..53 ..54 ..55 ..56 ..57 ..58 ..59 ..60 ..
## EM algorithm met convergence criteria at iteration   60

## Global deviance trend plotted.

## EM Trajectories plotted.

GAIC(m20,m21,m22,k=log(nrow(epil)))

##     df      AIC
## m20 33 1422.983
## m21 34 1428.355
## m22 63 1580.233

plotMP(m22$mass.points[,1], m22$mass.points[,2], m22$prob,theta=50,
       phi=30,col = "gold")

summary(m20)

## Warning in summary.gamlss(m20): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("NBI Mixture with NP", "Negative Binomial type I Mixture with NP" ) 
## 
## Call:  
## gamlssNP(formula = y ~ trt + lbase + visit, random = ~1 | subject,  
##     family = NBI, data = epil, K = 15, mixture = "np") 
## 
## Fitting method: RS()

## Warning in summary.gamlss(m20): observations with zero weight not used for
## calculating dispersion

## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.079e-01  2.325e+02   0.003 0.997227    
## trtprogabide -2.805e-01  7.307e-02  -3.838 0.000126 ***
## lbase         9.259e-01  4.838e-02  19.140  < 2e-16 ***
## visit        -2.733e-01  1.605e-01  -1.703 0.088661 .  
## MASS2         5.383e-05  2.327e+02   0.000 1.000000    
## MASS3         6.408e-05  2.325e+02   0.000 1.000000    
## MASS4         7.265e-05  2.325e+02   0.000 1.000000    
## MASS5         9.060e-05  2.325e+02   0.000 1.000000    
## MASS6         9.373e-01  2.325e+02   0.004 0.996783    
## MASS7         9.577e-01  2.325e+02   0.004 0.996713    
## MASS8         9.601e-01  2.325e+02   0.004 0.996705    
## MASS9         1.456e+00  2.325e+02   0.006 0.995002    
## MASS10        2.015e+00  2.325e+02   0.009 0.993086    
## MASS11        2.015e+00  2.325e+02   0.009 0.993086    
## MASS12        2.015e+00  2.325e+02   0.009 0.993086    
## MASS13        2.015e+00  2.325e+02   0.009 0.993086    
## MASS14        2.015e+00  2.331e+02   0.009 0.993105    
## MASS15        2.015e+00  5.749e+02   0.004 0.997204    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.1224     0.1929     -11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3540 
## Degrees of Freedom for the fit:  33
##       Residual Deg. of Freedom:  203 
##                       at cycle:  1 
##  
## Global Deviance:     1242.676 
##             AIC:     1308.676 
##             SBC:     1422.983 
## ******************************************************************

5 Hàm random() và re()

Khác với hàm gamlssNP(), hàm random() và re( ) sử dụng phép ước lượng likelihood cục bộ (local) còn gọi là phương pháp penalized quasi likelihood (PQL) để ước tính joint likelihood, lúc này sự biến thiên biến số kết quả sẽ được mô hình hóa thông qua xác suất có điều kiện tùy thuộc vào random effect.

Hàm random() trong gamlss là phiên bản của hàm cùng tên trong package gam của Tishibrani được cải biên bởi Rigby và Stasinopoulos để dùng cho package gamlss. Nó cho phép mô hình hóa random intercept có phân phối chuẩn.

Hàm re() là một cánh cửa liên thông với hàm lme() trong package nlme, một công cụ nổi tiếng cho mixed model (Pinheiro và Bates, 2000). Nó sẽ gọi package nlme() và kết hợp package này và gamlss để dựng mixed model với 3 dạng : random intercept và random slope, multilevel modelling và repeated measurement.

#Random function

m30<-gamlss(y~trt+lbase+visit+random(subject), 
           data=epil,
           family=NBI, 
           trace=FALSE)

m31<-gamlss(y~trt*lbase+visit+random(subject), 
            sigma.fo=~random(subject),
            data=epil,
            family=NBI, 
            trace=FALSE)

m32<-gamlss(y~trt*lbase+visit+lage+random(subject), 
            sigma.fo=~random(subject),
            data=epil,
            family=NBI, 
            trace=FALSE)

GAIC(m30,m31,m32,k=log(nrow(epil)))

##           df      AIC
## m30 48.25053 1390.694
## m31 65.63829 1424.132
## m32 65.52979 1425.443

summary(m30)

## ******************************************************************
## Family:  c("NBI", "Negative Binomial type I") 
## 
## Call:  gamlss(formula = y ~ trt + lbase + visit + random(subject),  
##     family = NBI, data = epil, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.82477    0.04801  38.007  < 2e-16 ***
## trtprogabide -0.32357    0.06576  -4.920 1.88e-06 ***
## lbase         0.99605    0.04322  23.045  < 2e-16 ***
## visit        -0.27075    0.14744  -1.836   0.0679 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.4855     0.2426  -10.24   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  236 
## Degrees of Freedom for the fit:  48.25053
##       Residual Deg. of Freedom:  187.7495 
##                       at cycle:  10 
##  
## Global Deviance:     1127.062 
##             AIC:     1223.563 
##             SBC:     1390.694 
## ******************************************************************

rdf=data.frame(mu=as.vector(m30$mu.s),
               muL=as.vector(m30$mu.s-(sqrt(m30$mu.var)/sqrt(236))),
               muU=as.vector(m30$mu.s+(sqrt(m30$mu.var)/sqrt(236))),
               Id=epil$subject)

rdf%>%mutate(.,Positive=ifelse(rdf$mu>0,"Positive","Negative"))%>%
  ggplot(aes(x=as.factor(Id),color=Positive))+
    geom_errorbar(aes(ymin=muL,ymax=muU))+
  geom_point(aes(y=mu))+
  geom_hline(yintercept = 0,color="blue",linetype=2)+
  theme_bw(5)+
  coord_flip()+
  scale_x_discrete("Idx")+
  scale_color_manual(values=c("#ff9c07","#d60a29"))

Hàm re() mang lại cho gamlss khả năng tương đương với package nlme(), tức là có thể mô hình hóa những random effect phức tạp mà không cần xét đến giả định phân phối chuẩn. Cách sử dụng hàm re cũn tương tự như package nlme, nếu bạn đã quen thuộc với mixed model thì sẽ thấy nó rất dễ dàng. Hàm re có thể áp dụng cho bất kì tham số phân phối nào, nó có những tùy chỉnh như sau;

fixed = khai báo nội dung cho công thức của fixed effect, nhưng không cần khai báo biến kết quả

random = công thức cho random effects

x1 + … xn |f1/… / fm

Với f1/…/fm là nested structure (khối), khi random effect nằm ở cấp độ từng quan sát, phần sau không cần

x1 + … xn |f

Một danh sách công thức:

List(f1=pdDiag(~ x1), f2=pdSymm(~ x2)) = random intercept và slope cho x1 ở cấp độ factor f1, với covariance structure pdDiag

pdSymm= giá trị mặc định, cấu trúc matrix xác định dương và đối xứng

pdBlocked: cấu trúc dạng khối – phân giác

random=list(f1=pdBlocked(list(pdIndent(_1),pdIdent(f2-1))))

pdCompSymm: cấu trúc đối xứng tổ hợp, thích hợp cho nested random effect (thí dụ factor f2 nested trong factor f1:

random=list(f1=pdCopmSymm(~f2-1))

pdDiag: cấu trúc matrix phân giác: giả định tính độc lập giữa random effects random=pdDiag(~x1+x2)

trong đó random slope cho x1 và x2 và intercept giả định như độc lập với nhau pdIdent: cấu trúc matrix Identity

khi cấp bậc của factor f1 có variance như nhau: random=pdIdent(~f1-1)

Correlation: xác định cấu trúc tương quan bên trong mỗi khối (nhóm) hay error variation

Đây là nơi chúng ta có thể đưa vào những cấu trúc tương quan đặc biệt của time series hay spatial structures

Về time series ta có:

corCompSymm: cấu trúc đối xứng tổ hợp, tương ứng với tương quan đồng nhất (uniform)

corAR1: autoregressive process cấp 1, hàm correlation = rho^k với k=0,1,2…

corCAR1: continuous autoregressive process AR1 process cho 1 hiệp biến số thời gian liên tục, rho^s với s > =0

corARMA= autoregressive moving average process, với arbitrary orders cho thành phần autoregressive và moving average (ARMA(p,q)) tương ứng với sự kết hợp giữa p autoregressive parameters và q moving average parameters

corExp: exponential spatial correlation

corGaus: Gaussian spatial correlation

corLin: linear spatial correlation;

corRatio: rational quadratic spatial correlation

corSpher: spherical spatial correlation

method =”ML” hay “REML”

#RE function

m40<-gamlss(y~trt+lbase+visit+
            re(random=~1|subject),
            family=NBI,
            data=epil,
            trace=FALSE)

m41<-gamlss(y~trt*lbase+visit+
              re(random=~1|subject),
            family=NBI,
            data=epil,
            trace=FALSE)

m42<-gamlss(y~trt*lbase+visit+
              re(random=~1|subject),
            sigma.fo=~re(random=~1|subject),
            family=NBI,
            data=epil,
            trace=FALSE)

GAIC(m40,m41,m42,k=log(nrow(epil)))

##           df      AIC
## m40 48.25031 1390.693
## m41 48.75377 1393.872
## m42 65.58082 1423.796

rdf=data.frame(mu=as.vector(m42$mu.s),
               sigma=as.vector(m42$sigma.s),
               Id=epil$subject)

rdf%>%gather(mu:sigma,key="Parameter",value="Effect")%>%
  mutate(.,Positive=ifelse(Effect>0,"Positive","Negative"))%>%
  ggplot(aes(x=as.factor(Id),y=Effect,color=Positive))+
  geom_point()+
  geom_hline(yintercept = 0,color="blue",linetype=2)+
  theme_bw(5)+
  coord_flip()+
  scale_x_discrete("Subjects")+
  facet_wrap(~Parameter,ncol=2)+
  scale_color_manual(values=c("#ff9c07","#d60a29"))

6 Marginal effects

term.plot(m40, pages=1)

term.plot(m30, pages=1)

epil%>%mutate(.,predict=predict(m30,type="response"))%>%
  ggplot(aes(x=trt,y=log(predict)))+
  geom_boxplot(aes(fill=trt),alpha=0.7)+
  my_theme(5)+facet_wrap(~period,ncol=1)+coord_flip()+
  scale_fill_manual(values=c("#ff9c07","#d60a29"))

7 Kết luận

Nhận xét: Nhi vừa trình bày một cách khái quát công dụng của 4 hàm trong gamlss cho phép mô hình hóa random effect hay thành phần ngẫu nhiên, bất định. Có thể thấy tuy gamlss đã hấp thu nội dung từ nhiều package có trước nó nhưng về mọi phương diện, quy trình dựng mixed model trong gamlss không tương thích với những công cụ kinh điển như nlme hay lme4. Tức là ta không thể trích xuất kết quả và suy diễn thống kê dễ dàng như output của package lme4 hay nlme. Đây là nhược điểm của gamlss.

Vậy khi nào ta nên dùng gamlss: chỉ khi bạn gặp vấn đề khó về phân phối (vì gamlss cho phép mô hình hóa nhiều họ phân phối hơn bất cứ package nào khác). Thực ra Nhi luôn dùng lme và nlme cho mixed model nhưng chỉ mới dùng đến gamlss 2 lần duy nhất, một lần là cho 1 longitudinal study cho 1 em bác sĩ nội trú bên VN vì biến số kết quả của em ấy có phân phối Gamma; một lần khác là cho một nghiên cứu lâm sàng với 10 biến số kết quả là count data có phân phối zero inflated negative binomial.

Xin cảm ơn và hẹn gặp lại.