Mô hình liên hợp

Phần 1: Trường phái RE Maximum Likelihood

1 Dẫn nhập

Đây là một bài thực hành với tham vọng cao, đó là giới thiệu đến các bạn sinh viên Y khoa và bác sĩ về mô hình liên hợp. Đây là một công cụ tinh tế, phức tạp và hữu ích với nhiều ứng dụng lâm sàng; tuy nhiên nó gần như chỉ lưu hành nội bộ trong giới thống kê chuyên nghiệp (cho các hãng dược phẩm chẳng hạn), chứ không được phổ biến cho người thực hành lâm sàng.

Phương pháp này được phát triển khoảng đầu thập niên 1990,và đặc biệt dành riêng cho ngành Y khoa (tất cả dataset thực hành đều là clinical trial), vì nó cho phép kiểm chứng những giả thuyết phức tạp, khai thác thông tin hữu ích từ các nghiên cứu thử nghiệm thuốc.

Cho đến năm 2015, lý thuyết về Joint model đã hoàn thiện và có thể tiếp cận trong R bằng tới 3 cách : trường phái Bayes, ước lượng hợp lý cực đại (RE Maximum likelihood) và Hierarchical likelihood.

Trong bài này, Nhi sẽ giới thiệu trước cách tiếp cận dùng REML, cách làm Bayes sẽ được hướng dẫn vào 1 dịp khác. Trong bài Nhi chủ yếu dùng package JM của tác giả Dimitris Rizopoulos (2012-2016).

Mục tiêu của bài thực hành gồm có:

1. Thăm dò dữ liệu của một Longitudinal và Time-Event data
2. Dựng Joint model từ Cox-PH model và Mixed linear model
3. Diễn giải kết quả mô hình
4. Kiểm tra phẩm chất mô hình
5. Ứng dụng mô hình để tiên lượng cho cá thể bệnh nhân

Do tính chất phức tạp của vấn đề, Nhi khuyến khích các bạn tìm đọc thêm tài liệu lý thuyết (bài này chỉ nhằm để thả thính…). Cụ thể là quyển sách :

2 Dataset AIDS

Trước hết ta tải 1 số package,dataset aids nằm trong package JM

library(JM) # package dựng Joint model
library(tidyverse) # data wrangling và ggplot
library(magrittr) # Toán tử pipe
library(ggkm) #Thăm dò Time event data

AIDS dataset có nguồn gốc từ một nghiên cứu (NC) thử nghiệm dược phẩm lâm sàng năm 1994 trên 467 bệnh nhân nhiễm HIV. Mục tiêu của NC này đó là so sánh hiệu quả điều trị giữa 2 loại thuốc antiretroviral khác nhau tên là Didanosine (mã = ddI) và zalcitabine (ddC). Outcome chính của NC này là biến cố tử vong với thời gian theo dõi kéo dài 2 năm. Tuy nhiên, NC còn có 1 outcome thứ hai là số lượng bạch cầu CD4 (gọi tắt là CD4). CD4 được khảo sát ở thời điểm đầu tiên, sau đó mỗi bệnh nhân được sắp xếp ngẫu nhiên vào 1 trong 2 phân nhóm trị liệu và CD4 tiếp tục được khảo sát định kỳ tại tháng thứ 2,6,12 và 18.

Vào lúc NC kết thúc, 188 bệnh nhân đã tử vong và người ta có tổng cộng 1405 giá trị CD4.

Như ta thấy, đây là một nghiên cứu có thiết kế phức tạp với 2 outcomes lần lượt là : 1 sự kiện (tử vong) và 1 biến định lượng (CD4) được theo dõi theo thời gian. Như vậy, mỗi bộ phận này tương ứng với một thiết kế nghiên cứu khác nhau : Longitudinal (cho CD4) và Time-Event (Survival) analysis.

Dĩ nhiên chúng ta có thể xử lý riêng từng bộ phận bằng 1 mô hình hồi quy hỗn hợp (Mixed model) (cho CD4) và 1 mô hình hồi quy Cox-PH (cho nguy cơ tử vong). Tuy nhiên, vấn đề không đơn giản như ta nghĩ. Để biết ta phải đương đầu với những khó khăn nào, Nhi sẽ thăm dò dữ liệu.

3 Thăm dò dữ liệu

Tổng số bệnh nhân là 467

aids$patient %>%
  unique() %>%
  length()

## [1] 467

Xem cấu trúc dữ liệu, thí dụ cho 1 bệnh nhân N°2

aids%>%filter(.$patient==2)%>%knitr::kable()

patient	Time	death	CD4	obstime	drug	gender	prevOI	AZT	start	stop	event
2	19	0	6.324555	0	ddI	male	noAIDS	intolerance	0	6	0
2	19	0	8.124038	6	ddI	male	noAIDS	intolerance	6	12	0
2	19	0	4.582576	12	ddI	male	noAIDS	intolerance	12	18	0
2	19	0	5.000000	18	ddI	male	noAIDS	intolerance	18	19	0

Thời gian theo dõi trung bình là 13.9 tháng, thấp nhất là khoảng 2 tuần, cao nhất là 21 tháng. Điều này chứng tỏ là tỉ lệ tử vong và mất mát dữ liệu khá cao (theo dự tính thời gian theo dõi tối đa 18 tháng.

summary(aids$Time)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.47   12.23   14.07   13.89   17.00   21.40

Phần 1: Diễn tiến CD4 theo thời gian (Thành tố Longitudinal)

Biểu đồ sau đây trình bày diễn tiến thay đổi của CD4 (ghi chú: đã được khai căn bậc 2) ở mỗi cá thể bệnh nhân và khuynh hướng trung bình giữa 2 phân nhóm điều trị, và giữa 2 nhóm có và không có tử vong.

Biểu đồ này cho chúng ta thấy 3 vấn đề:

1. Đầu tiên ta có thể thấy đó là không phải bất cứ bệnh nhân nào cũng có đủ 5 giá trị CD4. Vậy,có vấn đề thiếu sót dữ liệu và đây là vấn đề chung cho tất cả nghiên cứu lâm sàng. Nguyên nhân của sót dữ liệu có rất nhiều : bệnh nhân hủy bỏ cuộc hẹn, bệnh nhân rời bỏ nghiên cứu giữa chừng, mẫu bệnh phẩm bị thất lạc, thiết bị đo hư hỏng đúng vào lịch hẹn, mất hồ sơ, nhập liệu sai, và dĩ nhiên – bệnh nhân tử vong trước khi nghiên cứu kết thúc (Nhóm tử vong có số time point thấp , đa phần chấm dứt ở tháng thứ 6, là điểm hẹn thứ 3)

2. Một chi tiết khác, khó nhận ra hơn nhưng có thể tưởng tượng, đó là sai lệch ngẫu nhiên về thời điểm lấy mẫu CD4. Có thể bệnh nhân không đến đúng hẹn như kế hoạch, bác sĩ và họ có thể đã đặt hẹn ngẫu nhiên với sai lệch vài ngày đến 1-2 tuần trong cùng tháng, đơn vị mốc thời gian 2,6,12,18 do đó không chính xác tuyệt đối.

3. Vấn đề thứ 3 : Mỗi bệnh nhân là một « phân nhóm » đặc biệt trong nghiên cứu. Lộ trình diễn tiến của CD4 khác nhau giữa BN này và BN khác. Đó là nguyên nhân tại sao random effect và mixed model sẽ được sử dụng để xét đến yếu tố ngẫu nhiên cá thể.

library(colorRamps)
library(RColorBrewer) 

mypalette = length(unique(aids$patient))
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))

aids%>%mutate(.,Death=ifelse(.$death==1,"Death","Survived"))%>%
  ggplot()+
  geom_path(aes(x=obstime,y=CD4,color=patient),show.legend = F,alpha=0.5)+
  geom_smooth(aes(x=as.numeric(obstime),y=CD4,fill=drug,linetype=drug),
              color="black",
              method="lm",formula = y ~ splines::bs(x,3),alpha=0.4)+
  scale_color_manual(values = getPalette(mypalette))+
  facet_wrap(~Death)+
  theme_bw()

Có sự khác biệt rõ về khuynh hướng trung bình của CD4 giữa bệnh nhân sống sót và tử vong, cũng như giữa 2 phân nhóm điều trị, tuy nhiên theo thời gian CD4 giảm dần;

aids%>%filter(.$Time==.$stop)%>%
  mutate(.,Death=ifelse(.$death==1,"Death","Censored"))%>%
  ggplot()+
  geom_boxplot(aes(x=Death,y=CD4,fill=Death),alpha=0.7)+
  coord_flip()+
  facet_wrap(~drug,scales="free",ncol=1)+
  theme_bw()

aids%>%filter(.$Time==.$stop)%>%
  mutate(.,Death=ifelse(.$death==1,"Death","Censored"))%>%
  ggplot()+
  geom_boxplot(aes(x=drug,y=CD4,fill=drug))+
  coord_flip()+
  facet_wrap(~Death,scales="free",ncol=1)+
  theme_bw()

aids%>%
  ggplot(aes(x=drug,y=Time,fill=drug))+
  geom_boxplot(color="black")+
  geom_violin(alpha=0.3)+
  coord_flip()+
  theme_bw()

Nếu muốn giải quyết độc lập biến số CD4 theo thời gian, ta có thể dùng mô hình Mixed linear model, bao gồm mô hình Random intercept và random slope như sau:

mypalette = length(unique(aids$patient))
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))

rint=lme(CD4 ~ obstime, 
         random = ~ 1 | patient,
         data = aids)

aids%>%mutate(.,Fitted=rint$fitted[,1]+rint$fitted[,2])%>%
  ggplot()+
  geom_smooth(aes(x=as.numeric(obstime),y=Fitted,color=patient),
              method="lm",show.legend = F,se=F)+
  geom_smooth(aes(x=as.numeric(obstime),y=Fitted),
              method="lm",formula = y ~ x,show.legend = F,se=F,color="red")+
  scale_color_grey()+
  theme_bw()

rslope=lme(CD4 ~ obstime, 
           random = ~ obstime | patient,
                  data = aids)

aids%>%mutate(.,Fitted=rslope$fitted[,1]+rslope$fitted[,2])%>%
  ggplot()+
  geom_smooth(aes(x=as.numeric(obstime),y=Fitted,color=patient),
              method="lm",show.legend = F,se=F)+
  geom_smooth(aes(x=as.numeric(obstime),y=Fitted),
              method="lm",formula = y ~ x,show.legend = F,se=F,color="red")+
  scale_color_grey()+
  theme_bw()

rslope=lme(CD4 ~ obstime, 
           random = ~ obstime | patient,
           data = aids)

Liên hệ giữa diễn tiến CD4 và biến cố tử vong

Có sự liên hệ giữa CD4 và nguy cơ tử vong tại bất cứ thời điểm nào. Đây là một giả thuyết logic, vì lượng tế bào CD4 biểu thị cho khả năng miễn dịch của bệnh nhân. Suy giảm miễn dịch sẽ dẫn đến nhiễm trùng cơ hội và tăng nguy cơ tử vong.

binomial_smooth <- function(...) {
  geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}

aids%>%mutate(.,Death=ifelse(.$death==1,"Death","Survived"))%>%
  ggplot(aes(x=CD4,y=death,fill=drug,linetype=drug))+
  binomial_smooth(color="black",alpha=0.3)+
  facet_wrap(~obstime,ncol=2)+
  theme_bw()

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

mypalette = length(unique(aids$drug))
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))

aids%>%mutate(.,Death=ifelse(.$death==1,"Death","Survived"))%>%
  ggplot()+
  geom_boxplot(aes(x=as.factor(obstime),y=CD4,fill=drug),alpha=0.7)+
  scale_fill_manual(values = getPalette(mypalette))+
  facet_wrap(~Death,ncol=1)+
  coord_flip()+
  theme_bw()

aids%>%mutate(.,Death=ifelse(.$death==1,"Death","Survived"))%>%
  ggplot()+
  geom_histogram(aes(x=CD4,y=..density..,fill=Death,bins=5,color=Death),alpha=0.5)+
  scale_fill_manual(values = getPalette(mypalette))+
  scale_color_manual(values = getPalette(mypalette))+
  facet_grid(obstime~drug,scales = "free")+
  theme_bw()

## Warning: Ignoring unknown aesthetics: bins

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Có sự liên hệ giữa CD4 và thời gian sống của bệnh nhân

mypalette = length(unique(aids$drug))
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))

aids%>%
  ggplot(aes(x=CD4,y=Time))+
  geom_point(alpha=0.5,aes(color=drug))+
  geom_smooth(aes(x=CD4,y=Time,fill=drug,linetype=drug),
              color="black",
              method="lm",formula = y ~ splines::bs(x,3),alpha=0.4)+
  scale_color_manual(values = getPalette(mypalette))+
  theme_bw()

mypalette = 12
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))

aids%>%
  ggplot(aes(x=CD4,y=Time))+
  geom_boxplot(aes(group=cut_width(CD4,2.5),
                   fill=cut_width(CD4,2.5),
               color=cut_width(CD4,2.5)),
               alpha=0.5,
               show.legend = F)+
  geom_smooth(aes(x=CD4,y=Time,linetype=drug),
              color="black",
              method="lm",formula = y ~ splines::bs(x,3),alpha=0.4)+
  scale_fill_manual(values = getPalette(mypalette))+
  scale_color_manual(values = getPalette(mypalette))+
  theme_bw()

Tác động của điều trị lên nguy cơ tử vong

Sau đây là 1 vài biểu đồ trình bày 3 hàm Survival, Nguy cơ tích lũy và CDF, cho thấy thay đổi của nguy cơ tử vong theo thời gian.

ggplot(aids, 
       aes(time = Time, status = death, fill = drug, color = drug))+ 
  geom_km()+ 
  geom_kmband()+
  theme_bw()

ggplot(aids, 
       aes(time = Time, status = death, fill = drug, color = drug))+ 
  geom_km(trans = "cumhaz")+ 
  geom_kmband(trans = "cumhaz")+
  theme_bw()+ ylab("Cumulative hazard")

ggplot(aids, 
       aes(time = Time, status = death, fill = drug, color = drug))+ 
  geom_km(trans = "event")+ 
  geom_kmband(trans = "event")+
  theme_bw()+ ylab("Cumulative Events (CDF)")

library(survminer)

## Loading required package: ggpubr

## Warning: package 'ggpubr' was built under R version 3.4.1

library(survival)

fit <- survfit(Surv(Time, death) ~ drug, data = aids)

ggsurvplot(fit, 
           data = aids,
           size=1,
           conf.int = TRUE,
           pval = TRUE,
           risk.table = TRUE, 
           risk.table.col = "strata",
           risk.table.height = 0.3,
           ggtheme = theme_bw(),
           risk.table.y.text.col = T, 
           risk.table.y.text = FALSE
           )

fit <- survfit(Surv(Time, death) ~ drug, data = aids)

Một mô hình Cox đơn thuần không đủ khả năng để phát hiện sự khác biệt về hiệu ứng giữa 2 loại thuốc điều trị

coxphmod0<- coxph(Surv(Time, death) ~ drug, data = aids)

summary(coxphmod0)

## Call:
## coxph(formula = Surv(Time, death) ~ drug, data = aids)
## 
##   n= 1405, number of events= 412 
## 
##            coef exp(coef) se(coef)     z Pr(>|z|)
## drugddI 0.13022   1.13908  0.09859 1.321    0.187
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## drugddI     1.139     0.8779    0.9389     1.382
## 
## Concordance= 0.528  (se = 0.013 )
## Rsquare= 0.001   (max possible= 0.983 )
## Likelihood ratio test= 1.75  on 1 df,   p=0.1865
## Wald test            = 1.74  on 1 df,   p=0.1866
## Score (logrank) test = 1.75  on 1 df,   p=0.1863

ggcoxdiagnostics(coxphmod0)

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.064417

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.13087

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.5233e-030

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.017127

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.064417

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.13087

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 5.5233e-030

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 0.017127

ggforest(coxphmod0)

## Warning in .get_data(model, data = data): The `data` argument is not
## provided. Data will be extracted from model fit.

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

Khi biến số CD4 được đưa vào mô hình Cox, lúc này mọi thứ có vẻ khả quan hơn chút ít: Mô hình cho thấy CD4 thật sự có liên hệ với nguy cơ tử vong. sqrt(CD4) giảm 1 đơn vị làm tăng nguy cơ tử vong thêm 1.18 lần.

Tuy nhiên, mô hình này chưa tốt và độ phù hợp dữ liệu kém. Việc dùng trực tiếp CD4 như 1 hiệp biến số có thật sự đúng đắn chưa ? Câu trả lời là chưa ! Tuy CD4 có liên hệ với nguy cơ tử vong, nhưng liên hệ giữa chúng không thể được khảo sát hoàn toàn bởi mô hình Cox, vì mô hình Cox chỉ cho phép khảo sát 1 hiệp biến số độc lập, nằm bên ngoài thực thể bệnh lý (thí dụ phân nhóm điều trị, ô nhiễm môi trường, nhiễm trùng…) trong khi CD4 là một hiệp biến số nội tại, nó nằm ngay trong căn bệnh.

coxphmod<- coxph(Surv(Time, death) ~ drug+CD4, data = aids)

summary(coxphmod)

## Call:
## coxph(formula = Surv(Time, death) ~ drug + CD4, data = aids)
## 
##   n= 1405, number of events= 412 
## 
##             coef exp(coef) se(coef)       z Pr(>|z|)    
## drugddI  0.21579   1.24085  0.09881   2.184    0.029 *  
## CD4     -0.16292   0.84966  0.01463 -11.133   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## drugddI    1.2408     0.8059    1.0224    1.5060
## CD4        0.8497     1.1769    0.8256    0.8744
## 
## Concordance= 0.683  (se = 0.015 )
## Rsquare= 0.114   (max possible= 0.983 )
## Likelihood ratio test= 169.6  on 2 df,   p=0
## Wald test            = 126  on 2 df,   p=0
## Score (logrank) test = 138.6  on 2 df,   p=0

ggcoxdiagnostics(coxphmod)

ggforest(coxphmod)

## Warning in .get_data(model, data = data): The `data` argument is not
## provided. Data will be extracted from model fit.

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

4 Mô hình liên hợp: tóm tắt

Mô hình liên hợp (Joint model) là một phương pháp thống kê cho phép hợp nhất hai bộ phận : sự biến thiên phụ thuộc thời gian của một biến định lượng (longitudinal response ; thí dụ marker trong 1 bệnh lý) và nguy cơ xảy ra một sự kiện theo thời gian (survival analysis), nhằm khảo sát mối liên hệ giữa chúng.

Một cách đơn giản tối đa, nguyên tắc của mô hình liên hợp có thể được mô tả như sau :

1 mô hình liên hợp có 2 bộ phận bên trong nó (như 2 bánh xe răng cưa) : Một là mô hình phân tích sự kiện (Cox-PH model), và 1 mô hình phân tích phương sai phụ thuộc thời gian với yếu tố ngẫu nhiên (Mixed model).

1. Bộ phận Survival có mục tiêu đo lường mối liên hệ giữa diễn tiến của 1 biến số định lượng (marker M) và nguy cơ xảy ra 1 biến cố (E) theo thời gian. Lưu ý 1 chi tiết quan trọng, đó là giá trị của M tại thời điểm t là 1 giá trị Hư mà Thực (có thực, nhưng không thể hoặc chưa thể xác định được). Điều này mới nghe có vẻ kì cục, nhưng có thể hiểu được, thí dụ nếu thời điểm t rơi vào khoảng giữa 2 cuộc hẹn trong lịch, thì ta không thể nào biết marker khi đó có giá trị bao nhiêu, nhưng trong cơ thể của bệnh nhân khi đó vẫn diễn ra sự thay đổi. Do đó M được khảo sát như 1 bệnh sử tính đến thời điểm t.

Bộ phận survival có vai trò chủ chốt trong mô hình, vì nó cho phép xác định tham số đo lường hiệu ứng của longitudinal outcome M lên nguy cơ xảy ra E.

Trong phần survival còn có hàm nguy cơ gốc h0(t), có thể xem như cơ số cho hàm EXP đi sau. Hàm h0 trong mô hình Cox cũng có nhưng không được xác định.Trong Joint model, h0 phải được xác định, có thể hình dung là chia nhỏ khoảng thời gian khảo sát ra nhiều đoạn thường bằng Piece wise, B-spline, cubic splines … Việc làm này cũng cần phải cẩn trọng vì có nguy cơ gây overfitting. Do đó, có quy tắc đặt ra là giới hạn tổng số biến (bao gồm bậc của hàm spline và biến độc lập trong mô hình) từ 1/10 đến 1/20 tổng số sự kiện trong mẫu.

2. Bộ phận thứ 2 là mô hình Longitudinal. Nhiệm vụ của nó là ước lượng càng chính xác càng tốt lịch sử diễn tiến của biến M (marker). Đây là móc nối giữa 2 mô hình bộ phận.

M được ước tính dựa trên : yếu tố thời gian, matrix mô hình, hiệu ứng cố định (beta) và hiệu ứng ngẫu nhiên (cá thể), và sai số ngẫu nhiên (giả định có phân phối chuẩn với trung bình =0 và phương sai sigma^2).

Mô hình longitudinal có bản chất là 1 linear mixed model có chứa biến số thời gian. Sau này, các tác giả còn mở rộng phần Longitidinal bằng cách đưa thêm vào các thành phần bất định khác.

Ta thử nhìn lại 1 lần nữa 2 quá trình song song: nguy cơ tử vong tích lũy theo thời gian, và CD4 thay đổi theo thời gian:

lp=aids%>%mutate(.,Death=ifelse(.$death==1,"Death","Survived"))%>%
  ggplot()+
  geom_smooth(aes(x=as.numeric(Time),y=CD4,fill=drug,linetype=drug),
              color="black",alpha=0.4)+
  scale_color_manual(values = getPalette(mypalette))+
  theme_bw()

hz=ggplot(aids, 
       aes(time = Time, status = death, fill = drug, color = drug))+ 
  geom_km(trans = "cumhaz")+ 
  geom_kmband(trans = "cumhaz")+
  theme_bw()+ ylab("Cumulative hazard")

gridExtra::grid.arrange(hz,lp,ncol=1)

## `geom_smooth()` using method = 'loess'

5 Dựng Joint model bằng package JM

Đầu tiên, ta dựng 2 mô hình Longitudinal và Cox-PH trong package JM:

# Tạo dataset cho mô hình Cox

aids%>%.[!duplicated(.$patient),]->aids.id

#mixed longitudinal submodel

longcomp=lme(CD4 ~ obstime+ obstime:drug,
             random = ~ obstime | patient, 
             data = aids)

longcomp%>%summary()

## Linear mixed-effects model fit by REML
##  Data: aids 
##        AIC      BIC    logLik
##   7147.575 7184.295 -3566.788
## 
## Random effects:
##  Formula: ~obstime | patient
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 4.5901582 (Intr)
## obstime     0.1738082 -0.155
## Residual    1.7497905       
## 
## Fixed effects: CD4 ~ obstime + obstime:drug 
##                     Value  Std.Error  DF  t-value p-value
## (Intercept)      7.188833 0.22215874 936 32.35899  0.0000
## obstime         -0.163451 0.02080804 936 -7.85519  0.0000
## obstime:drugddI  0.028272 0.02970929 936  0.95163  0.3415
##  Correlation: 
##                 (Intr) obstim
## obstime         -0.160       
## obstime:drugddI  0.000 -0.682
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.32530054 -0.41785802 -0.04720642  0.40631129  4.32727623 
## 
## Number of Observations: 1405
## Number of Groups: 467

# Time event submodel (CoxPH)

coxcomp=coxph(Surv(Time, death) ~ drug, 
              data = aids.id, 
              x = TRUE)

coxcomp%>%summary()

## Call:
## coxph(formula = Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
## 
##   n= 467, number of events= 188 
## 
##           coef exp(coef) se(coef)     z Pr(>|z|)
## drugddI 0.2102    1.2339   0.1462 1.437    0.151
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## drugddI     1.234     0.8104    0.9264     1.643
## 
## Concordance= 0.531  (se = 0.019 )
## Rsquare= 0.004   (max possible= 0.99 )
## Likelihood ratio test= 2.07  on 1 df,   p=0.15
## Wald test            = 2.07  on 1 df,   p=0.1506
## Score (logrank) test = 2.07  on 1 df,   p=0.1498

Hàm jointModel sẽ hợp nhất 2 mô hình bộ phận lại với nhau, với hàm nguy cơ gốc h0 được ước tính bằng phương pháp Spline, Proportional hazard và GH.Quy trình này sẽ mất khá nhiều thời gian, bạn nên kiên nhẫn chờ đến khi mô hình converged.

Package JM tương thích với nhiều hàm trong base R, thí dụ: summary( ) để xem nội dung model.

jointmod=jointModel(lmeObject=longcomp, 
                    survObject=coxcomp,
                         timeVar = "obstime",
                         method = "spline-PH-GH")

summary(jointmod)

## 
## Call:
## jointModel(lmeObject = longcomp, survObject = coxcomp, timeVar = "obstime", 
##     method = "spline-PH-GH")
## 
## Data Descriptives:
## Longitudinal Process     Event Process
## Number of Observations: 1405 Number of Events: 188 (40.3%)
## Number of Groups: 467
## 
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Relative risk model with spline-approximated
##      baseline risk function
## Parameterization: Time-dependent 
## 
##    log.Lik      AIC      BIC
##  -4334.627 8705.254 8779.888
## 
## Variance Components:
##              StdDev    Corr
## (Intercept)  4.5345  (Intr)
## obstime      0.1685 -0.0589
## Residual     1.8752        
## 
## Coefficients:
## Longitudinal Process
##                   Value Std.Err z-value p-value
## (Intercept)      7.1944  0.1346 53.4509 <0.0001
## obstime         -0.1872  0.0211 -8.8893 <0.0001
## obstime:drugddI  0.0192  0.0272  0.7055  0.4805
## 
## Event Process
##           Value Std.Err z-value p-value
## drugddI  0.3520  0.1578  2.2299  0.0258
## Assoct  -0.2957  0.0380 -7.7885 <0.0001
## bs1     -3.8251  0.5726 -6.6806 <0.0001
## bs2     -1.0674  0.5815 -1.8357  0.0664
## bs3     -3.9097  0.6275 -6.2310 <0.0001
## bs4     -1.2185  0.3767 -3.2351  0.0012
## bs5     -2.7493  0.4508 -6.0981 <0.0001
## bs6     -2.0172  0.7200 -2.8019  0.0051
## bs7     -2.7028  1.4141 -1.9113  0.0560
## bs8     -1.2243  2.7616 -0.4433  0.6575
## bs9     -6.7377  7.3450 -0.9173  0.3590
## 
## Integration:
## method: Gauss-Hermite
## quadrature points: 15 
## 
## Optimization:
## Convergence: 0

Diễn giải kết quả Joint model

Đầu tiên, kết quả cho biết giá trị AIC, BIC và Log likelihood của model

Phần thứ nhất là Fixed effect, khảo sát sự thay đổi của CD4 theo thời gian (phần Longitudinal response), và liệu có sự tương tác với thuốc điều trị (nhóm DDi) hay không ?. Câu trả lời là : CD4 giảm theo thời gian, trung bình -0.187 đơn vị mỗi tháng, tương tác thuốc và thời gian không có ý nghĩa.

Phần thứ hai là khảo sát về biến cố tử vong (phần Survival). Những tham số hồi quy có thể được diễn giải như 1 mô hình Cox cổ điển. Phần này cho ta biết nhiều thông tin thú vị mà bản thân mô hình Cox trước đó không làm được. Đầu tiên, đó là yếu tố Association, hay liên hệ giữa bệnh sử CD4 và thuốc điều trị đối với nguy cơ tử vong :

Ta có thể diễn giải : tại bất cứ thời điểm nào, trên bệnh nhân AIDS sử dụng thuốc DDI , mỗi đơn vị CD4 giảm đi sẽ làm tăng nguy cơ tử vong lên 34% = exp(0.2957)-1

Mô hình sử dụng B-spline bậc 9 (cắt khoảng thời gian 21 tháng ra làm 9 đoạn), tuy nhiên kết quả chỉ có ý nghĩa thống kê tại ps1,ps3,ps4,bs5 và bs6

6 Suy diễn thống kê và kiểm tra mô hình

Ta có thể tính thủ công khoảng tin cậy 97.5% cho các tham số hồi quy (hiệu ứng) quan trọng

# Confidence intervals for the regression coefficients
confint(jointmod, parm = "Longitudinal")

##                       2.5 %        est.      97.5 %
## (Intercept)      6.93054678  7.19435267  7.45815857
## obstime         -0.22851664 -0.18723420 -0.14595175
## obstime:drugddI -0.03412007  0.01919035  0.07250077

# Confidence intervals for the hazard ratios 

confint(jointmod, parm = "Event")%>%exp()

##             2.5 %      est.    97.5 %
## drugddI 1.0435252 1.4218586 1.9373579
## Assoct  0.6906397 0.7439959 0.8014741

Ta có thể làm 1 kiểm định Wald đa biến cho toàn bộ các tham số:

# Multivariate Wald test for simultaneously testing all 

anova(jointmod,process="both")

## 
## Marginal Wald Tests Table
## 
## Longitudinal Process
##                 Chisq df Pr(>|Chi|)
## obstime      117.4799  2     <1e-04
## drug           0.4978  1     0.4805
## obstime:drug   0.4978  1     0.4805
## 
## Event Process
##          Chisq df Pr(>|Chi|)
## drug    4.9724  1     0.0258
## Assoct 60.6612  1     <1e-04

Để kiểm tra phần phân tích sự kiện, ta có thể làm 1 phân tích phương sai : ta dựng1 mô hình Cox tối giản rồi cũng làm joint model với nó, sau đó dùng ANOVA so sánh 2 mô hình joint model cũ và mới

#Anova for Event process

cox2=coxph(Surv(Time, death) ~ 1, data = aids.id, x = TRUE)

jointmod2=jointModel(lmeObject=longcomp, 
                    survObject=cox2,
                    timeVar = "obstime",
                    method = "spline-PH-GH")

anova(jointmod2, jointmod)

## 
##               AIC     BIC  log.Lik  LRT df p.value
## jointmod2 8708.33 8778.82 -4337.16                
## jointmod  8705.25 8779.89 -4334.63 5.08  1  0.0243

Tiếp theo, ta có thể kiểm tra hiệu ứng liên hợp giữa Thuốc và CD4 bằng 1 scoretest

MLmod=update(longcomp, method = "ML")

pwc=piecewiseExp.ph(coxcomp)

init=list(betas = fixef(MLmod), 
             sigma = MLmod$sigma,
             D = getVarCov(MLmod),
             gammas = coef(pwc)[1], alpha = 0,
             xi = exp(coef(pwc)[-1]))

JMScoreTest=jointModel(MLmod, 
                       coxcomp,
                          timeVar = "obstime", 
                          method = "spline-PH-GH", 
                          init = init,
                          only.EM = TRUE)

score.vector=JMScoreTest$Score

inv.Hessian=vcov(JMScoreTest)

ScoreStat=c(t(score.vector) %*% inv.Hessian %*% score.vector)

cbind("Statistic" = ScoreStat, "df" = 1,
      "p-value" = pchisq(ScoreStat, 1, lower.tail = FALSE))

##      Statistic df   p-value
## [1,] 0.8131939  1 0.3671777

Sai số và Marginal effect

plot(jointmod)

#Standardized residuals

resMarg=residuals(jointmod, 
                  process = "Longitudinal",
                           type = "Marginal")

fitMarg=fitted(jointmod, 
                  process = "Longitudinal",
                  type = "Marginal")

myColor=rev(RColorBrewer::brewer.pal(11, "Spectral"))
myColor_scale_fill=scale_fill_gradientn(colours = myColor)
myColor_scale_col=scale_color_gradientn(colours = rev(RColorBrewer::brewer.pal(11, "Spectral")))

aids%>%mutate(.,Fitted=fitMarg,Residual=resMarg)%>%
  ggplot()+
  geom_point(aes(x=Fitted,y=resMarg,color=Fitted))+
  geom_smooth(aes(x=Fitted,y=resMarg),se=F,
              method="lm",formula = y ~ splines::bs(x,3),color="black")+
  geom_hline(yintercept = 0.0,color="blue",linetype=2)+
  theme_bw()+
  myColor_scale_col

aids%>%mutate(.,Fitted=fitMarg,Residual=resMarg)%>%
  ggplot()+
  geom_smooth(aes(x=obstime,y=fitMarg,color=patient),se=F,
              method="lm",formula = y ~ x,show.legend = F)+
  facet_wrap(~drug)+
  theme_bw()

aids%>%mutate(.,Fitted=fitMarg,Residual=resMarg)%>%
  ggplot()+
  geom_smooth(aes(x=obstime,y=fitMarg,color=drug),se=F,
              method="lm",formula = y ~ splines::bs(x))+
  theme_bw()

# Event time marginal effect

EventRes=residuals(jointmod, process = "Event")
mit=fitted(jointmod, process = "Longitudinal",
               type = "EventTime")


aids%>%mutate(.,MartingaleRes=EventRes,ssfvlo=mit)%>%
  ggplot(aes(x=ssfvlo,y=MartingaleRes))+
  geom_point(aes(color=ssfvlo),alpha=0.5,show.legend = F)+
  geom_smooth(se=F,method="lm",formula = y ~ splines::bs(x),color="black")+
  myColor_scale_col+
  geom_hline(yintercept = 0.0,color="blue",linetype=2)+
  theme_bw()+labs(y="Martingale residual",
                  x="Subject specific fitted values longitudinal outcome")

aids%>%mutate(.,MartingaleRes=EventRes,ssfvlo=mit)%>%
  ggplot(aes(x=ssfvlo,y=MartingaleRes))+
  geom_point(aes(color=drug),alpha=0.5,show.legend = F)+
  geom_smooth(se=F,method="lm",formula = y ~ splines::bs(x),color="black")+
  geom_hline(yintercept = 0.0,color="blue",linetype=2)+
  theme_bw()+labs(y="Martingale residual",
                  x="Subject specific fitted values longitudinal outcome")+
  facet_wrap(~drug)

#Residual of event CoxSnell

resCST=residuals(jointmod, 
                 process = "Event",
                    type = "CoxSnell")

sfit=survfit(Surv(resCST,death)~1,data=aids.id)

plot(sfit, mark.time = FALSE, conf.int = TRUE,col=c("red","orange","red4"),
     xlab = "Cox-Snell Residuals", ylab = "Survival Probability",
     main = "Survival Function of Cox-Snell Residuals")

curve(exp(-x), from = 0, to = max(aids.id$Time), add = TRUE,
        col = "blue", lwd = 2)

#Plot2

sfit2=survfit(Surv(resCST, death) ~ drug, data = aids.id)

plot(sfit2, mark.time = FALSE, conf.int = FALSE,col=c("red","blue"),
       xlab = "Cox-Snell Residuals", ylab = "Survival Probability",
       main = "Survival Function of Cox-Snell Residuals")

curve(exp(-x), from = 0, 
      to = max(aids.id$Time), 
      add = TRUE,
      col = "purple", lwd = 2)

7 Ứng dụng: Tiên lượng cho cá thể

Một trong những ứng dụng tuyệt vời của Joint model là để tiên lượng sống còn cho cá thể từng bệnh nhân dựa vào bệnh sử của họ, thí dụ ở tháng thứ 6, dựa vào kết quả CD4 ở 3 lần trước, ta có thể tiên lượng khả năng sống của bệnh nhân đó cho tương lai nhiều tháng sau

Nhi chọn bất kì bệnh nhân nào (dĩ nhiên là phải censored) , thí dụ số 99

#Prediction with CI by MonteCarlo simulation

set.seed(123)

aids[aids$patient=="99",]

##     patient  Time death      CD4 obstime drug gender prevOI         AZT
## 292      99 13.03     0 4.898979       0  ddC   male noAIDS intolerance
## 293      99 13.03     0 1.732051       6  ddC   male noAIDS intolerance
##     start  stop event
## 292     0  6.00     0
## 293     6 13.03     0

#From the last time point

set.seed(333)

survPrbs=survfitJM(jointmod, idVar = "patient",
                   newdata = aids[aids$patient=="99",],
                   last.time = "obstime")

survPrbs

## 
## Prediction of Conditional Probabilities for Event
##  based on 200 Monte Carlo samples
## 
## $`99`
##      times   Mean Median  Lower  Upper
## 1   6.0000 1.0000 1.0000 1.0000 1.0000
## 1   6.0368 0.9987 0.9988 0.9973 0.9995
## 2   6.6553 0.9771 0.9787 0.9519 0.9913
## 3   7.2738 0.9554 0.9584 0.9083 0.9833
## 4   7.8924 0.9323 0.9375 0.8667 0.9746
## 5   8.5109 0.9060 0.9131 0.8076 0.9644
## 6   9.1294 0.8750 0.8834 0.7434 0.9513
## 7   9.7479 0.8384 0.8472 0.6789 0.9402
## 8  10.3665 0.7965 0.8048 0.6113 0.9252
## 9  10.9850 0.7514 0.7615 0.5437 0.9066
## 10 11.6035 0.7070 0.7155 0.4800 0.8859
## 11 12.2221 0.6679 0.6745 0.4244 0.8685
## 12 12.8406 0.6369 0.6435 0.3738 0.8547
## 13 13.4591 0.6125 0.6175 0.3321 0.8434
## 14 14.0776 0.5895 0.5958 0.2965 0.8357
## 15 14.6962 0.5645 0.5691 0.2615 0.8236
## 16 15.3147 0.5361 0.5449 0.2142 0.8151
## 17 15.9332 0.5084 0.5176 0.1689 0.8069
## 18 16.5518 0.4824 0.4874 0.1304 0.7983
## 19 17.1703 0.4570 0.4567 0.1007 0.7872
## 20 17.7888 0.4330 0.4326 0.0810 0.7777
## 21 18.4074 0.4063 0.4015 0.0655 0.7657
## 22 19.0259 0.3764 0.3619 0.0460 0.7559
## 23 19.6444 0.3533 0.3424 0.0284 0.7470
## 24 20.2629 0.3365 0.3352 0.0133 0.7418
## 25 20.8815 0.3172 0.3147 0.0028 0.7398
## 26 21.5000 0.2430 0.2210 0.0000 0.6935

survPrbs[1]$summaries$`99`%>%as_tibble()->tdf
rbind(tdf,c(0,rep(1,4)))->tdf
  tdf%>%
  ggplot(aes(x=times))+
  geom_ribbon(aes(ymin=Lower,ymax=Upper),alpha=0.5,fill="gold")+
  geom_line(aes(y=Median),color="red3",size=1,linetype=1)+
  geom_line(aes(y=Mean),color="blue",size=1,linetype=1)+
  geom_line(aes(y=Lower),color="black",size=1,linetype=2)+
  geom_line(aes(y=Upper),color="black",size=1,linetype=2)+
  theme_bw()+
  labs(title="Patient N°99",x="Time",y="Predicted survival probability",caption="Joint Model")

Hoặc bệnh nhân số 300

survPrbs3=survfitJM(jointmod, idVar = "patient",
                      newdata = aids[aids$patient=="300",],
                      last.time = "obstime")
  
  survPrbs3

## 
## Prediction of Conditional Probabilities for Event
##  based on 200 Monte Carlo samples
## 
## $`300`
##      times   Mean Median  Lower  Upper
## 1   2.0000 1.0000 1.0000 1.0000 1.0000
## 1   2.3256 0.9867 0.9874 0.9732 0.9950
## 2   2.9441 0.9587 0.9611 0.9176 0.9845
## 3   3.5626 0.9291 0.9336 0.8589 0.9730
## 4   4.1812 0.9003 0.9065 0.7991 0.9617
## 5   4.7997 0.8736 0.8814 0.7496 0.9509
## 6   5.4182 0.8493 0.8573 0.7018 0.9407
## 7   6.0368 0.8273 0.8367 0.6576 0.9322
## 8   6.6553 0.8067 0.8189 0.6159 0.9244
## 9   7.2738 0.7866 0.8027 0.5763 0.9171
## 10  7.8924 0.7656 0.7847 0.5309 0.9093
## 11  8.5109 0.7421 0.7640 0.4825 0.9007
## 12  9.1294 0.7147 0.7349 0.4287 0.8865
## 13  9.7479 0.6829 0.7086 0.3703 0.8717
## 14 10.3665 0.6471 0.6664 0.3116 0.8536
## 15 10.9850 0.6094 0.6279 0.2586 0.8391
## 16 11.6035 0.5734 0.5943 0.2101 0.8264
## 17 12.2221 0.5427 0.5619 0.1582 0.8154
## 18 12.8406 0.5193 0.5372 0.1248 0.8063
## 19 13.4591 0.5015 0.5188 0.1065 0.7985
## 20 14.0776 0.4854 0.4976 0.0915 0.7912
## 21 14.6962 0.4683 0.4774 0.0743 0.7837
## 22 15.3147 0.4492 0.4618 0.0565 0.7765
## 23 15.9332 0.4304 0.4417 0.0452 0.7689
## 24 16.5518 0.4124 0.4222 0.0361 0.7613
## 25 17.1703 0.3946 0.4035 0.0284 0.7542
## 26 17.7888 0.3784 0.3848 0.0222 0.7484
## 27 18.4074 0.3613 0.3648 0.0138 0.7441
## 28 19.0259 0.3423 0.3381 0.0079 0.7401
## 29 19.6444 0.3266 0.3227 0.0030 0.7347
## 30 20.2629 0.3132 0.3074 0.0019 0.7225
## 31 20.8815 0.2894 0.2959 0.0001 0.6893
## 32 21.5000 0.2337 0.2218 0.0000 0.6419

  survPrbs3[1]$summaries$`300`%>%as_tibble()->tdf
  rbind(tdf,c(0,rep(1,4)))->tdf
  tdf%>%
    ggplot(aes(x=times))+
    geom_errorbar(aes(ymin=Lower,ymax=Upper,color=reorder(times,Median)),show.legend = F)+
    geom_point(aes(y=Median,fill=reorder(times,Median)),shape=21,color="black",show.legend = F)+
    theme_bw()+
    labs(title="Patient N°300",x="Time",y="Predicted survival probability",caption="Joint Model")

Vì chứa cả 2 thành phần longitudinal và survival, một mô hình Joint model còn cho phép tiên lượng giá trị marker (CD4) theo thời gian:

Thí dụ ta chọn ngẫu nhiên bệnh nhân số 55. Ta có thể tiên lượng CD4 cho tương lai trong 4 hoàn cảnh khác nhau: ngay từ lần hẹn đầu tiên, vào thời điểm tháng thứ 2, 6, và 12

aids[aids$patient==55,]->newdat

longPreds=vector("list", nrow(newdat))

for (i in 1:nrow(newdat)) {
set.seed(123)
    longPreds[[i]]=predict(jointmod, 
                           newdata = newdat[1:i, ],
                           type = "Subject", 
                           interval = "confidence", 
                           returnData = TRUE,
                           idVar="patient")
    longPreds[[i]]$FollowUp=round(max(newdat[1:i, "obstime"]), 1)
  }

longPreds%<>%do.call(rbind,.)%>%as_tibble()

longPreds%>%
  ggplot(aes(x=obstime))+
  geom_ribbon(aes(ymin=low,ymax=upp),alpha=0.5,fill="gold")+
  geom_line(aes(y=pred),color="red3",size=1)+
  geom_line(aes(y=low),color="black",size=1,linetype=2)+
  geom_line(aes(y=upp),color="black",size=1,linetype=2)+
  theme_bw()+
  facet_wrap(~FollowUp)+
  labs(title="Patient N°55",x="Time",y="Predicted CD4",caption="Joint Model")

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

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

Cho bệnh nhân số 55 này, ta có thể kiểm tra độ chính xác của tiên lượng trong 4 hoàn cảnh: vào tháng thứ 2,6,12 và 18

# Performance evaluation

set.seed(123)

newdat$t0=as.numeric(newdat$obstime == 0)

rocaids= rocJM(jointmod, 
               dt = c(2,6,12,18), 
               data = newdat, 
               optThr="youden",
               M=1000,burn.in=500,
               idVar="patient")

plot(rocaids, legend = TRUE)

rocaids

## 
## Areas under the time-dependent ROC curves
## 
## Estimation: Monte Carlo (500 samples)
## Difference: absolute, lag = 1 (0)
## Thresholds range: (-1.41, 40.89)
## 
## Case: 55 
## Recorded time(s): 0, 2, 6, 12 
##  dt t + dt    AUC    Cut
##   2     14 0.7988 5.3541
##   6     18 0.8281 5.5233
##  12     24 0.8594 7.2154
##  18     30 0.8817 8.7383

8 Kết luận:

Trong bài này, các bạn đã làm quen với Joint model, một trong những công cụ cần thiết cho nghiên cứu lâm sàng.

Mô hình liên hợp (Joint model) cho phép chúng ta hợp nhất hai bộ phận tưởng chừng như độc lập trong một thử nghiệm lâm sàng, đó là bài toán Longitudinal analysis và bài toán Survival analysis. Việc hợp nhất này cho phép giải quyết hàng loạt câu hỏi, như :

• Chứng minh về liên hệ giữa bệnh sử thông qua 1 marker và nguy cơ tử vong,

• chứng minh 1 surrogate marker có thể dùng để đánh giá hiệu quả điều trị của thuốc,

• tiên lượng diễn tiến của marker và nguy cơ tử vong cho từng bệnh nhân riêng biệt một khi đã biết bệnh sử của cá nhân này.

HY vọng bài viết tạo được cảm hứng cho các bạn bắt đầu tìm hiểu về phương pháp này.

Một bài khác sẽ được thực hiện để đi theo con đường Bayes.

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