真实世界大数据分析系列联合模型的R语言实现与解读——以分析PBC2数据

真实世界⼤数据分析系列联合模型的R语⾔实现与解读——以分析PBC2数据集

为例

公众号：医学⼤数据挖掘分析

本系列上⼀期我们已经简单介绍了联合模型的相关原理（往期推⽂：真实世界⼤数据分析系列|联合模型的简单介绍），

本期我们从实例出发，看如何使⽤R语⾔实现联合模型，如何解读联合模型的结果。

数据来源

数据为JM程序包中的1974年到1984年Mayo原发性胆汁肝硬化（Mayo Clinic Primary Biliary Cirrhosis Data, PBC2）数

据集，包含了312例患者，总计1945个观测值，包含患者编号、死亡或删失时间、⽣存状态、⽤药情况、年龄、性别、

观测时间、⾎清胆红素等20个变量。其中140例患者死亡，29例接受移植⼿术，143例患者存活。由于死亡和删失，该

研究对患者平均进⾏6.2次随访，中位⽣存时间为9.43年。具体数据如下图所⽰。

本研究的⽬的为：探讨纵向指标rBilir与肝硬化患者死亡风险之间的关系。下⾯构建联合模型。

⼀、构建联合模型

library(JMbayes)

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

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

（1）构建纵向⼦模型：

为了满⾜拟合线性混合模型的因变量分布（服从正态分布）要求，对⾎清胆红素进⾏对数变换得到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）联合模型的参数估计：

构建联合模型：

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

summary(1)

我们看联合模型的结果，纵向⼦模型部分，观测时间对log(rBilir)的变化差异有统计学意义（P<0.001），随着观测时

间的增加，log(rBilir)也逐渐增加（β=0.2313，P<0.001）；⽣存⼦模型部分，患者年龄差异有统计学意义

（P<0.001），患者年龄每增加⼀个单位，患者死亡风险增加1.0358倍（95% CI: 1.0089-1.0618）。另外，联合模型的

共享参数（Assoct）有统计学意义，表明log(rBilir)每纵向增加1个单位，患者发⽣死亡风险增加3.8675倍（95% CI:

3.2619-4.6856）。

⼆、模型预测效能与评价

我们可以计算出模型的AIC、BIC、LogLikelihood值，分别为4761.51、4855.935、-2350.755。这些评价指标可⽤于模

型间的对⽐，我们挑选值AIC/BIC最⼩或对数似然函数值最⼤的模型为最佳模型。

利⽤到起始时间点的可⽤纵向信息，基于联合模型可以计算指定时间点的ROC曲线和AUC估计值。下图为起始时间为

5，预测到时间点为13的ROC曲线及其对应的AUC值，AUC值越接近1预测效果越好。

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

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

三、个体化动态预测

根据建⽴的联合模型，我们可以对任⼀患者的⽣存概率和纵向结果进⾏动态预测，如对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号患者在随访期间的动态⽣存概率。垂直虚线表⽰最后⼀次⾎清胆红素测定的时间点。垂线左边描述了拟合的

纵向轨迹，垂线右边的实线为⽣存概率中位数估计量，虚线对应其95%置信区间。我们观测到，在第三次测量后，2号

患者的对数⾎清胆红素⽔平明显增加，⽣存概率的下降速率变得更为陡峭，表明她的病情正在恶化。

参考⽂献：

[1] 翟映红, 陈琪, 韩贺东, 赵欣欣, ⾼永晴, 周详, 贺佳. 联合模型介绍及在医学研究中的应⽤. 中华流⾏病学杂志, 2019,

40(11): 1456-1460.

[2] 肖宇飞,冯佳宁,王晓璇,⽑倩,⽯福艳,王素珍, XIAO Yufei, FENG Jianing, WANG Xiaoxuan, MAO Qian, SHI Fuyan,

WANG Suzhen. 利⽤数据库数据采⽤联合模型动态预测312例肝硬化患者预后的观察分析[J]. ⼭东⼤学学报：医学版,

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.

更多内容，请关注“医学⼤数据挖掘分析”公众号，欢迎留⾔联系～

邮箱：*********************

地址：⼴州市天河区珠江东路⾼德置地秋⼴场F座