R数据分析:生存数据的预测模型建立方法与评价

之前写了生存分析列线图的做法,列线图作为一个预测模型可视化工具,我们使用它的过程其实就是一个给新数据做预测的过程,其内在本身的模型就是我们基于现有数据训练的一个预测模型,今天也算是接着上一篇文章继续写生存分析的预测模型的效果评价。

生存数据预测模型和我们之前写的连续变量结局和分类结局的预测模型不同的地方就在于我们得考虑生存数据的删失和时间因素,通过这么一个预测模型,我们期望的目标是帮助临床医生去回答特定病人在某个时间的生存概率。从这个角度讲我们对模型评估和评价的标准就有一个锚定了。

Thus, survival prediction models differ from traditional prediction models for continuous or binary outcomes by appropriately accommodating censoring that is present in time-to-event data.to answer questions such as “What is the probability that this patient will be alive in 5 years, given their baseline covariate information?” This predicted probability can then be used by clinicians to make important decisions regarding patient care

比如说我收集了某个癌症病人的很大的有代表性的数据集,我通过我的数据学习出来了一个预测模型,再来一个新的癌症病人,模型可以告诉我,这个病人能活多久。

如果这个新数据本身有标签,我们通过对比实际标签(具体时间的生存概率)和模型预测结果(预测的具体时间的生存概率),就可以来评价模型优劣。逻辑和常规(分类结局和连续结局)的预测模型还是一样的。

先回顾生存分析

依然是先回顾下生存分析中的常见的术语:

[图片上传失败...(image-690b11-1681819253492)]

我们的结局变量两个水平,一个是发生事件,另外一个是删失;同时这个结局还依赖于一个时间变量。

刚刚写了,我们做生存数据的预测,回答的是we are predicting the probability that an event happens at a particular time.----某个时刻发生事件的概率。所以这个时候常规的评价模型的指标都不好使了。

Due to the presence of the censoring in survival data, the standard evaluation metrics for regression such as root of mean squared error and ܴ R2 are not suitable for measuring the performance in survival analysis.

对于生存数据的预测模型,此时的评估模型的指标有下面3个:Concordance index (C-index),Brier score,Mean absolute error。今天的任务就是一个一个带大家捋一遍,希望能帮助大家理解在“个体特定时间的生存概率这个锚定标准下,这些指标为什么可以用来评价模型。”

Concordance index (C-index)

首先看C指数,这个一致性指数在分类结局的预测模型中给大家提到过,就是ROC曲线下面积,对于生存数据的预测模型来讲,这个指数和灵敏度特异度就没关系了,它比的是实际值和预测值的排序一致不一致。理解方法可以参考秩和检验。

For a binary outcome, C-index is identical to the area under the ROC curve (AUC).

The concordance index or C-index is a generalization of the area under the ROC curve (AUC) that can take into account censored data. It represents the global assessment of the model discrimination power.

其逻辑在于:每个个案通过模型都给它一个风险分,如果模型表现好,那么风险分高的个案应该会先发生事件,按照这个逻辑,然后我们用模型给每个个案都赋一个风险分,形成很多个可以对比的组(2个为一组):在组内确实满足刚刚讲的“风险分越大,事件越先发生”那么这个组就是一致性的组,否则就是不一致的组,这样一致性的组占所有对比组的比例就是C-index:

指数的计算方法如下:

[图片上传失败...(image-1a262b-1681819253492)]

其中,分子上是一致性的组,分母是所有组。那么这个值就是越大越好。

[图片上传失败...(image-e81d88-1681819253492)]

上面就是生存分析预测模型评价中C-index的内在逻辑,大家作为应用型科研工作者关注逻辑就好,请自动忽略掉数学表达。

Brier score

再来看第二个评价指标,叫做Brier score。这个Brier score是个案在t时间的生存状态减去t时间的预测生存概率的差的平方的均值

其可以用来评价模型的逻辑在于:如果我的模型真的可以很好的预测特定时刻的生存概率,那么对于某个时刻我的生存状态确实是1,那么模型应该说我此时的生存概率无限大;反之模型应该说我的生存概率无限小。

因为牵扯到具体时间,这个指标只能截一个时间点去看,其算法如下:

[图片上传失败...(image-732428-1681819253492)]

we found that BS depends on the selection of time point t. Generally, the median of the observation time is selected as the time point.

就是个案在t时间的生存状态减去t时间的预测生存概率的差,比如在t时间个案实际观测到的是死亡(取0),那么这时候模型预测的生存概率应该越小越好;t时间个案实际观测到的是存活(取1),那么这个时候模型预测的生存概率应该是越大越好;肯定是减了之后的差越小越好嘛,也就是这个Brier score越小越好,并且得小于0.25才说明这个模型好过瞎猜。但是这指标只能看某个时间点模型的预测准确性。

上面就是生存分析预测模型评价中Brier score的内在逻辑,大家作为应用型科研工作者关注逻辑就好,自动忽略掉数学表达。

Mean absolute error

MAE这个指标在连续变量结局的预测模型中也有的,指的是预测值和实际值的差的绝对值的和,在生存分析的预测模型中,就是实际生存时间和模型预测生存时间的差的绝对值的和。算法如下:

[图片上传失败...(image-c55d1a-1681819253492)]

这个指标只考虑了非删失数据,实际中就用得比较少。基本不用管。

模型评价实操

解释完指标之后我们再看实操做法,依然我们选取 JAMA Surg.的文章做参考,文章名如下:

Hyder O, Marques H, Pulitano C, et al. A Nomogram to Predict Long-term Survival After Resection for Intrahepatic Cholangiocarcinoma: An Eastern and Western Experience. JAMA Surg. 2014;149(5):432–438. doi:10.1001/jamasurg.2013.5168

文章中对模型评估的方法学介绍如下:

[图片上传失败...(image-bacd5b-1681819253492)]

可以看到,这篇文章报告了C指数,用自助抽样样本画了校准曲线,还进行了模型的验证。我们首先来看C指数的做法,文章中报告了C指数的值和置信区间:

Predictive accuracy (discrimination) of the final model was measured by calculating the Harrell C index, which was 0.692 (95% CI, 0.624-0.762).

如果你是用coxph函数跑模型,那么模型的输出结果中自动会出来C指数C指数的标准误的,如下图:

[图片上传失败...(image-76415f-1681819253492)]

比如我们就单独想要这个指数,可以直接运行下面代码:

cindex(formula, data)

要得到C指数的置信区间的话就得求助concordance.index函数,代码如下:

concordance.index(predict(c),surv.time = dt,surv.event = e,method = "noether")

输出如下,有C指数,标准误和对应置信区间上下限:

[图片上传失败...(image-1f08ae-1681819253492)]

看完了C指数的操作我们再看校准曲线的画法,论文中给到的校准曲线长这样:

[图片上传失败...(image-b15934-1681819253492)]

首先我们来理解什么是校准曲线,上图中横轴是模型预测的生存概率,纵轴是实际的生存概率。图中还有一条灰色的虚线,代表预测概率和实际的生存概率一致,最理想的情况下校准曲线是一条对角线(预测概率等于实际概率),我们实际写文章的时候只要看着不要偏太离谱就行。

Calibration plot is a visual tool to assess the agreement between predictions and observations in different percentiles (mostly deciles) of the predicted values.

还要理解的是我们本身生存概率的分布是连续的,而图中只是画了3个点,这是因为算法将数据进行了分箱处理,上图中就是将原始数据分成了3组,这个操作使用calibrate函数中的参数m进行控制

For survival models, "predicted" means predicted survival probability at a single time point, and "observed" refers to the corresponding Kaplan-Meier survival estimate, stratifying on intervals of predicted survival。

同时,对于生存数据我们本身做预测的时候也是需要限定时间的,所以需要设定参数u。

比如我们要自助抽样20次,数据分箱,每箱200个,做时间点6的校准曲线的示例代码如下:

cal <- calibrate(f, u=6, cmethod='KM', m=200, B=20)
plot(cal)

关于模型验证的结果,论文中通过报告重复抽样验证结果中训练数据和测试数据的C指数说明了模型并没有过拟合,原文如下:

Bootstrap validation of the model with 300 iterations revealed minimal evidence of model overfit. The training data set C statistic was 0.699, and the testing data set C statistic was 0.706, which represented the bias-corrected estimate of model performance in the future.

此部分的实现代码如下:

validate(f, B=300) 

通过输出结果中便可推算出相应数据集的C指数。

[图片上传失败...(image-fc62b0-1681819253492)]

D_{xy} are equal to 2 * (C - 0.5)where C is the C-index or concordance probability

然后我们再对比下训练数据和测试数据的C指数的差异就可以去得到我们自己模型的结论。

好了,到这儿本期按照JAMA surgery文章给大家写的生存数据预测模型的做法与评价方法就给大家写完了,其实生存数据的预测模型还有别的评价方法比如时间依赖的ROC,决策曲线等等,安排在下期,请持续关注。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,362评论 5 477
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,330评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,247评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,560评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,580评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,569评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,929评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,587评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,840评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,596评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,678评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,366评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,945评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,929评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,165评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 43,271评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,403评论 2 342

推荐阅读更多精彩内容