Mô hình liên hợp

Phần 2: Bayesian Joint model

1 Dẫn nhập

Trong bài trước, chúng ta đã làm quen với package JM cho phép dựng mô hình Joint model bằng phương pháp REML. Như đã hứa, Nhi sẽ giới thiệu tiếp một package khác (của cùng tác giả Rizoupoulos) cho phép dựng Joint model bằng phương pháp Bayes.

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.

Để tóm tắt lại về nội dung mô hình Joint model, những thành phần sau đây được xét:

1. Ti = Thời gian theo dõi đến khi sự kiện xảy ra cho một cá thể (i) xác định

2. Ci = Thời gian cho đến khi cá thể (i) được loại khỏi vòng theo dõi

3. Trạng thái sự kiện xảy ra là 1 biến nhị phân (0,1) được xác định bằng 1 hàm với quy luật: kết quả = 1 nếu Ti < hoặc = Ci, và =0 nếu Ti > Ci

4. Bệnh sử cá thể (yi) là một vector chứa những giá trị của marker M theo thời gian cho một cá thể xác định. Mỗi phần tử yij trong vector này tương ứng với 1 giá trị của M cho cá thể i tại thời điểm tj (j là số lần khảo sát)

Thành phần longitudinal

Mục tiêu đầu tiên đó là ước tính giá trị Mu của yi cho một cá thể xác định, tại 1 thời điểm t xác định. Điều này có thể được thực hiện thông qua một mô hình GLM mở rộng. Nội dung của mô hình này gồm có :

1. biến kết quả Y(it) được mô tả bằng 1 phân phối xác suất có điều kiện phụ thuộc vào 2 yếu tố: thời gian (t) và hiệu ứng ngẫu nhiên (bi). Vì đây là 1 mô hình GLM, ta có thể dùng một hàm link function G để liên kết yi(t) với design matrix và tham số hồi quy của mô hình.

2. Vế bên phải của mô hình gồm có 2 phần : hiệu ứng chính (fixed effect) gồm design matrix (các biến độc lập, phụ thuộc thời gian) xi(t) + vector tham số hồi quy beta, và hiệu ứng ngẫu nhiên gồm design matrix Zi(t) với các yếu tố ngẫu nhiên và vector random effect bi.

3. Hiệu ứng ngẫu nhiên (vector bi) được giả định có phân phối chuẩn với trung bình =zero và matrix cariance covariance D.

Thành phần survival

Mục tiêu tiếp theo, đó là ước tính nguy cơ xảy ra sự kiện (tử vong) tại thời điểm t xác định và cho một cá thể i xác định, biết rằng nguy cơ này phụ thuộc vào đồng thời : một tập hợp hiệp biến số (covariate) Wi và một vector bệnh sử H.

H ở đây chính là kết quả ước tính giá trị marker cho cá thể đó trong khoảng thời gian s từ 0 cho đến thời điểm t.

Hàm nguy cơ (hi) có thể được phân tích thành 1 hàm nguy cơ gốc h0 cho thời gian và một hàm exponential chứa lần lượt : vector tham số hồi quy Gamma tương ứng với tập hợp hiệp biến số Wi và một hàm f biểu thị cho hiệu ứng « liên hệ » giữa giá trị marker (Mu, ước tính qua mô hình Longitidinal hay mixed model), hiệu ứng ngẫu nhiên bi và vector tham số hồi quy Alpha, cho hiệu ứng quan hệ (association) giữa marker và covariates khác trong mô hình.

Hàm nguy cơ gốc h0 thường được ước lượng với B-spline (là một thủ thuật để mô hình hóa quan hệ phi tuyến tính trong mô hình)

Như vậy có thể thấy kết quả chính của mô hình liên hợp gồm 1 số vector tham số hồi ququan trọng y như alpha (hiệu ứng liên hợp, association, cho bộ phận survival), beta (hiệu ứng chính, fixed effect cho bộ phận longitudinal) , gamma (hiệu ứng độc lập, cho bộ phận survival).

Về cơ bản, package JMBayes có cấu trúc method giống như package JM (sử dụng REML). Phương pháp Bayes chỉ áp dụng trên joint model, nhưng 2 mô hình thành phần Longtudinal và Survival vẫn được dựng theo cách truyền thống từ bên trong hoặc bên ngoài (JM và JMbayes tương thích với mixed model object của package nlme (2014) hay MASS và Coxph model của packae survival (2014). Các hàm lme, coxph() và cú pháp của chúng được tích hợp sẵn trong JM và JM Bayes. Hàm jointModelBayes() sẽ dựng joint model.

Khi sử dụng cách tiếp cận theo trường phái Bayes, ta muốn xác định phân phối hậu nghiệm của toàn bộ các vector tham số trong joint model, bao gồm cả vector random effect bi.

Một cách đơn giản, Nhi chỉ có thể giải thích là tác giả Rizopoulos đã dùng những chuỗi Markov Monte Carlo (MCMC) với algorithm Metropolis Hastings. Những chuỗi MCMC này có thể được trích xuất sau đó để khảo sát.

Ngoài ra JPBayes có một số chứng năng nâng cao so với JM, nó cho phép mở rộng hiệu ứng liên hệ trong joint model với những tham số pho tuyến tính bổ sung, nó cũng cho phép làm Bayesian model Averaging (BMA) để kết hợp nhiều mô hình Joint model khác nhau cho việc tiên lượng.

Chi tiết kỹ thuật của package có thể tham khảo trong bài báo:

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.409.4524&rep=rep1&type=pdf

Tác giả cũng từng làm 1 webminar năm 2014:

https://www.youtube.com/watch?v=vqGpQmnup0Y

Bài thực hành này chỉ trình bày một thí dụ cơ bản với 1 covariate và Association effect đơn giản.

2 Minh họa: Prothrombin và Prednisone/ Bệnh Xơ gan

Thí dụ minh họa trong bài sử dụng dataset prothro. Đây là một nghiên cứu thử nghiệm lâm sàng trên 488 bệnh nhân Xơ gan. Có hai biến số outcome, một là maker prothrombin (một yếu tố đông máu), được khảo sát tại thời điểm zero, sau đó lặp lại mỗi tháng, hai là sự kiện tử vong. Mỗi bệnh nhân có thể được đo prothrombin tối đa đến 17 lần. Yếu tố can thiệp là 2 phương pháp điều trị bằng prednisone và placebo. Như bài trước, mục tiêu của chúng ta là khảo sát liên hệ giữa sự thay đổi của marker prothrombin theo thời gian và yếu tố điều trị đối với nguy cơ tử vong của bệnh nhân xơ gan.

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

library(JMbayes) # package dựng Joint model bằng phương pháp Bayes
library(tidyverse) # data wrangling và ggplot
library(magrittr) # Toán tử pipe
library(ggkm) #Thăm dò Time event data

3 Thăm dò dữ liệu

Trước hết, chúng ta viết một vài đoạn code để mô tả dữ liệu bằng những con số. Một số câu hỏi có thể đặt ra :

1. Cấu trúc dữ liệu ? Có 2 định dạng khác nhau của dataset, 1 cho longitudinal study và 1 cho survival study.

Ghi chú: id là mã số định danh cho bệnh nhân; pro là giá trị prothrombin, treat là phân nhóm điều trị, time là thời điểm mỗi visit tính từ bắt đầu nghiên cứu; Time là thời gian theo dõi tối đa đến khi bệnh nhân tử vong hoặc censored, start là mốc thời gian khởi đầu mỗi phân đoạn theo dõi, stop là mốc thời gian kết thúc mỗi phân đoạn, death là biến số xác nhận bệnh nhân tử vong hoặc censored, event là trạng thái sự kiện tại mỗi thời điểm.

prothro%>%filter(id==c("2","8"))%>%knitr::kable()

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

id	pro	time	treat	Time	start	stop	death	event
2	73	0.6872194	prednisone	6.754463	0.6872194	0.9610119	1	0
2	64	1.1882598	prednisone	6.754463	1.1882598	1.4428869	1	0
2	58	1.7139415	prednisone	6.754463	1.7139415	1.9959479	1	0
2	75	2.7214982	prednisone	6.754463	2.7214982	3.7728617	1	0
2	45	4.7503012	prednisone	6.754463	4.7503012	5.7167890	1	0
2	99	6.6997043	prednisone	6.754463	6.6997043	6.7544628	1	1
8	108	0.5777023	placebo	7.584054	0.5777023	1.0842186	0	0
8	100	1.3032527	placebo	7.584054	1.3032527	1.6071624	0	0
8	120	1.8042931	placebo	7.584054	1.8042931	2.0534443	0	0
8	100	3.0500493	placebo	7.584054	3.0500493	4.0165371	0	0
8	100	5.0514730	placebo	7.584054	5.0514730	6.0152229	0	0
8	130	6.9926624	placebo	7.584054	6.9926624	7.5840543	0	0

prothros%>%head()%>%knitr::kable()

id	Time	death	treat
1	0.4134268	1	prednisone
2	6.7544628	1	prednisone
3	13.3939328	0	prednisone
4	0.7939985	1	prednisone
5	0.7501917	1	prednisone
6	0.7693571	1	prednisone

2. Tổng số bệnh nhân ? Con số này là 488

prothro$id%>%
  unique() %>%
  length()

## [1] 488

3. Có bao nhiêu bệnh nhân đã tử vong ? Thời gian theo dõi trung bình của mỗi bệnh nhân là bao nhiêu:

Kết quả cho thấy thời gian trung bình đến khi mất theo dõi hay kết thúctheo dõi là 4.8 tháng (cao nhất 13.39, thấp nhất 0.06), thời gian trung bình đến khi tử vong là 2.87 tháng (0.01 đến 10.33 tháng) Tỉ lệ tử vong/censored là 292/196

prothros%>%mutate(.,Death=ifelse(.$death==1,"Dead","Censored"))%>%
{psych::describeBy(.$Time,group=.$Death)}

## 
##  Descriptive statistics by group 
## group: Censored
##    vars   n mean   sd median trimmed  mad  min   max range skew kurtosis
## X1    1 196  4.8 3.68    5.2    4.63 5.53 0.06 13.39 13.33 0.21    -1.29
##      se
## X1 0.26
## -------------------------------------------------------- 
## group: Dead
##    vars   n mean  sd median trimmed  mad  min   max range skew kurtosis
## X1    1 292 2.87 2.7   2.05    2.51 2.43 0.01 10.33 10.32 0.91    -0.24
##      se
## X1 0.16

4. Mỗi phân nhóm điều trị có bao nhiêu bệnh nhân ? Họ được theo dõi trong bao lâu ?

Nhóm điều trị với Prednisone có 237 bệnh nhân, nhóm Placebo có 251 bệnh nhân. Thời gian theo dõi trung bình tương ứng là 3.63 và 3.66 tháng

prothros%>%mutate(.,Death=ifelse(.$death==1,"Dead","Censored"))%>%
  {psych::describeBy(.$Time,.$treat)}

## 
##  Descriptive statistics by group 
## group: placebo
##    vars   n mean   sd median trimmed  mad  min   max range skew kurtosis
## X1    1 251 3.63 3.35    2.6     3.3 3.41 0.01 12.18 12.17 0.62    -0.97
##      se
## X1 0.21
## -------------------------------------------------------- 
## group: prednisone
##    vars   n mean   sd median trimmed  mad  min   max range skew kurtosis
## X1    1 237 3.66 3.19   2.64    3.33 3.17 0.04 13.39 13.36 0.76    -0.45
##      se
## X1 0.21

5. Phân bố các trường hợp tử vong ?

prothros%>%mutate(.,Death=ifelse(.$death==1,"Dead","Censored"))%>%
  xtabs(data=.,~Death+treat)

##           treat
## Death      placebo prednisone
##   Censored     109         87
##   Dead         142        150

6. Giá trị prothrombin cơ bản lúc bắt đầu nghiên cứu là bao nhiêu ?

paste("Prothrombin T0: Giữa 2 phân nhóm điều trị")

## [1] "Prothrombin T0: Gi<U+1EEF>a 2 phân nhóm di<U+1EC1>u tr<U+1ECB>"

prothro%>%filter(.$time==0)%>%
{psych::describeBy(.$pro,.$treat)}

## 
##  Descriptive statistics by group 
## group: placebo
##    vars   n  mean    sd median trimmed   mad min max range skew kurtosis
## X1    1 251 68.65 23.88     67   68.03 26.69  12 134   122 0.31    -0.38
##      se
## X1 1.51
## -------------------------------------------------------- 
## group: prednisone
##    vars   n  mean    sd median trimmed  mad min max range skew kurtosis
## X1    1 237 68.83 22.14     68   68.75 25.2  16 126   110 0.07    -0.68
##      se
## X1 1.44

paste("Prothrombin T0: Giữa Có/Không có tử vong")

## [1] "Prothrombin T0: Gi<U+1EEF>a Có/Không có t<U+1EED> vong"

prothro%>%filter(.$time==0)%>%
  mutate(.,Death=ifelse(.$death==1,"Dead","Censored"))%>%
  {psych::describeBy(.$pro,.$Death)}

## 
##  Descriptive statistics by group 
## group: Censored
##    vars   n  mean    sd median trimmed   mad min max range  skew kurtosis
## X1    1 196 75.21 22.29   73.5    75.7 24.46  12 134   122 -0.09     -0.1
##      se
## X1 1.59
## -------------------------------------------------------- 
## group: Dead
##    vars   n mean    sd median trimmed   mad min max range skew kurtosis
## X1    1 292 64.4 22.52     61   63.35 23.72  16 134   118 0.43    -0.41
##      se
## X1 1.32

Bây giờ ta sẽ trình bày thông tin bằng biểu đồ:

Hình 1) Số lần visit ở mỗi bệnh nhân có thể dao động từ 1 đến 17 lần

prothro%>%mutate(.,Death=ifelse(.$death==1,"Dead","Censored"))%>%
  ggplot(aes(x=id,fill=Death))+
  stat_bin(binwidth=1,show.legend = F)+
  stat_bin(binwidth=1, geom="text", aes(label=..count..),size=3)+
  theme_bw()+
  coord_flip()+
  facet_grid(treat~Death)+
  labs(title="Visits count",x="Patient Id",y="Number of Visits")

Hình 2: Diễn tiến của prothrombin rất khác nhau ở mỗi bệnh nhân, nhưng nhìn chung : những bệnh nhân tử vong có prothrombin thấp hơn, phân nhóm prednisone có prothrombin thấp hơn phân nhóm placebo

library(colorRamps)
library(RColorBrewer) 

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

prothro%>%mutate(.,Death=ifelse(.$death==1,"Dead","Censored"))%>%
  ggplot()+
  geom_path(aes(x=time,y=pro,color=as.factor(id)),show.legend = F,alpha=0.6)+
  geom_smooth(aes(x=as.numeric(time),y=pro,fill=treat,linetype=treat),
              color="black",
              method="lm",formula = y ~ splines::bs(x,3),alpha=0.4)+
  scale_color_manual(values = getPalette(mypalette))+
  facet_wrap(~Death)+
  theme_bw()

Hình 3a,b): Giả sử ta chia đều thời gian theo dõi thành 11 vị trí: Giá trị prothrombin trung bình tại mỗi thời điểm là như sau:

prothro%>%mutate(.,Death=ifelse(.$death==1,"Death","Survived"))%>%
  ggplot(aes(y=pro,x=death,fill=treat))+
  geom_boxplot(alpha=0.5)+
  coord_flip()+
  facet_wrap(~round(time,0),ncol=3)+
  theme_bw()

Như vậy tại mỗi thời điểm, có vẻ như không có sự tương phản rõ rệt giữa 2 phân nhóm điều trị

prothro%>%mutate(.,Death=ifelse(.$death==1,"Death","Survived"))%>%
  ggplot(aes(x=pro,fill=treat))+
  geom_density(alpha=0.5)+
  facet_wrap(~round(time,0),scales="free_y",ncol=3)+
  theme_bw()

Hình 4a,b) Mô hình random intercept và random slope với natural cubic spline:

Nếu xét riêng bài toán longitudinal cho prothrombin, ta có thể dùng mô hình mixed model (random intercept hoặc random slope). Trong mô hình này biến kết quả là prothrombinở thang đo logarit, phụ thuộc vào thời gian (hay nói cách khác: thời điểm lấy mẫu = biến time). Mặt khác, theo quan sát ban đầu ta có thể nhận ra là diễn tiến của prothrombin (biến kết quả) theo thời gian đi theo một lộ trình phi tuyến tính, do đó ta có thể “uốn cong” đồ thị của mô hình với một hàm natural cubic spline (hàm này được tích hợp sẵn trong package JMbayes). Kết quả của mô hình như sau

Mô hình random intercept với natural cubic spline

rint=lme(log(pro) ~ ns(time,3), 
         random = ~ 1 | id,
         data = prothro)

prothro%>%mutate(.,Fitted=predict(rint))%>%
  ggplot()+
  geom_path(aes(x=as.numeric(time),y=Fitted,color=as.factor(id)),
            show.legend = F,
            alpha=0.1)+
  geom_smooth(aes(x=as.numeric(time),y=Fitted),
              method="lm",alpha=0.5,
              formula=y~ns(x,3),se=F,
            show.legend = F,color="red")+
  scale_color_grey()+
  theme_bw()+
  labs(title="Random Intercept model: Y~ns(time,3)+1|Id",x="Time",y="Predicted prothromb (Log scale)")

Mô hình random slope với natural cubic spline

rslope=lme(log(pro) ~ ns(time,3), 
         random = ~ ns(time,3) | id,
         data = prothro)

prothro%>%mutate(.,Fitted=predict(rslope))%>%
  ggplot()+
  geom_path(aes(x=as.numeric(time),
                y=Fitted,color=as.factor(id)),
            show.legend = F,
            alpha=0.1)+
  geom_smooth(aes(x=as.numeric(time),y=Fitted),
              method="lm",alpha=0.5,
              formula=y~ns(x,3),se=F,
            show.legend = F,color="red")+
  scale_color_grey()+
  theme_bw()+
  labs(title="Random Slope model: Y~ns(time,3)+ ns(time,3)|Id",x="Time",y="Predicted prothromb (Log scale)")

HÌnh 5 a,b,c,d):

Tiếp theo chúng ta xét đến bài toán về nguy cơ tử vong (survival study); 3 hình thức khác nhau của biểu đồ Kaplan-Meier được trình bày trong hình 5, chúng mô tả hàm survival, hàm nguy cơ và hàm CDF theo thời gian phân biệt cho 2 nhóm điều trị.

Nếu không xét đến yếu tố prothrombin, nguy cơ tử vong có vẻ cao hơn ở nhóm prednisone so với nhóm placebo

library(survminer)

## Loading required package: ggpubr

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

library(survival)

fit <- survfit(Surv(Time, death) ~ treat, data = prothro)

ggsurvplot(fit, 
           data = prothro,
           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
)

# KM plots

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

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

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

Hình 6a) Ta cũng có thể so sánh nguy cơ tử vong giữa 2 phân nhóm có giá trị Prothrombin ở thời điểm T0 cao hơn và thấp hơn ngưỡng trung vị trong mẫu nghiên cứu:

Các bệnh nhân có prothrombin ban đầu (T0) thấp hơn trung vị thì có nhiều nguy cơ tử vong hơn.

prothro%>%filter(time==0.0)%>%mutate(.,Pro=ifelse(.$pro>=median(.$pro),"High","Low"))%>%
ggplot(aes(time = Time, status = death, fill = Pro, color = Pro))+ 
  geom_km(trans = "cumhaz")+ 
  geom_kmband(trans = "cumhaz")+
  theme_bw()+
  facet_wrap(~treat)+
  ylab("Cumulative Hazard")

Tương tự, nếu chỉ xét giá trị prothrombin ở thời điểm cuối cùng (ngay trước khi censore hay tử vong): giá trị prothrombin thấp cũng có liên hệ với nguy cơ tử vong:

prothro%>%filter(Time==stop)%>%mutate(.,Pro=ifelse(.$pro>=median(.$pro),"High","Low"))%>%
  ggplot(aes(time = Time, status = death, fill = Pro, color = Pro))+ 
  geom_km(trans = "cumhaz")+ 
  geom_kmband(trans = "cumhaz")+
  theme_bw()+
  facet_wrap(~treat)+
  ylab("Cumulative Hazard")

Tuy nhiên, ta không thể chỉ khảo sát marker Prothrombin duy nhất tại 2 thời điểm Đầu và Cuối, mà ta muốn khảo sát toàn bộ diễn tiến của marker này và kiểm tra liệu có liên hệ giữa diễn tiến này với yếu tố điều trị gây thay đổi nguy cơ tử vong hay không ?Do đó ta sẽ sử dụng Joint model

4 Dựng Joint Model theo phương pháp Bayes

Mô hình bộ phận thứ nhất: Longitudinal

Đầu tiên, chúng ta dựng một mô hình tuyến tính hỗn hợp (Mixed model) cho diễn tiến của prothrombin theo thời gian. Biến kết quả của mô hình là giá trị prothrombin ở thang đo logarit. Như đã nói ở trên, sự thay đổi theo thời gian ở đa số bệnh nhân là không tuyến tính, do đó Nhi sẽ sử dụng hàm natural cubic spline cho cả 2 hiệu ứng: Fixed và random trong mô hình. Như vậy đây chính là 1 mô hình random slope.

lmeFit=lme(log(pro)~ns(time, 3),
                    data = prothro,
                    random = ~ ns(time, 3) | id)

summary(lmeFit)

## Linear mixed-effects model fit by REML
##  Data: prothro 
##        AIC      BIC    logLik
##   1735.742 1825.657 -852.8711
## 
## Random effects:
##  Formula: ~ns(time, 3) | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##              StdDev    Corr                
## (Intercept)  0.3013339 (Intr) n(,3)1 n(,3)2
## ns(time, 3)1 0.6194410 -0.212              
## ns(time, 3)2 0.4032879 -0.142  0.439       
## ns(time, 3)3 0.1879770  0.182 -0.387  0.447
## Residual     0.2495787                     
## 
## Fixed effects: log(pro) ~ ns(time, 3) 
##                  Value  Std.Error   DF   t-value p-value
## (Intercept)   4.184412 0.01676300 2477 249.62190  0.0000
## ns(time, 3)1 -0.104434 0.04902390 2477  -2.13027  0.0332
## ns(time, 3)2  0.281738 0.03982575 2477   7.07428  0.0000
## ns(time, 3)3  0.308066 0.04887171 2477   6.30356  0.0000
##  Correlation: 
##              (Intr) n(,3)1 n(,3)2
## ns(time, 3)1 -0.068              
## ns(time, 3)2 -0.379  0.031       
## ns(time, 3)3 -0.077 -0.397  0.659
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -8.33432522 -0.41654029  0.07603526  0.51689202  3.78977554 
## 
## Number of Observations: 2968
## Number of Groups: 488

plot(lmeFit)

mypalette = length(unique(prothro$id))
getPalette = colorRampPalette(brewer.pal(9,"RdBu"))

prothro%>%mutate(.,Fitted=predict(lmeFit))%>%
  ggplot()+
  geom_point(aes(x=as.numeric(time),
                y=log(pro),color=as.factor(id)),alpha=0.2,show.legend = F)+
  geom_path(aes(x=as.numeric(time),
                y=Fitted,color=as.factor(id)),
            show.legend = F,
            alpha=0.2)+
  geom_smooth(aes(x=as.numeric(time),y=Fitted),
              method="lm",alpha=0.5,
              formula=y~ns(x,3),se=F,
            show.legend = F,color="black")+
  scale_color_manual(values = rev(getPalette(mypalette)))+
  theme_bw()+
  labs(title="Random Slope model: Y~ns(time,3)+ ns(time,3)|Id",x="Time",y="Predicted prothromb (Log scale)")

Thành phần thứ hai trong Joint model là một mô hình CoxPH, chỉ chứa yếu tố phân nhóm điều trị.

coxFit=coxph(Surv(Time,death) ~ treat, 
                  data = prothros,
                  x = TRUE)

summary(coxFit)

## Call:
## coxph(formula = Surv(Time, death) ~ treat, data = prothros, x = TRUE)
## 
##   n= 488, number of events= 292 
## 
##                    coef exp(coef) se(coef)     z Pr(>|z|)
## treatprednisone 0.09989   1.10504  0.11720 0.852    0.394
## 
##                 exp(coef) exp(-coef) lower .95 upper .95
## treatprednisone     1.105     0.9049    0.8782      1.39
## 
## Concordance= 0.51  (se = 0.016 )
## Rsquare= 0.001   (max possible= 0.999 )
## Likelihood ratio test= 0.73  on 1 df,   p=0.394
## Wald test            = 0.73  on 1 df,   p=0.3941
## Score (logrank) test = 0.73  on 1 df,   p=0.3939

Sau đó, ta sẽ dựng Joint Model cơ bản với package JMbayes với 50 ngàn lượt lấy mẫu cho chuỗi MCMC

jointfit=jointModelBayes(lmeFit,coxFit,timeVar = "time",n.iter=50000)

## 
##  MCMC iterations:
## 
## 
  |                                                        
  |                                                  |   0%
  |                                                        
  |                                                  |   1%
  |                                                        
  |+                                                 |   1%
  |                                                        
  |+                                                 |   2%
  |                                                        
  |+                                                 |   3%
  |                                                        
  |++                                                |   3%
  |                                                        
  |++                                                |   4%
  |                                                        
  |++                                                |   5%
  |                                                        
  |+++                                               |   5%
  |                                                        
  |+++                                               |   6%
  |                                                        
  |+++                                               |   7%
  |                                                        
  |++++                                              |   7%
  |                                                        
  |++++                                              |   8%
  |                                                        
  |++++                                              |   9%
  |                                                        
  |+++++                                             |   9%
  |                                                        
  |+++++                                             |  10%
  |                                                        
  |+++++                                             |  11%
  |                                                        
  |++++++                                            |  11%
  |                                                        
  |++++++                                            |  12%
  |                                                        
  |++++++                                            |  13%
  |                                                        
  |+++++++                                           |  13%
  |                                                        
  |+++++++                                           |  14%
  |                                                        
  |+++++++                                           |  15%
  |                                                        
  |++++++++                                          |  15%
  |                                                        
  |++++++++                                          |  16%
  |                                                        
  |++++++++                                          |  17%
  |                                                        
  |+++++++++                                         |  17%
  |                                                        
  |+++++++++                                         |  18%
  |                                                        
  |+++++++++                                         |  19%
  |                                                        
  |++++++++++                                        |  19%
  |                                                        
  |++++++++++                                        |  20%
  |                                                        
  |++++++++++                                        |  21%
  |                                                        
  |+++++++++++                                       |  21%
  |                                                        
  |+++++++++++                                       |  22%
  |                                                        
  |+++++++++++                                       |  23%
  |                                                        
  |++++++++++++                                      |  23%
  |                                                        
  |++++++++++++                                      |  24%
  |                                                        
  |++++++++++++                                      |  25%
  |                                                        
  |+++++++++++++                                     |  25%
  |                                                        
  |+++++++++++++                                     |  26%
  |                                                        
  |+++++++++++++                                     |  27%
  |                                                        
  |++++++++++++++                                    |  27%
  |                                                        
  |++++++++++++++                                    |  28%
  |                                                        
  |++++++++++++++                                    |  29%
  |                                                        
  |+++++++++++++++                                   |  29%
  |                                                        
  |+++++++++++++++                                   |  30%
  |                                                        
  |+++++++++++++++                                   |  31%
  |                                                        
  |++++++++++++++++                                  |  31%
  |                                                        
  |++++++++++++++++                                  |  32%
  |                                                        
  |++++++++++++++++                                  |  33%
  |                                                        
  |+++++++++++++++++                                 |  33%
  |                                                        
  |+++++++++++++++++                                 |  34%
  |                                                        
  |+++++++++++++++++                                 |  35%
  |                                                        
  |++++++++++++++++++                                |  35%
  |                                                        
  |++++++++++++++++++                                |  36%
  |                                                        
  |++++++++++++++++++                                |  37%
  |                                                        
  |+++++++++++++++++++                               |  37%
  |                                                        
  |+++++++++++++++++++                               |  38%
  |                                                        
  |+++++++++++++++++++                               |  39%
  |                                                        
  |++++++++++++++++++++                              |  39%
  |                                                        
  |++++++++++++++++++++                              |  40%
  |                                                        
  |++++++++++++++++++++                              |  41%
  |                                                        
  |+++++++++++++++++++++                             |  41%
  |                                                        
  |+++++++++++++++++++++                             |  42%
  |                                                        
  |+++++++++++++++++++++                             |  43%
  |                                                        
  |++++++++++++++++++++++                            |  43%
  |                                                        
  |++++++++++++++++++++++                            |  44%
  |                                                        
  |++++++++++++++++++++++                            |  45%
  |                                                        
  |+++++++++++++++++++++++                           |  45%
  |                                                        
  |+++++++++++++++++++++++                           |  46%
  |                                                        
  |+++++++++++++++++++++++                           |  47%
  |                                                        
  |++++++++++++++++++++++++                          |  47%
  |                                                        
  |++++++++++++++++++++++++                          |  48%
  |                                                        
  |++++++++++++++++++++++++                          |  49%
  |                                                        
  |+++++++++++++++++++++++++                         |  49%
  |                                                        
  |+++++++++++++++++++++++++                         |  50%
  |                                                        
  |+++++++++++++++++++++++++                         |  51%
  |                                                        
  |++++++++++++++++++++++++++                        |  51%
  |                                                        
  |++++++++++++++++++++++++++                        |  52%
  |                                                        
  |++++++++++++++++++++++++++                        |  53%
  |                                                        
  |+++++++++++++++++++++++++++                       |  53%
  |                                                        
  |+++++++++++++++++++++++++++                       |  54%
  |                                                        
  |+++++++++++++++++++++++++++                       |  55%
  |                                                        
  |++++++++++++++++++++++++++++                      |  55%
  |                                                        
  |++++++++++++++++++++++++++++                      |  56%
  |                                                        
  |++++++++++++++++++++++++++++                      |  57%
  |                                                        
  |+++++++++++++++++++++++++++++                     |  57%
  |                                                        
  |+++++++++++++++++++++++++++++                     |  58%
  |                                                        
  |+++++++++++++++++++++++++++++                     |  59%
  |                                                        
  |++++++++++++++++++++++++++++++                    |  59%
  |                                                        
  |++++++++++++++++++++++++++++++                    |  60%
  |                                                        
  |++++++++++++++++++++++++++++++                    |  61%
  |                                                        
  |+++++++++++++++++++++++++++++++                   |  61%
  |                                                        
  |+++++++++++++++++++++++++++++++                   |  62%
  |                                                        
  |+++++++++++++++++++++++++++++++                   |  63%
  |                                                        
  |++++++++++++++++++++++++++++++++                  |  63%
  |                                                        
  |++++++++++++++++++++++++++++++++                  |  64%
  |                                                        
  |++++++++++++++++++++++++++++++++                  |  65%
  |                                                        
  |+++++++++++++++++++++++++++++++++                 |  65%
  |                                                        
  |+++++++++++++++++++++++++++++++++                 |  66%
  |                                                        
  |+++++++++++++++++++++++++++++++++                 |  67%
  |                                                        
  |++++++++++++++++++++++++++++++++++                |  67%
  |                                                        
  |++++++++++++++++++++++++++++++++++                |  68%
  |                                                        
  |++++++++++++++++++++++++++++++++++                |  69%
  |                                                        
  |+++++++++++++++++++++++++++++++++++               |  69%
  |                                                        
  |+++++++++++++++++++++++++++++++++++               |  70%
  |                                                        
  |+++++++++++++++++++++++++++++++++++               |  71%
  |                                                        
  |++++++++++++++++++++++++++++++++++++              |  71%
  |                                                        
  |++++++++++++++++++++++++++++++++++++              |  72%
  |                                                        
  |++++++++++++++++++++++++++++++++++++              |  73%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++             |  73%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++             |  74%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++             |  75%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++            |  75%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++            |  76%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++            |  77%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++           |  77%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++           |  78%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++           |  79%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++          |  79%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++          |  80%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++          |  81%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++         |  81%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++         |  82%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++         |  83%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++        |  83%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++        |  84%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++        |  85%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++       |  85%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++       |  86%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++       |  87%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++      |  87%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++      |  88%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++      |  89%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++     |  89%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++     |  90%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++     |  91%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++    |  91%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++    |  92%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++    |  93%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++++   |  93%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++++   |  94%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++++   |  95%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++  |  95%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++  |  96%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++  |  97%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++++++ |  97%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++++++ |  98%
  |                                                        
  |+++++++++++++++++++++++++++++++++++++++++++++++++ |  99%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++++|  99%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%

summary(jointfit)

## 
## Call:
## jointModelBayes(lmeObject = lmeFit, survObject = coxFit, timeVar = "time", 
##     n.iter = 50000)
## 
## Data Descriptives:
## Longitudinal Process     Event Process
## Number of Observations: 2968 Number of Events: 292 (59.8%)
## Number of subjects: 488
## 
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Relative risk model with penalized-spline-approximated 
##      baseline risk function
## Parameterization: Time-dependent value 
## 
##       LPML      DIC       pD
##  -3788.693 7155.456 1931.453
## 
## Variance Components:
##               StdDev    Corr                
## (Intercept)   0.4717  (Intr)  n(,3)1  n(,3)2
## ns(time, 3)1  0.7863 -0.0297                
## ns(time, 3)2  0.7370 -0.0572 -0.0743        
## ns(time, 3)3  1.3559 -0.0658  0.0789 -0.1094
## Residual      0.2372                        
## 
## Coefficients:
## Longitudinal Process
##                Value Std.Err Std.Dev    2.5%   97.5%      P
## (Intercept)   4.1845  0.0006  0.0228  4.1412  4.2290 <0.001
## ns(time, 3)1 -0.1209  0.0015  0.0454 -0.2146 -0.0359   0.01
## ns(time, 3)2  0.2265  0.0013  0.0415  0.1423  0.3044 <0.001
## ns(time, 3)3  0.3770  0.0032  0.0765  0.2314  0.5336 <0.001
## 
## Event Process
##                    Value Std.Err  Std.Dev    2.5%    97.5%      P
## treatprednisone  -0.1617  0.0045   0.1219 -0.3977   0.0769  0.197
## Assoct           -2.4639  0.0293   0.1646 -2.7848  -2.1415 <0.001
## tauBs           255.3340 11.6116 151.4681 57.4600 638.9329     NA
## 
## MCMC summary:
## iterations: 50000 
## adapt: 3000 
## burn-in: 3000 
## thinning: 25 
## time: 12.1 min

jointfit$mcmc%>%names()

## [1] "betas"     "sigma"     "b"         "D"         "gammas"    "Bs.gammas"
## [7] "tauBs"     "alphas"

Bản tóm tắt nội dung mô hình Joint model cung cấp các thông tin bao gồm:

Thống kê về tính phù hợp của mô hình, sử dụng Log pseudo marginal likelihood value (LPML), deviance information criterion (DIC) và kích thước của thành phần tham số trong DIC (pD).

Trung bình của phân phối hậu nghiệm của tất cả các tham số hồi quy, kèm theo sai số chuẩn (SE) và khoảng tin cậy 95%. Những kết quả này được ước tính dựa vào chuỗi MCMC (được khảo sát như một time series vector). Một điểm thú vị đó là tuy sử dụng phương pháp ước lượng tham số hồi quy theo Bayes, với phân phối hậu nghiệm, nhưng package JMbayes lại trình bày kết quả suy diễn thống kê dưới hình thức p-value (tail probability với null hypothesis testing như trong phương pháp frequentist): p value ở đây chính là xác suất tham số hồi quy lớn hơn hoặc nhỏ hơn zero.

P value = 2 * min { Pr(b>0),Pr(b<0)} với b là tham số hồi quy cần kiểm tra

Thực ra cách suy diễn này cũng tương đồng như phương pháp của John Kruschke dựa vào 1 ngưỡng ý nghĩa hay 1 khoảng vô nghĩa thực dụng (ROPE), chỉ khác là Kruschke dùng mật độ xác suất (bao nhiêu % phân phối hậu nghiệm nằm trong/ngoài ROPE hay nằm bên trái/phải ngưỡng ý nghĩa) còn Rizopoulos dùng xác suất p.

Phân phối hậu nghiệm của tham số hồi quy được lưu trữ trong các chuỗi MCMC, có tất cả 8 chuỗi tương ứng với 8 thành phần trong joint model. Ở đây ta chỉ khảo sát những thông tin quan trọng nhất gồm có:

Chuỗi MCMC Beta

Chuỗi MCMC beta chứa tham số hồi quy cho thành phần Longitudinal (mô hình random slope ước tính log(prothrombin) tại thời điểm t cho một cá thể nhất định).

jointfit$mcmc$betas%>%head()%>%knitr::kable()

(Intercept)	ns(time, 3)1	ns(time, 3)2	ns(time, 3)3
4.166755	-0.1371939	0.2371638	0.4604426
4.232075	-0.0993499	0.2963757	0.3017229
4.198056	-0.0616516	0.2462436	0.3165887
4.219529	-0.1473808	0.2748667	0.2276635
4.234687	-0.1174135	0.2467452	0.3722705
4.195657	-0.0815266	0.1961530	0.3373747

betadf=jointfit$mcmc$betas%>%as_tibble()

names(betadf)=c("Intercept","NS_Time1","NS_Time2","NS_Time3")

betadf%>%gather(Intercept:NS_Time3,key="Components",value="Coefficient")%>%
  ggplot()+
  geom_density(aes(x=Coefficient,fill=Components),alpha=0.6)+
  theme_bw()+
  facet_wrap(~Components,scales="free",ncol=2)

betadf%>%
  mutate(.,Iter=as.numeric(rownames(.)))%>%
  gather(Intercept:NS_Time3,key="Components",value="Coefficient")%>%
  ggplot()+
  geom_path(aes(x=Iter,y=Coefficient,color=Components,group=1),alpha=0.7)+
  theme_bw()+
  facet_wrap(~Components,scales="free",ncol=1)

Chuỗi MCMC sigma

Chuỗi MCMC sigma chứa phân phối hậu nghiệm của Residual trong mô hình random slope:

jointfit$mcmc$sigma%>%head()

##          sigma
## [1,] 0.2464715
## [2,] 0.2351639
## [3,] 0.2356347
## [4,] 0.2416009
## [5,] 0.2392787
## [6,] 0.2360419

sigmadf=jointfit$mcmc$sigma%>%as_tibble()

p1=sigmadf%>%ggplot()+
  geom_density(aes(x=sigma),alpha=0.6,fill="red")+
  theme_bw()

p2=sigmadf%>%
  mutate(.,Iter=as.numeric(rownames(.)))%>%
  ggplot()+
  geom_path(aes(x=Iter,y=sigma),alpha=0.7,color="red")+
  theme_bw()

gridExtra::grid.arrange(p1,p2,ncol=1)

Chuỗi MCMC b

Chuỗi b chứa hiệu ứng ngẫu nhiên cá thể cho Intercept và Slope, đây là 1 matrix rất lớn, với kích thước = số lượt lấy mẫu (khoảng 8000) x số tham số trong mô hình x số bệnh nhân.

Chuỗi gammas

Chuỗi gamma chứa phân phối hậu nghiệm của hiệu ứng riêng phần của yếu tố điều trị trong mô hình Survival (treatprednisone).

gammadf=jointfit$mcmc$gammas%>%as_tibble()

p1=gammadf%>%ggplot()+
  geom_density(aes(x=treatprednisone),alpha=0.6,fill="purple")+
  theme_bw()

p2=gammadf%>%
  mutate(.,Iter=as.numeric(rownames(.)))%>%
  ggplot()+
  geom_path(aes(x=Iter,y=treatprednisone),alpha=0.7,color="purple")+
  theme_bw()

gridExtra::grid.arrange(p1,p2,ncol=1)

Chuỗi BS.gamma

Chuỗi BS.gamma chứa phân phối hậu nghiệm cho tham số hồi quy của 17 phân đoạn thời gian trong mô hình Survival, do kích thước của nó quá lớn nên chúng ta cũng không khảo sát sâu hơn

bsgam=jointfit$mcmc$Bs.gammas%>%as_tibble()

bsgam%>%head()%>%knitr::kable()

Bs.gammas1	Bs.gammas2	Bs.gammas3	Bs.gammas4	Bs.gammas5	Bs.gammas6	Bs.gammas7	Bs.gammas8	Bs.gammas9	Bs.gammas10	Bs.gammas11	Bs.gammas12	Bs.gammas13	Bs.gammas14	Bs.gammas15	Bs.gammas16	Bs.gammas17
8.142618	8.028174	7.852684	7.735790	7.625893	7.529252	7.512875	7.545106	7.594168	7.777030	7.977989	8.043824	8.146587	8.484370	8.790008	9.097839	9.537896
8.046355	7.934682	7.774667	7.744428	7.659609	7.573838	7.497213	7.470633	7.551092	7.768141	7.976529	8.085545	8.260737	8.603179	8.977115	9.376499	9.913058
8.286810	8.126373	7.938846	7.903592	7.838276	7.818553	7.752809	7.698992	7.779366	7.979521	8.203416	8.300705	8.451332	8.760036	9.113401	9.450458	9.935655
8.442114	8.085402	7.845935	7.738835	7.641463	7.608680	7.572032	7.568777	7.690309	7.932411	8.226173	8.405535	8.524631	8.795154	9.110450	9.388181	9.702182
8.558670	8.198795	7.973443	7.842892	7.697231	7.637369	7.614329	7.616467	7.749046	7.997117	8.236845	8.400628	8.536905	8.839088	9.211345	9.564557	9.933011
8.811566	8.385017	8.160656	8.021760	7.886372	7.867034	7.842348	7.831169	7.895563	8.041047	8.227116	8.363656	8.448777	8.754916	9.155494	9.567143	9.987988

Hmisc::describe(bsgam)

## bsgam 
## 
##  17  Variables      2000  Observations
## ---------------------------------------------------------------------------
## Bs.gammas1 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.934   0.7412    7.821    8.046 
##      .25      .50      .75      .90      .95 
##    8.504    8.952    9.379    9.751   10.034 
## 
## lowest :  6.937464  7.016839  7.091582  7.147383  7.162695
## highest: 10.690605 10.692979 10.722805 10.752265 10.806587
## ---------------------------------------------------------------------------
## Bs.gammas2 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.801   0.7434    7.711    7.929 
##      .25      .50      .75      .90      .95 
##    8.342    8.814    9.251    9.621    9.868 
## 
## lowest :  6.844597  7.015248  7.037638  7.053279  7.132792
## highest: 10.479946 10.501547 10.566833 10.589966 10.685109
## ---------------------------------------------------------------------------
## Bs.gammas3 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.668   0.7611    7.559    7.791 
##      .25      .50      .75      .90      .95 
##    8.188    8.679    9.150    9.509    9.730 
## 
## lowest :  6.668485  6.928073  6.942118  6.955891  6.984306
## highest: 10.447261 10.477611 10.573543 10.575178 10.668534
## ---------------------------------------------------------------------------
## Bs.gammas4 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.542    0.772    7.406    7.661 
##      .25      .50      .75      .90      .95 
##    8.052    8.559    9.034    9.387    9.624 
## 
## lowest :  6.500742  6.789214  6.816522  6.827519  6.849389
## highest: 10.352992 10.376645 10.450199 10.467366 10.564317
## ---------------------------------------------------------------------------
## Bs.gammas5 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.453   0.7784    7.319    7.569 
##      .25      .50      .75      .90      .95 
##    7.956    8.466    8.949    9.328    9.558 
## 
## lowest :  6.428777  6.592291  6.639888  6.674039  6.678085
## highest: 10.315809 10.336024 10.343200 10.353307 10.386162
## ---------------------------------------------------------------------------
## Bs.gammas6 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.399    0.791    7.245    7.505 
##      .25      .50      .75      .90      .95 
##    7.893    8.409    8.892    9.292    9.538 
## 
## lowest :  6.392566  6.500920  6.512691  6.587005  6.600791
## highest: 10.280576 10.298094 10.303442 10.340657 10.369492
## ---------------------------------------------------------------------------
## Bs.gammas7 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.392   0.7961    7.258    7.513 
##      .25      .50      .75      .90      .95 
##    7.886    8.406    8.881    9.295    9.538 
## 
## lowest :  6.323360  6.391773  6.478027  6.515396  6.536359
## highest: 10.282893 10.291667 10.352427 10.396060 10.402068
## ---------------------------------------------------------------------------
## Bs.gammas8 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.419   0.7986    7.292    7.530 
##      .25      .50      .75      .90      .95 
##    7.913    8.437    8.922    9.319    9.551 
## 
## lowest :  6.345751  6.398935  6.472127  6.529794  6.545229
## highest: 10.238440 10.308912 10.338838 10.384987 10.418001
## ---------------------------------------------------------------------------
## Bs.gammas9 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.491   0.8058    7.367    7.601 
##      .25      .50      .75      .90      .95 
##    7.972    8.504    9.003    9.404    9.620 
## 
## lowest :  6.411109  6.467243  6.538231  6.622581  6.638399
## highest: 10.306940 10.385985 10.401643 10.404050 10.423646
## ---------------------------------------------------------------------------
## Bs.gammas10 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.606   0.8119    7.468    7.714 
##      .25      .50      .75      .90      .95 
##    8.081    8.615    9.118    9.523    9.757 
## 
## lowest :  6.593275  6.601119  6.695314  6.704427  6.713113
## highest: 10.493302 10.497455 10.506366 10.518864 10.524146
## ---------------------------------------------------------------------------
## Bs.gammas11 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.757   0.8146    7.607    7.845 
##      .25      .50      .75      .90      .95 
##    8.221    8.769    9.264    9.661    9.915 
## 
## lowest :  6.786049  6.821886  6.875645  6.879372  6.881041
## highest: 10.624128 10.655827 10.661483 10.682578 10.693685
## ---------------------------------------------------------------------------
## Bs.gammas12 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    8.923    0.819    7.747    7.993 
##      .25      .50      .75      .90      .95 
##    8.401    8.938    9.442    9.825   10.074 
## 
## lowest :  6.946539  6.993020  7.022604  7.059983  7.087073
## highest: 10.775010 10.796138 10.796169 10.853785 10.903102
## ---------------------------------------------------------------------------
## Bs.gammas13 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1      9.1   0.8276    7.917    8.163 
##      .25      .50      .75      .90      .95 
##    8.569    9.102    9.636   10.002   10.270 
## 
## lowest :  7.041233  7.207530  7.263250  7.283345  7.293830
## highest: 10.988641 10.994570 11.000694 11.071627 11.116787
## ---------------------------------------------------------------------------
## Bs.gammas14 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1     9.29   0.8472    8.099    8.314 
##      .25      .50      .75      .90      .95 
##    8.734    9.292    9.843   10.206   10.477 
## 
## lowest :  7.030630  7.439889  7.465634  7.476441  7.490832
## highest: 11.288478 11.307226 11.312833 11.349971 11.448798
## ---------------------------------------------------------------------------
## Bs.gammas15 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    9.471   0.8698    8.244    8.463 
##      .25      .50      .75      .90      .95 
##    8.933    9.464   10.032   10.426   10.701 
## 
## lowest :  7.046545  7.532705  7.553363  7.553397  7.566575
## highest: 11.535976 11.538103 11.551505 11.608862 11.627421
## ---------------------------------------------------------------------------
## Bs.gammas16 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1     9.65   0.9065    8.327    8.606 
##      .25      .50      .75      .90      .95 
##    9.071    9.644   10.222   10.655   10.944 
## 
## lowest :  7.138890  7.624017  7.654898  7.672475  7.705988
## highest: 11.793102 11.802715 11.813050 11.877193 11.878617
## ---------------------------------------------------------------------------
## Bs.gammas17 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     2000        0     1982        1    9.822   0.9768    8.423    8.683 
##      .25      .50      .75      .90      .95 
##    9.184    9.840   10.417   10.919   11.209 
## 
## lowest :  7.196583  7.718103  7.722063  7.767341  7.772541
## highest: 12.137313 12.161188 12.185025 12.192673 12.310141
## ---------------------------------------------------------------------------

Chuỗi Alphas

Chuỗi alphas chứa thông tin quan trọng nhất trong mô hình joint model mà ta quan tâm, đó là hệu ứng tương tác giữa bệnh sử diễn tiến của prothrombin và phân nhóm điều trị (association effect).

jointfit$mcmc$alphas%>%head()

##         Assoct
## [1,] -2.253908
## [2,] -2.239581
## [3,] -2.306681
## [4,] -2.292939
## [5,] -2.328565
## [6,] -2.336292

alphadf=jointfit$mcmc$alphas%>%as_tibble()

p1=alphadf%>%ggplot()+
  geom_density(aes(x=Assoct),alpha=0.6,fill="gold")+
  theme_bw()

p2=alphadf%>%
  mutate(.,Iter=as.numeric(rownames(.)))%>%
  ggplot()+
  geom_path(aes(x=Iter,y=Assoct),alpha=0.7,color="orange")+
  theme_bw()

gridExtra::grid.arrange(p1,p2,ncol=1)

5 Ứng dụng của mô hình

Tiên lượng cho cá thể

Chúng ta dùng mô hình này để tiên lượng cho một bệnh nhân trong mẫu khảo sát, thí dụ số 9.

ND =prothro[prothro$id == 9, ]

longfit=predict(jointfit, 
                newdata = ND, 
                type = "Subject",
                interval = "confidence", 
                return = TRUE)

survfit=survfitJM(jointfit,newdata = ND)

survfit[1]$summaries$`9`%>%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()+
  geom_vline(xintercept=max(longfit$Time),color="blue",linetype=2)+
  labs(title="Patient N°9",x="Time",y="Predicted survival probability",caption="Joint Model")

plot(survfit, estimator = "mean", include.y = TRUE,
     conf.int = TRUE, fill.area = TRUE, col.area = "gold")

longfit%>%
  ggplot(aes(x=time))+
  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()+
  geom_vline(xintercept=max(longfit$Time),color="blue",linetype=2)+
  labs(title="Patient N°9",x="Time",y="Predicted prothromb (Log scale)",caption="Joint Model")

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

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

6 Kết luận

Joint model là một phương pháp hữu ích, và package JMbayes là một cách tiếp cận thú vị, với phương pháp Bayes (MCMC) lai với null hypothesis testing và log pseudo likelihood, nhằm ước tính hệ số hồi quy cho Joint model.

Nếu bạn tìm hiểu sâu hơn nội dung package, sẽ thấy khả năng của nó rộng hơn rất nhiều so với phương pháp REML trong bài trước, và thí dụ cơ bản mà ta vừa thực hiện.

Chúc các bạn thành công.

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