真實世界大數據分析系列聯(lián)合模型的R語言實現與解讀——以分析PBC2數據

真實世界?數據分析系列聯(lián)合模型的R語?實現與解讀——以分析PBC2數據集

為例

公眾號：醫(yī)學?數據挖掘分析

本系列上?期我們已經簡單介紹了聯(lián)合模型的相關原理（往期推?：真實世界?數據分析系列|聯(lián)合模型的簡單介紹），

本期我們從實例出發(fā)，看如何使?R語?實現聯(lián)合模型，如何解讀聯(lián)合模型的結果。

數據來源

數據為JM程序包中的1974年到1984年Mayo原發(fā)性膽汁肝硬化（Mayo Clinic Primary Biliary Cirrhosis Data, PBC2）數

據集，包含了312例患者，總計1945個觀測值，包含患者編號、死亡或刪失時間、?存狀態(tài)、?藥情況、年齡、性別、

觀測時間、?清膽紅素等20個變量。其中140例患者死亡，29例接受移植?術，143例患者存活。由于死亡和刪失，該

研究對患者平均進?6.2次隨訪，中位?存時間為9.43年。具體數據如下圖所?。

本研究的?的為：探討縱向指標rBilir與肝硬化患者死亡風險之間的關系。下?構建聯(lián)合模型。

?、構建聯(lián)合模型

library(JMbayes)

pbc2$status2 <- c(pbc2$status != "alive")

$status2 <- c($status != "alive")

（1）構建縱向?模型：

為了滿?擬合線性混合模型的因變量分布（服從正態(tài)分布）要求，對?清膽紅素進?對數變換得到log(rBilir)，并以

log(rBilir)為因變量，選擇觀測時間（year）、?藥情況（drug）、觀測時間和?藥情況的交互作?為固定效應部分，

患者標識影響下的觀測時間作為隨機效應，構建線性混合模型：

1 <- lme(log(rBilir) ~ drug * year, data = pbc2,

random = ~ year | id, method = "ML")

summary(1)

（2）構建?存?模型：

以?藥情況、年齡、?藥情況和年齡的交互作?為基線變量，構建cox模型：

1 <- coxph(Surv(years, status2) ~ drug * age, data = , x = TRUE)

summary(1)

（3）聯(lián)合模型的參數估計：

構建聯(lián)合模型：

1 <- jointModelBayes(1, 1, timeVar = "year")

summary(1)

我們看聯(lián)合模型的結果，縱向?模型部分，觀測時間對log(rBilir)的變化差異有統(tǒng)計學意義（P<0.001），隨著觀測時

間的增加，log(rBilir)也逐漸增加（β=0.2313，P<0.001）；?存?模型部分，患者年齡差異有統(tǒng)計學意義

（P<0.001），患者年齡每增加?個單位，患者死亡風險增加1.0358倍（95% CI: 1.0089-1.0618）。另外，聯(lián)合模型的

共享參數（Assoct）有統(tǒng)計學意義，表明log(rBilir)每縱向增加1個單位，患者發(fā)?死亡風險增加3.8675倍（95% CI:

3.2619-4.6856）。

?、模型預測效能與評價

我們可以計算出模型的AIC、BIC、LogLikelihood值，分別為4761.51、4855.935、-2350.755。這些評價指標可?于模

型間的對?，我們挑選值AIC/BIC最?或對數似然函數值最?的模型為最佳模型。

利?到起始時間點的可?縱向信息，基于聯(lián)合模型可以計算指定時間點的ROC曲線和AUC估計值。下圖為起始時間為

5，預測到時間點為13的ROC曲線及其對應的AUC值，AUC值越接近1預測效果越好。

plot(rocJM(1, pbc2, Tstart = 5, Thoriz = 8))

aucJM(1, pbc2, Tstart = 5, Thoriz = 8)

三、個體化動態(tài)預測

根據建?的聯(lián)合模型，我們可以對任?患者的?存概率和縱向結果進?動態(tài)預測，如對2號患者進?預測：

ND <- pbc2[pbc2$id == 2, ] # the data of Subject 2

survPreds <- vector("list", nrow(ND))

for (i in 1:nrow(ND)) {

survPreds[[i]] <- survfitJM(1, newdata = ND[1:i, ])

}

par(mfrow = c(2, 2), oma = c(0, 2, 0, 2))

for (i in c(1,3,5,7)) {

plot(survPreds[[i]], estimator = "median", = TRUE,

include.y = TRUE, main = paste("Follow-up time:",

round(survPreds[[i]]$, 1)), ylab = "", ylab2 = "")

}

mtext("log rum bilirubin", side = 2, line = -1, outer = TRUE)

mtext("Survival Probability", side = 4, line = -1, outer = TRUE)

上圖為2號患者在隨訪期間的動態(tài)?存概率。垂直虛線表?最后?次?清膽紅素測定的時間點。垂線左邊描述了擬合的

縱向軌跡，垂線右邊的實線為?存概率中位數估計量，虛線對應其95%置信區(qū)間。我們觀測到，在第三次測量后，2號

患者的對數?清膽紅素?平明顯增加，?存概率的下降速率變得更為陡峭，表明她的病情正在惡化。

參考?獻：

[1] 翟映紅, 陳琪, 韓賀東, 趙欣欣, ?永晴, 周詳, 賀佳. 聯(lián)合模型介紹及在醫(yī)學研究中的應?. 中華流?病學雜志, 2019,

40(11): 1456-1460.

[2] 肖宇飛,馮佳寧,王曉璇,?倩,?福艷,王素珍, XIAO Yufei, FENG Jianing, WANG Xiaoxuan, MAO Qian, SHI Fuyan,

WANG Suzhen. 利?數據庫數據采?聯(lián)合模型動態(tài)預測312例肝硬化患者預后的觀察分析[J]. ?東?學學報：醫(yī)學版,

2020, 58(9):7.

[3] Dimitris, Rizopoulos, Laura, et al. Combining Dynamic Predictions From Joint Models for Longitudinal and Time-to-

Event Data Using Bayesian Model Averaging[J]. Journal of the American Statistical Association, 2014, 109(508):1385-

1397.

[4] Rizopoulos D. Joint models for longitudinal and time-to-event data:with applications in R[M]. Boca Raton: CRC

Press, 2012.

[5] Kan L , Sheng L . Dynamic predictions in Bayesian functional joint models for longitudinal and time-to-event data:

An application to Alzheimer's dia[J]. Statistical Methods in Medical Rearch, 2017, 28(2):217.

更多內容，請關注“醫(yī)學?數據挖掘分析”公眾號，歡迎留?聯(lián)系～

郵箱：*********************

地址：?州市天河區(qū)珠江東路?德置地秋?場F座