系列文章目录

R语言–Cox模型校准曲线原理（一）数据

文章目录

系列文章目录
前言
一、回顾校准曲线绘图数据
二、计算
- 1.分组
- 2.mean.predicted、n、events
- 3.KM
- 4.std.err
总结

前言

上一节给自己留了个坑，这赶紧过来填上。众所周知Cox模型的校准曲线数据是rms包中的calibrate函数计算而来的。那里面的数据到底是如何计算来的？搞清楚了这个，或许对进一步了解cox模型有一定帮助。

一、回顾校准曲线绘图数据

第二排数据的来源，也就是mean.predicted，KM，std.err是如何计算的呢？这几个数据是不需要迭代就能算出来的。剩余的数据和1000次bootstrap（自助法抽样）后计算相关。

KM： 生存分析中的Kaplan Meier法简称KM法。calibrate函数计算的结果（以下简称cal）中的KM数据得就是通过Kaplan Meier法计算得到的实际*年生存率
mean.predicted： 就是通过模型计算而来的平均*年生存率。
std.err： 是预测生存率的标准误。

二、计算

为了更清楚解决这个问题，采用R语言中lung数据框中的数据，方便大家复现以及发现问题。

1.分组

在通过KM法构建COX模型的校准曲线时，会根据病例的风险进行分组，然后计算每组的平均值，关于这里为什么分组，是因为，KM法计算生存概率时对病例数有一定的要求，通俗来说就是，如果只有一两个病人，是无法计算他的*年生存率的，因此分组数不能太多导致各组病例数较少，大部分文献中分3组至4组就差不多了。

这里先弄清楚两个问题：

根据什么分组： 根据*年预测生存概率
怎么计算*年生存概率预测值： 山人自有妙计

##随便构建一个cox模型并绘制一个默认的校准曲线
set.seed(1234)
ind <- sample(1:2,size = nrow(a),prob = c(0.7,0.3),replace =T)
train <- a[ind==1,]
test <- a[ind==2,]
dd=datadist(train)
options(datadist="dd")
#设置时间为240，即240天
f<-cph(Surv(time,status)~sex+ph.ecog,time.inc=240,data = train,x=TRUE,y=TRUE,surv=TRUE)
#下面为默认校准曲线代码
cal <- calibrate(f, cmethod="KM", method="boot", u=240, m=  round(nrow(train)/4), B=1000)
plot(cal,xlim=c(0,1),ylim=c(0,1))

#####接下来根据cox模型手动分组，拟划分为4组
g=4
####训练集240天生存概率计算：（可以同理预测测试集的生存概率）
survival <- f$surv.summary[2,1,1]^exp(f$linear.predictors)
###分组（根据预测的生存概率大小）
cuts <- unique(quantile(c(0, 1, survival), seq(0, 1, length = g + 1), na.rm = TRUE))
###通过KM法计算实际生存概率（这里指分组后各组的20年生存概率，不可能计算个体的20年生存率，只有群体才有这个性质）
real_surv <- groupkm(survival,Srv = f$y,u = 240,cuts = cuts)
##对比结果和上述代码的校准曲线cal的7-10列数据。
real_surv
cal

可以发现real_surv结果中的x为mean.predicted，n为每一个分组中病例数，events为每一个分组中发生结局事件的例数。

那么问题又来了，不是分4组吗，为什么这里第一组是中的实际生存概率以及预测生存概率是NA，病例数为0呢。主要原因在于：由于变量数较少导致预测的概率分布较集中（比如例子中，sex为2分类，ph.ecog为2分类，那么预测出来的概率共为4中情况）。内部分组函数使用cut2。

cuts
cut2(survival,cuts)
table(cut2(survival,cuts))

2.mean.predicted、n、events

这个值简而言之就是*年生存率预测值的平均值。这里再深究以下它的计算过程。
由于第1组（共4组）无病例数，因此这里以第2组为例子。

#上一部分代码已经给出分组情况
cut2(survival,cuts)
group=unclass(cut2(survival,cuts))
#可以发现结果就是real_surv数据中的”x“的第二个数据，即real_surv[2,'x']
mean(survival[group=="2",])
real_surv[2,'x']##n
table(f$y[group=="2",2])

3.KM

real_surv中的KM值是根据Kaplan Meier法计算的实际生存概率：
survfitKM函数会根据生存时间、缺失、结局事件发生计算生存概率，那么我们只需要找到结局时间240对应的生存概率即可，生存时间储存在km$time中，生存概率储存在km$surv中。

#通过KM法，240天生存率，详见相关统计教材
km = survfitKM(as.factor(as.factor(rep("1", sum(group=='2')))), f$y)
##
kmp = km$surv[max((1:length(c(km$time)))[c(km$time) <= 240])]
kmp
real_surv
#结果一致

4.std.err

关于这个标准误的计算，和上述km中实际生存概率一样，只需要找到结局时间240对应的std.err即可，生存时间储存在km$time中，标准误储存在km$std.err中。

#通过KM法，240天生存率，详见相关统计教材
kmstderr = km$std.err[max((1:length(c(km$time)))[c(km$time) <= 240])]
kmstderr
real_surv
#结果一致

总结

文章纯属心血来潮，加深自己对生存分析的理解，留个坑，下一次说一下怎样通过多次的bootstrap法对计算cox模型的校准值的。

本人小白，如有错误，请批评指正。

R语言--Cox模型校准曲线原理（二）算法相关推荐

R语言cox回归模型案例（绘制列线图、校正曲线）：放疗是否会延长胰脏癌手术患者的生存时间
R语言cox回归模型案例(绘制列线图.校正曲线):放疗是否会延长胰脏癌手术患者的生存时间目录
R语言结构方程模型（SEM）在生态学领域中的应用
前言:结构方程模型(Sructural Equation Model)是一种建立.估计和检验研究系统中多变量间因果关系的模型方法,它可以替代多元回归.因子分析.协方差分析等方法,利用图形化模型方式清晰 ...
R语言BIOMOD2模型的物种分布模拟
随着生物多样性全球大会的举办,不论是管理机构及科研单位.高校都在积极准备,根据国家林草局最新工作指示,我国将积极整合.优化自然保护地,加快推进国家公园体制试点,构建以国家公园为主体的自然保护地体系.针 ...
R语言结构方程模型（SEM）在生态学领域中的实践
结构方程模型(Sructural Equation Model)是一种建立.估计和检验研究系统中多变量间因果关系的模型方法,它可以替代多元回归.因子分析.协方差分析等方法,利用图形化模型方式清晰展示研 ...
R语言︱机器学习模型评价指标+（转）模型出错的四大原因及如何纠错
笔者寄语:机器学习中交叉验证的方式是主要的模型评价方法,交叉验证中用到了哪些指标呢? 交叉验证将数据分为训练数据集.测试数据集,然后通过训练数据集进行训练,通过测试数据集进行测试,验证集进行验证. 模 ...
R语言结构方程模型（SEM）教程
详情点击链接:R语言结构方程模型(SEM)在生态学应用结构方程模型(Sructural Equation Model)是一种建立.估计和检验研究系统中多变量间因果关系的模型方法,它可以替代多元回归. ...
R语言分类模型：逻辑回归模型LR、决策树DT、推理决策树CDT、随机森林RF、支持向量机SVM、Rattle可视化界面数据挖掘、分类模型评估指标（准确度、敏感度、特异度、PPV、NPV）
R语言分类模型:逻辑回归模型LR.决策树DT.推理决策树CDT.随机森林RF.支持向量机SVM.Rattle可视化界面数据挖掘.分类模型评估指标(准确度.敏感度.特异度.PPV.NPV) 目录
R语言回归模型构建、回归模型基本假设（正态性、线性、独立性、方差齐性）、回归模型诊断、car包诊断回归模型、特殊观察样本分析、数据变换、模型比较、特征筛选、交叉验证、预测变量相对重要度
R语言回归模型构建.回归模型基本假设(正态性.线性.独立性.方差齐性).回归模型诊断.car包诊断回归模型.特殊观察样本分析.数据变换.模型比较.特征筛选.交叉验证.预测变量相对重要度目录
R语言KNN模型数据分类实战
R语言KNN模型数据分类实战目录 R语言KNN模型数据分类实战 #KNN基本原理 #案例分析 #数据集说明
R语言xgboost模型构建:基于prima糖尿病数据集
R语言xgboost模型构建:基于prima糖尿病数据集目录 R语言xgboost模型构建:基于prima糖尿病数据集

R语言--Cox模型校准曲线原理（二）算法