写在前面
PSYCH
上次分享中，我们主要介绍了生存分析的概念，以及各项内在定义，数据收集格式。本期内容主要为如何进行生存分析，以及分析相关的危险因素等。
Vol.1
生存分析的目的与方法
1、描述生存过程：
绘制K-M曲线（kaplan-meier曲线），主要用于描述生存率、死亡速度、中位生存时间。它是一种无参数法（non-parametric）来从观察对象的生存时间，估计生存概率的方法。对于研究中的第n个时间点tn，生存概率可以计算为：
而其中，S（tn-1）是指tn-1的时候生存概率是多少，dn是指在时间点tn所发生的事件数，rn指快要到时间点tn时还存货的人，则t0=0, S(0)=1。
2、比较生存过程：
展示组间存在的差异，进行组间比较。基于生存曲线上每个点的数据，合成一个全面的组间比较结果。需要观察生存曲线的状态，来决定用何种统计方式来进行分组比较。（对于前期差异较大的数据使用Breslow进行比较，对于远期差异较大的数据使用LogRank进行比较。）
LogRank：对每一个点赋予的权重值是相同的，但对于结局事件的远期差别赋予更大的权重。所以对远期差异更敏感。
Breslow：给随访时期的前期的时间点权重更大，所以对近期差异更敏感。
3、计算矫正后的风险
挖掘影响生存结局的潜在因素，并得出矫正后的风险比（HR）。
Vol.2
Cox比例风险模型的构建
1、目的：
用于寻找阳性时间发生速度（风险值）的影响因素（危险因素/保护因素）
2、定义：
当进行生存分析模型构建时（K-M），属于单变量的分析过程，此时忽略了其他变量共同的影响，和变量之间的相互作用。为了解决以上问题，需要构建cox比例风险回归模型。
其中，h0(t)是基准风险方程，可以是任意一个针对时间t的非负方程；Xi是实例i的特征向量；β是参数向量，该向量通过最大化cox部分似然得到。
3、方法：
①筛选有差异的变量，使用单因素cox分析，或者其他筛选变量的方法（如lasso，随机森林等等），再将单因素分析具有统计学意义的变量纳入多因素cox回归中。
②在多因素cox回归中，相当于进行变量之间的矫正，分析独立危险因素或保护因素。
Vol.3
生存分析的R语言实践
首先，根据分析目的，加载相应的R包：
library(ggplot2)
library(survminer)
library(survival)
library(readxl)
其中，survminer与survival为生存分析的主要R包，其余两个包为辅助作图与分析数据。
载入数据：PRDATA<- read_excel("D:/R/PRdata.xlsx")
第一步：进行无分组的生存分析：
#用survival包的Surv函数对时间及生存状态进行拟合（K-M）
fit<- survfit(Surv(time,status)~1,data= PRDATA)
#绘图
ggsurvplot(fit,palette="#C96868")
得到一个带有下降过程的生存曲线，其中纵坐标代表生存概率，横坐标代表由各个样本的存活时间。
第二步：计算分组的生存分布情况
#使用性别进行分组
fit1<-survfit(Surv(time,status)~sex,data= PRDATA)
ggsurvplot(fit1, censor.shape="|", censor.size=4, linetype="strata", conf.int=TRUE,
pval=TRUE, data= PRDATA)
得到数据形式如下图所示，由于基于性别进行分组cox分析时，P值大于0.05，生存曲线相交，代表两组实际上的生存曲线没有统计学差异，因此性别不是影响此疾病生存概率的影响因素。
第三步：计算累计风险曲线
#累计风险曲线
ggsurvplot(fit1, data= PRDATA, conf.int= TRUE, fun="cumhaz")
两组人群，均随着时间的增加，死亡风险在不断上升。
#分组与整体联合绘制
ggsurvplot(fit1, data=PRDATA, conf.int=TRUE,
pval=TRUE,surv.median.line = "hv",add.all= TRUE
)#添加中位生存时间与总生存曲线。
#添加生存表
ggsurvplot(fit1, data=PRDATA, conf.int=TRUE,
pval=TRUE,surv.median.line = "hv",
risk.table= TRUE,
risk.table.col="strata",
legend.labs=c("Male","Female"),
risk.table.height=0.25,
ggtheme = theme_bw())#添加中位生存时间与总生存曲线。
第四步：计算COX模型
##单变量Cox回归
res.cox<- coxph(Surv(time,status)~sex,data= PRDATA,method="efron")
res.cox
此结果与上述分组生存分析计算结果大致相同。
同时可以采用：summary(res.cox)来产生更完整的报告。
如果将单变量coxph函数同时用于多个协变量：
covariates<- c(“age”,”sex”, “...”)
univ_formulas<- sapply(covariates, function(x) as.formula(paste(‘Surv(time,status)~, x)))
univ_models<- lapply( univ_formula, function(x){coxph(x,data = PRDATA0)})
最后需要进行多因素cox回归：
将所有单因素cox回归中有显著差异的变量输入多变量cox回归中：
res.cox<- coxph(Surv(time, status)~ age+sex+..., data= PRDATA)
summary(res.cox)
在summary结果中，可以看到最下方的三行中输出检验结果，在似然比估计，wald检验以及log-rank检验中是否具有统计学意义，并可以求出各回归中变量的系数、95%置信区间以及P值。
在构建模型后，可以对模型进行假设有效性检验：
测试比例危险假设，
test.ph<- cox.zph(res.cox)
test.ph
当输出的协变量P值均大于0.05是，我们可以认为假设成立，不拒绝原假设，原比例风险模型成立。
本期主要讲述了生存分析在R上的实践过程，下一期我们会主要介绍生存分析在SPSS上的实践过程。
PSYCH统计实验室
通知公告
网络分析课程目前开放视频课啦！
单次课200元/讲（学生），250元/讲（非学生）
共有四讲内容：
①横断面网络分析简介与基础
②网络分析与因子分析
③交叉滞后网络分析
④时间序列网络分析
购买后开放视频权限14天，可多次申请。
并赠送所有课程相关资料（无PPT）
如果想申请购买，请联系M18812507626
更多资讯
关注我们
文稿：莲花清瘟
排版：Peruere
责编：Wink
审核：摘星
本文由“Psych统计自习室”课题组原创，欢迎转发至朋友圈。如需转载请联系后台，征得作者同意后方可转载。
