R语言实战——Cox 比例风险回归模型

简介: COX比例风险模型(cox proportional-hazards model)是英国统计学家D.R.COX于1972年提出的一种半参数回归模型,它可同时研究多个风险因素和事件结局发生情况、发生时间的关系,从而克服了简单生存分析中单因素限制的不足。

定义

COX比例风险模型(cox proportional-hazards model)是英国统计学家D.R.COX于1972年提出的一种半参数回归模型,它可同时研究多个风险因素和事件结局发生情况、发生时间的关系,从而克服了简单生存分析中单因素限制的不足。

COX回归分析特点

鉴于临床数据的特殊性,COX回归比起一般的多重线性回归Logistic回归在临床研究中具有更为广泛的应用。为了突出COX回归分析的特点,将以上三种回归分析方法作一次比较:

COX回归与另外两者的共同点:


(1)分析目的:均可用于研究自变量影响程度,校正混杂因素以及作预测分析。

(2)自变量类型:均可为连续型数值变量,或离散型分类变量,顺序变量。

(3)自变量筛选:均可采用逐步回归方法筛选变量。

COX回归与另外两者的相异点:

(1)因变量类型:COX回归因变量为生存资料(包含二分类结局变量和连续型生存时间变量),而多重线性回归因变量为数值型变量,Logistic回归因变量为分类或顺序型变量。

(2)因变量数据分布:COX回归对于数据分布不作要求,而多重线性回归,Logistic回归则要求数据分布分别近似正态分布和二项分布。

(3)模型参数解释:当自变量增大一个单位时,对于COX回归,改变了风险比HR的自然对数值。对于Logistic回归,改变了优势比OR的自然对数值。而对于多重线性回归,改变的即是Y值本身。

(4)是否允许数据删失:COX回归允许删失值,而多重线性回归,Logistic回归不允许。

常用:

COX回归分析原理简介

作为生存分析的子项目,COX回归分析的掌握有赖于对生存概率,风险概率,累积风险概率等基础知识的掌握。

(COX回归模型,又称“比例风险回归模型 (proportional hazards model,简称Cox模型)”,是由英国统计学家D.R.Cox (1972)年提出的一种半参数回归模型。 该模型以生存结局和生存时间为因变量,可同时分析众多因素对生存期的影响,能分析带有截尾生存时间的资料,且不要求估计资料的生存分布类型。)

在简单生存分析中,由于仅考虑单个影响因素(且为分类型或顺序型变量),故采取的是直接绘制生存曲线并作 Log Rank检验(https://zhuanlan.zhihu.com/p/392104512)来判断影响因素和生存概率的方法。

而在COX回归分析中,需要同时考虑多个影响因素(可为分类/顺序型变量,也可为数值型变量),故而绘制生存曲线的方式显然不合适,此时就需通过建模的方法来进一步分析。

COX回归模型推导

在多因素情况下,风险概率的计算需要同时考虑生存时间T和自变量X,因此用 h(t,x) 来表示 t 时刻的风险函数,若设定自变量取值为0,则称 h(t,0) 为 t 时刻的基准风险函数。

固定时间点t,取风险函数和基准风险函数之比可得到t时刻下的风险比值HR,该HR是关于自变量X的函数,且不再依赖于时间T,故称其为比例风险模型

若对上述模型做简单变换,即得到我们常见的COX回归模型(也称COX比例风险模型),其中前半部分为非参数模型,后半部分为参数模型,故COX回归模型为一种半参数模型。关于模型详细推导,https://www.zhihu.com/question/343779367/answer/1493205766

COX回归模型系数意义

在COX回归模型中,取某一自变量系数为e的幂数,得到的值即为HR值。考虑HR值在临床研究中的实际意义,则当系数大于0(HR>1)时,该自变量为危险因素;当系数小于0(HR<1)时,该自变量为保护因素。

COX回归模型假设条件

(1)比例风险假设:又称PH假设,指模型自变量系数为固定值,不随时间T的变化而变化。

(2)对数线性假设: 模型中对数风险比值Ln(HR)与自变量X呈线性关系。

补充:

做临床研究,学过统计的应该都知道

OR,RR,HR(统称为“3R”)

代表什么意思,不过虽然很多人知道这3个之间的区别和意义,但是在做Meta分析的时候还是会遇到下面的问题:

纳入的原始研究有的用RR,有的用OR,有的用HR,那么在Meta分析(https://baike.baidu.com/item/Meta-Analysis/811753)合并效应值的时候到底用哪个呢?

在说这个之前,我们再回顾下这3R的定义和差别:

OR

(odds ratio):

又叫比值比、优势比,指病例组中暴露人数与非暴露人数的比值除以对照组中暴露人数与非暴露人数的比值,是病例对照研究中的一个常用指标。

OR 的计算公式是:

OR=(病例组暴露人数 /非暴露人数) /(对照组暴露人数 /非暴露人数)

RR:(Risk Ratio):

又叫相对危险度、危险比,实质就是两个率的比,是两组真实发病率、患病率或死亡率的比值,是队列研究的常用指标。

RR 的计算公式是:RR=暴露组的发病或死亡率 /非暴露组的发病或死亡率

从公式即可看出OR或者RR的值越大,表明暴露的效应越大 ,暴露与结局关联的强度也就越大。以 1 为界限:OR或者RR=1,暴露因素与疾病之间无关联;OR或者RR>1,

暴露因素是疾病的危险因素(正相关),认为暴露与疾病呈 "正"关联;OR或者RR<1,

暴露因素是疾病的保护因素(负相关) ,认为暴露与疾病呈 "负"关联。

这里举个简单的例子

以吸烟与肺癌关系的研究为例,两组患者发生肺癌的数量如四格表所示:

OR

=(A/B)÷(C/D)

*A/B 指吸烟人群中发生肺癌的比例, C/D 指非吸烟人群中肺癌的比例。

RR

=[A/(A+B)]÷ [C/(C+D)]

*A/(A+B) 即吸烟组肺癌的发生率, C/(C+D) 为不吸烟组肺癌的发生率。

上面的公式差异主要是分母,即 A/(A+B) 和 A/B 的区别,如果A的数值非常小,则 A/(A+B)≈A/B。

这有什么意义呢?

就是如果疾病发生率非常低(比如一些罕见病),OR和RR的大小是近似的。因此这类疾病,

由于队列研究实施比较困难,但是病例对照研究中的OR值比较易于获得,因此我们常常通过计算OR来替代RR值。从临床解释上来说RR值比OR值解释起来更容易,比如RR值是1.67,那么你就可以说吸烟的人肺癌发生率是不吸烟的1.67倍。

那么有同学就会问为什么论文里看到OR值比较多,而RR值非常少呢?

这就是因为RR值需要得到两个率才能计算,但是基于医院的研究往往很难得到率,基于医院的临床研究中常常是只能得到分子而无法获得分母的大小(就比如上面的例子,A容易获得,但是B去哪找呢?人家没得肺癌也不会来医院呀?),因此基于医院的研究获得RR值比较困难。但是这不代表医院的什么研究都得不到RR值,比如设计这样一个队列研究,探讨肺癌术后化疗和不化疗的患者肺癌复发情况,这时我们就可以计算复发率,在比较化疗的效应的时候我们就可以计算化疗效应的RR值。

HR

(Risk Ratio):

风险函数比,是生存分析资料中用于估计因为某种因素的存在而使死亡/缓解/复发等风险改变的倍数。其计算公式为:

HR=暴露组的风险函数 h1(t)/非暴露组的风险函数 h2(t)(t 指在相同的时间点上,风险函数指危险率函数、条件死亡率、瞬时死亡率)   HR主要通过 COX 回归分析得出,需要用软件来算。

HR 与 RR的区别

1)两者均用于前瞻性研究,HR与RR 意思差不多,但HR 还考虑了时间因素,包含了时间效应的 RR 就是 HR;

2)从终点时间的角度来看, RR 考虑了终点事件的差异,HR 不仅考虑了终点事件的有无,还考虑了到达终点所用的时间及截尾数据。

前面都是跟大家回顾下3R的定义和区别,下面回归主题。

如果你主要获得是OR值,OR值主要来源于病例对照研究,在做Meta分析的时候可以把OR转换为RR进行合并。

从计算公式来讲,当研究结局的发生率比较低时,OR≈RR, RR和OR两者之间可以互换

,然后进行合并。那么如果发生率比较大,就需要进行转换,怎么转换呢?主要用下面的公式:OR转换为RR   RR=OR/[(1–P1)+(P1×OR)], 其中P1是对照组结局指标的发生率。

95%CI的计算

SElog(RR)=SElog(OR)×log(RR)/log(OR)


如果你获得的是RR值或者HR值,RR和HR都是前瞻性研究的效应值,前面讲了两个意思实际上差不多,只不过HR引入了时间因素,在做Meta分析时,HR可以直接转换为RR,然后进行合并。

所以从上面来看,选择RR在一般情况下都是通用的,其实最简单的办法就是,如果在纳入的文献中既有OR也有RR,HR,可以把这些都放在一起合并,然后再以3个R为分组因素做亚组分析。

MSE MAE

MSE

定义:MSE(均方误差)函数一般用来检测模型的预测值和真实值之间的偏差。MSE是真实值与预测值的差值的平方然后求和平均。通过平方的形式便于求导,所以常被用作线性回归的损失函数。

MAE

定义:MAE(Mean Absolute Error)平均绝对误差。是绝对误差的平均值,可以更好地反映预测值误差的实际情况。


代码

#install.packages("survival")#引用包library(survival)
trainFile="TCGA.expTime.txt"#train组输入文件setwd("文件路径")                        #设置工作目录rt=read.table(trainFile, header=T, sep="\t", row.names=1,check.names=F)    #读取train组输入文件#rt[,3:ncol(rt)]=log2(rt[,3:ncol(rt)]+1)rt[,"futime"]=rt[,"futime"]/365#COX模型构建multiCox=coxph(Surv(futime, fustat) ~., data=rt)
multiCox=step(multiCox,direction="both")
multiCoxSum=summary(multiCox)
#输出模型的公式outMultiTab=data.frame()
outMultiTab=cbind(
coef=multiCoxSum$coefficients[,"coef"],
HR=multiCoxSum$conf.int[,"exp(coef)"],
HR.95L=multiCoxSum$conf.int[,"lower .95"],
HR.95H=multiCoxSum$conf.int[,"upper .95"],
pvalue=multiCoxSum$coefficients[,"Pr(>|z|)"])
outMultiTab=cbind(id=row.names(outMultiTab),outMultiTab)
outMultiTab=outMultiTab[,1:2]
write.table(outMultiTab,file="multiCox.txt",sep="\t",row.names=F,quote=F)
#输出train组的风险文件(TCGA)trainScore=predict(multiCox,type="risk",newdata=rt)
coxGene=rownames(multiCoxSum$coefficients)
coxGene=gsub("`","",coxGene)
outCol=c("futime","fustat",coxGene)
Risk=as.vector(ifelse(trainScore>median(trainScore),"high","low"))
outTab=cbind(rt[,outCol],riskScore=as.vector(trainScore), Risk)
write.table(cbind(id=rownames(outTab),outTab),file="risk.TCGA.txt",sep="\t",quote=F,row.names=F)












相关文章
|
2月前
【R语言实战】——带有高斯新息的金融时序的GARCH模型拟合预测及VAR/ES风险度量
【R语言实战】——带有高斯新息的金融时序的GARCH模型拟合预测及VAR/ES风险度量
|
2月前
【R语言实战】——带有新息为标准学生t分布的金融时序的GARCH模型拟合预测
【R语言实战】——带有新息为标准学生t分布的金融时序的GARCH模型拟合预测
|
2月前
|
数据可视化 算法
【R语言实战】——kNN和朴素贝叶斯方法实战
【R语言实战】——kNN和朴素贝叶斯方法实战
|
2月前
【R语言实战】——Logistic回归模型
【R语言实战】——Logistic回归模型
|
2月前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
2月前
|
机器学习/深度学习 数据可视化
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为
|
2月前
|
数据采集 数据可视化
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
|
2月前
利用R语言进行典型相关分析实战
利用R语言进行典型相关分析实战
|
2月前
|
Web App开发 数据可视化 数据挖掘
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
|
2月前
|
机器学习/深度学习 算法
R语言分类回归分析考研热现象分析与考研意愿价值变现
R语言分类回归分析考研热现象分析与考研意愿价值变现