全文链接：http://tecdat.cn/?p=2858
本文的目的是对如何在R中进行生存分析进行简短而全面的评估。关于该主题的文献很广泛，仅涉及有限数量的（常见）问题。可用的R包数量反映了对该主题的研究范围。（点击文末“阅读原文”获取完整代码数据）。
视频：R语言生存分析原理与晚期肺癌患者分析案例
拓端
，赞4
R语言生存分析Survival analysis原理与晚期肺癌患者分析案例
R包
可以使用各种R包来解决特定问题。以下是本次用于读取，管理，分析和显示数据的软件包。 运行以下行以安装和加载所需的包。
代码语言：javascript复制if (!require(pacman)) install.packages("pacman")pacman::p_load(tidyverse, survival )**
数据
该评价将基于orca数据集，数据集包含1985年1月1日至2005年12月31日期间芬兰最北部省份诊断为口腔鳞状细胞癌（OSCC）的338名患者的一部分。患者的随访始于癌症诊断之日，并于2008年12月31日死亡，迁移或随访截止日期结束。死亡原因分为两类：（1） ）OSCC死亡; （2）其他原因造成的死亡。
数据集包含以下变量：
id=序号，
sex=性别，类别1 =“女性”，2 =“男性”，
age=诊断癌症日期的年龄（年），
stage=肿瘤的TNM分期（因子）：1 =“I”，...， 4 =“IV”，5 =“unkn”
time=自诊断至死亡或审查的随访时间（以年为单位），
event=结束随访的事件（因子）：1 =正常，2 =口腔癌死亡， 3 =其他原因造成的死亡。
将数据从URL加载到R中。
代码语言：javascript复制head(orca) id sex age stage time event1 1 Male 65.42274 unkn 5.081 Alive2 2 Female 83.08783 III 0.419 Oral ca. death3 3 Male 52.59008 II 7.915 Other death4 4 Male 77.08630 I 2.480 Other death5 5 Male 80.33622 IV 2.500 Oral ca. death6 6 Female 82.58132 IV 0.167 Other deathsummary(orca) id sex age stage time event Min. : 1.00 Female:152 Min. :15.15 I :50 Min. : 0.085 Alive :109 1st Qu.: 85.25 Male :186 1st Qu.:53.24 II :77 1st Qu.: 1.333 Oral ca. death:122 Median :169.50 Median :64.86 III :72 Median : 3.869 Other death :107 Mean :169.50 Mean :63.51 IV :68 Mean : 5.662 3rd Qu.:253.75 3rd Qu.:74.29 unkn:71 3rd Qu.: 8.417 Max. :338.00 Max. :92.24 Max. :23.258 生存数据分析
生存分析侧重于事件数据的时间。在我们的例子中，是诊断后的死亡时间。
为了定义失效时间随机变量，我们需要： 1。时间起源（诊断OSCC）， 2。时间尺度（诊断后的年数，年龄）， 3。事件的定义。我们将首先考虑总死亡率 。
图1：转换的框图。
代码语言：javascript复制 Alive Oral ca. death Other death 109 122 107 FALSE TRUE 109 229 以图形方式显示观察到的随访时间对于生存数据的分析非常有帮助。
OSCC死亡更有可能在诊断后早期发生，而不是其他原因引起的死亡。类型怎么样？
代码语言：javascript复制 'Surv' num [1:338, 1:2] 5.081 0.419 7.915 2.480 2.500 0.167 5.925 1.503 13.333 7.666 ... - attr(*, "dimnames")=List of 2 ..$ : NULL ..$ : chr [1:2] "time" "status" - attr(*, "type")= chr "right"然后将创建的生存对象用作生存分析的其他特定函数中的因变量。
相关视频
拓端
，赞11
拓端
，赞4
估计生存函数
非参数估计
我们将首先介绍一类非参数估计 。
Kaplan–Meier
生存曲线基于每个死亡时间的风险数量和事件数量。包的survfit()创建（估计）生存曲线 。
Call: survfit(formula = Surv(time, all) ~ 1, data = orca) n events *rmean *se(rmean) median 0.95LCL 0.95UCL 338.000 229.000 8.060 0.465 5.418 4.331 6.916 * restricted mean with upper limit = 23.3 函数返回估计的生存曲线的摘要。
代码语言：javascript复制 time n.risk n.event n.censor surv std.err upper lower1 0.085 338 2 0 0.9940828 0.004196498 1.0000000 0.98594012 0.162 336 2 0 0.9881657 0.005952486 0.9997618 0.97670413 0.167 334 4 0 0.9763314 0.008468952 0.9926726 0.96025924 0.170 330 2 0 0.9704142 0.009497400 0.9886472 0.95251755 0.246 328 1 0 0.9674556 0.009976176 0.9865584 0.94872286 0.249 327 1 0 0.9644970 0.010435745 0.9844277 0.9449699ggsurvplot()的survminer提供了估计的生存曲线的信息性说明。
默认的KM图表显示了生存函数。
点击标题查阅往期内容
【视频】分类模型评估:精确率、召回率、ROC曲线、AUC与R语言生存分析时间依赖性ROC实现
左右滑动查看更多
01
02
03
04
生存曲线估算
生存曲线在精算师和人口统计学中非常普遍。它特别适用于分组数据。
为了在实际示例中显示此方法，我们首先需要创建聚合数据，即将后续分组并在每个层中计算风险。
基于分组的数据，我们估计会用生存曲线。
代码语言：javascript复制 nsubs nlost nrisk nevent surv pdf hazard se.surv se.pdf se.hazard0-1 338 0 338.0 64 1.0000 0.1893 0.2092 0.0000 0.0213 0.02601-2 274 4 272.0 41 0.8107 0.1222 0.1630 0.0213 0.0179 0.02542-3 229 9 224.5 21 0.6885 0.0644 0.0981 0.0252 0.0136 0.02143-4 199 12 193.0 20 0.6241 0.0647 0.1093 0.0265 0.0140 0.02444-5 167 9 162.5 13 0.5594 0.0448 0.0833 0.0274 0.0121 0.02315-6 145 14 138.0 13 0.5146 0.0485 0.0989 0.0279 0.0131 0.02746-7 118 5 115.5 8 0.4662 0.0323 0.0717 0.0283 0.0112 0.02547-8 105 8 101.0 9 0.4339 0.0387 0.0933 0.0286 0.0126 0.03118-9 88 7 84.5 1 0.3952 0.0047 0.0119 0.0288 0.0047 0.01199-10 80 4 78.0 8 0.3905 0.0401 0.1081 0.0288 0.0137 0.038210-11 68 4 66.0 5 0.3505 0.0266 0.0787 0.0291 0.0116 0.0352Nelson-Aalen估计
图形比较
可以绘制不同的生存函数估计值来评估潜在的差异。
可以从估计的生存曲线导出诸如分位数的集中趋势的度量。
代码语言：javascript复制 q km.quantile km.lower km.upper fh.quantile fh.lower fh.upper25 0.25 1.333 1.084 1.834 1.333 1.084 1.74750 0.50 5.418 4.331 6.916 5.418 4.244 6.91375 0.75 13.673 11.748 16.580 13.673 11.748 15.833估计半数人的寿命超过5.4年。 第一个四分之一的人在1.3年内死亡，而前四分之三的人的寿命超过1.3岁。 前三分之三的人在13.7年内死亡，而前四分之一的人死亡时间超过13.7岁。
估计量的图形表示（基于使用KM的生存曲线）
参数估算
我们将考虑三种常见的选择：指数，Weibull和log-logistic模型。
代码语言：javascript复制flexsurvreg(formula = su_obj ~ 1, data = orca, dist = "exponential")Estimates: est L95% U95% se rate 0.11967 0.10513 0.13621 0.00791N = 338, Events: 229, Censored: 109Total time at risk: 1913.673Log-likelihood = -715.1802, df = 1AIC = 1432.36同样，可以用非参数估计图形地比较不同的方法
生存曲线的比较
例如，肿瘤阶段是癌症存活研究中的重要预后因素。我们可以估计和绘制不同颜色的不同组（阶段）的生存曲线。
代码语言：javascript复制 stage D Y x pt rate lower upper conf.level1 I 25 336.776 25 336.776 0.07423332 0.04513439 0.1033322 0.952 II 51 556.700 51 556.700 0.09161128 0.06646858 0.1167540 0.953 III 51 464.836 51 464.836 0.10971611 0.07960454 0.1398277 0.954 IV 57 262.552 57 262.552 0.21709985 0.16073995 0.2734597 0.955 unkn 45 292.809 45 292.809 0.15368380 0.10878136 0.1985862 0.95通常，与具有高阶段肿瘤的患者相比，具有较低阶段肿瘤的诊断患者具有较低的（死亡率）。可以使用survfit()函数执行生存函数的整体比较。
Call: survfit(formula = su_obj ~ stage, data = orca) n events median 0.95LCL 0.95UCLstage=I 50 25 10.56 6.17 NAstage=II 77 51 7.92 4.92 13.34stage=III 72 51 7.41 3.92 9.90stage=IV 68 57 2.00 1.08 4.82stage=unkn 71 45 3.67 2.83 8.17由于低肿瘤阶段的发病率较低，因此肿瘤分期增加的中位生存时间也会减少。可以观察到相同的行为，分别针对不同的肿瘤阶段绘制KM生存曲线。
也可以为每个阶段级别构建整个生存表。这里是每个肿瘤阶段生存表的前3行。
代码语言：javascript复制# Groups: strata [5] time n.risk n.event n.censor surv std.err upper lower strata <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> 1 0.17 50 1 0 0.98 0.0202 1 0.942 I 2 0.498 49 1 0 0.96 0.0289 1 0.907 I 3 0.665 48 1 0 0.94 0.0357 1 0.876 I 4 0.419 77 1 0 0.987 0.0131 1 0.962 II 5 0.498 76 1 0 0.974 0.0186 1 0.939 II 6 0.665 75 1 0 0.961 0.0229 1 0.919 II 7 0.167 72 1 0 0.986 0.0140 1 0.959 III 8 0.249 71 1 0 0.972 0.0199 1 0.935 III 9 0.413 70 1 0 0.958 0.0246 1 0.913 III 10 0.085 68 2 0 0.971 0.0211 1 0.931 IV 11 0.162 66 1 0 0.956 0.0261 1 0.908 IV 12 0.167 65 1 0 0.941 0.0303 0.999 0.887 IV 13 0.162 71 1 0 0.986 0.0142 1 0.959 unkn 14 0.167 70 2 0 0.958 0.0249 1 0.912 unkn 15 0.17 68 1 0 0.944 0.0290 0.999 0.892 unkn arrange_ggsurvplots(glist, print = TRUE, ncol = 2, nrow = 1)
Mantel-Haenszel logrank测试
默认参数rho = 0实现log-rank或Mantel-Haenszel测试。
Call:survdiff(formula = su_obj ~ stage, data = orca) N Observed Expected (O-E)^2/E (O-E)^2/Vstage=I 50 25 39.9 5.573 6.813stage=II 77 51 63.9 2.606 3.662stage=III 72 51 54.1 0.174 0.231stage=IV 68 57 33.2 16.966 20.103stage=unkn 71 45 37.9 1.346 1.642 Chisq= 27.2 on 4 degrees of freedom, p= 2e-05 Peto＆Peto Gehan-Wilcoxon测试
代码语言：javascript复制survdiff(formula = su_obj ~ stage, data = orca, rho = 1) N Observed Expected (O-E)^2/E (O-E)^2/Vstage=I 50 14.5 25.2 4.500 7.653stage=II 77 29.3 39.3 2.549 4.954stage=III 72 30.7 33.8 0.284 0.521stage=IV 68 40.3 22.7 13.738 21.887stage=unkn 71 32.0 25.9 1.438 2.359 Chisq= 30.9 on 4 degrees of freedom, p= 3e-06 不同的测试使用不同的权重来比较生存函数。在实际例子中，他们给出了可比较的结果，表明不同肿瘤阶段的生存函数是不同的。
建模生存数据
当比较因子水平的生存函数时，非参数检验特别可行。它们非常强大，高效，通常简单/直观。 然而，随着感兴趣因素的数量增加，非参数测试变得难以进行和解释。相反，回归模型对于探索生存与预测因子之间的关系更为灵活。
我们将介绍两种不同的广泛模型：半参数（即比例风险）和参数模型。
CoxPH模型
在我们的例子中，我们将考虑将死亡时间建模为性别，年龄和肿瘤阶段的函数。
可以使用coxph()功能来建立Cox比例风险模型survival。
summary(m1)Call:coxph(formula = su_obj ~ sex I((age - 65)/10) stage, data = orca) n= 338, number of events= 229 coef exp(coef) se(coef) z Pr(>|z|)sexMale 0.35139 1.42104 0.14139 2.485 0.012947I((age - 65)/10) 0.41603 1.51593 0.05641 7.375 1.65e-13stageII 0.03492 1.03554 0.24667 0.142 0.887421stageIII 0.34545 1.41262 0.24568 1.406 0.159708stageIV 0.88542 2.42399 0.24273 3.648 0.000265stageunkn 0.58441 1.79393 0.25125 2.326 0.020016 exp(coef) exp(-coef) lower .95 upper .95sexMale 1.421 0.7037 1.0771 1.875I((age - 65)/10) 1.516 0.6597 1.3573 1.693stageII 1.036 0.9657 0.6386 1.679stageIII 1.413 0.7079 0.8728 2.286stageIV 2.424 0.4125 1.5063 3.901stageunkn 1.794 0.5574 1.0963 2.935Concordance= 0.674 (se = 0.02 )Rsquare= 0.226 (max possible= 0.999 )Likelihood ratio test= 86.76 on 6 df, p=<2e-16Wald test = 80.5 on 6 df, p=3e-15Score (logrank) test = 82.86 on 6 df, p=9e-16我们可以检查数据是否与每个变量的比例风险假设分别和全局一致。
代码语言：javascript复制 rho chisq psexMale -0.00137 0.000439 0.983I((age - 65)/10) 0.07539 1.393597 0.238stageII -0.04208 0.411652 0.521stageIII -0.06915 1.083755 0.298stageIV -0.10044 2.301780 0.129stageunkn -0.09663 2.082042 0.149GLOBAL NA 4.895492 0.557显然没有找到违反比例假设的证据。
Cox模型的结果表明性别，年龄和阶段的显着影响。特别是，每增加10年，死亡率就会增加50％。与男性和女性相比，全因死亡率的HR为1.42。此外，估计数中第一阶段和第二阶段之间未发现任何差异。因此，谨慎的做法是将这些主题从数据中排除，并将前两个阶段组合为一个。
代码语言：javascript复制round(ci.exp(m2), 4) exp(Est.) 2.5% 97.5%sexMale 1.3284 0.9763 1.8074I((age - 65)/10) 1.4624 1.2947 1.6519st3III 1.3620 0.9521 1.9482st3IV 2.3828 1.6789 3.3818显示和图形化比较多变量Cox模型的结果的便捷方式是通过森林图。
让我们逐步绘制预测的生存曲线，根据拟合的模型确定性别和年龄的值
代码语言：javascript复制newd sex age st3 id1 Male 40 I II 12 Female 40 I II 23 Male 80 I II 34 Female 80 I II 45 Male 40 III 56 Female 40 III 67 Male 80 III 78 Female 80 III 89 Male 40 IV 910 Female 40 IV 1011 Male 80 IV 1112 Female 80 IV 12AFT模型
参数模型假设生存时间的分布。
代码语言：javascript复制Call:flexsurvreg(formula = Surv(time, all) ~ sex I((age - 65)/10) st3, data = orca2, dist = "weibull")Estimates: data mean est L95% U95% se exp(est) L95% shape NA 0.93268 0.82957 1.04861 0.05575 NA NAscale NA 13.53151 9.97582 18.35456 2.10472 NA NAsexMale 0.53184 -0.33905 -0.66858 -0.00951 0.16813 0.71245 0.51243I((age - 65)/10) -0.15979 -0.41836 -0.54898 -0.28773 0.06665 0.65813 0.57754st3III 0.26966 -0.32567 -0.70973 0.05839 0.19595 0.72204 0.49178st3IV 0.25468 -0.95656 -1.33281 -0.58030 0.19197 0.38421 0.26374 U95% shape NAscale NAsexMale 0.99053I((age - 65)/10) 0.74996st3III 1.06012st3IV 0.55973N = 267, Events: 184, Censored: 83Total time at risk: 1620.864Log-likelihood = -545.858, df = 6AIC = 1103.716可以证明，假设指数或威布尔分布的AFT模型可以重新参数化为比例风险模型。
显示eha。
Call:weibreg(formula = Surv(time, all) ~ sex I((age - 65)/10) st3, data = orca2)Covariate Mean Coef Exp(Coef) se(Coef) Wald psex Female 0.490 0 1 (reference) Male 0.510 0.316 1.372 0.156 0.043 I((age - 65)/10) -0.522 0.390 1.477 0.062 0.000 st3 I II 0.551 0 1 (reference) III 0.287 0.304 1.355 0.182 0.095 IV 0.162 0.892 2.440 0.178 0.000 log(scale) 2.605 13.532 0.156 0.000 log(shape) -0.070 0.933 0.060 0.244 Events 184 Total time at risk 1620.9 Max. log. likelihood -545.86 LR test statistic 68.7 Degrees of freedom 4 Overall p-value 4.30767e-14系数的（指数）具有与Cox比例模型的系数的等效解释。
通过将参数提供fn给summary或plot方法，可以汇总或绘制拟合模型的参数的任何函数。例如，Weibull模型下的中位存活率可以概括为
newd <- data.frame(sex = c("Male", "Female"), age = 65, st3 = "I II")summary(m2w, newdata = newd, fn = median.weibull, t = 1, B = 10000)sex=Male, I((age - 65)/10)=0, st3=I II time est lcl ucl1 1 6.507834 4.898889 8.631952sex=Female, I((age - 65)/10)=0, st3=I II time est lcl ucl1 1 9.134466 6.801322 12.33771将结果与Cox模型的结果进行比较。
代码语言：javascript复制survfit(m2, newdata = newd)Call: survfit(formula = m2, newdata = newd) n events median 0.95LCL 0.95UCL1 267 184 7.00 5.25 10.62 267 184 9.92 7.33 13.8泊松回归
可以证明，Cox模型在数学上等效于对数据的特定变换的泊松回归模型。
我们首先定义观察事件（all == 1）的唯一时间，并使用包中的survSplit()函数survival来分割数据。
head(orca_splitted, 15)拟合条件泊松回归，其中时间的影响（作为因子变量）可以被边缘化（不估计来提高计算效率）。
mod_poi <- gnm(all ~ sex I((age-65)/10) st3, data = orca_splitted, family = poisson, eliminate = factor(time))summary(mod_poi)将从条件Poisson获得的估计值与cox比例风险模型进行比较。
代码语言：javascript复制round(data.frame(cox = ci.exp(m2), poisson = ci.exp(mod_poi)), 4) cox.exp.Est.. cox.2.5. cox.97.5. poisson.exp.Est.. poisson.2.5. poisson.97.5.sexMale 1.3284 0.9763 1.8074 1.3284 0.9763 1.8074I((age - 65)/10) 1.4624 1.2947 1.6519 1.4624 1.2947 1.6519st3III 1.3620 0.9521 1.9482 1.3620 0.9521 1.9482st3IV 2.3828 1.6789 3.3818 2.3828 1.6789 3.3818如果我们想要估计基线风险，我们还需要估计泊松模型中时间的影响。
代码语言：javascript复制orca_splitted$dur <- with(orca_splitted, time - tstart)mod_poi2 <- glm(all ~ -1 factor(time) sex I((age-65)/10) st3, data = orca_splitted, family = poisson, offset = log(dur))基线风险包括阶梯函数，其中速率在每个时间间隔内是恒定的。
代码语言：javascript复制newd <- data.frame(time = cuts, dur = 1, sex = "Female", age = 65, st3 = "I II")blhaz <- data.frame(ci.pred(mod_poi2, newdata = newd))ggplot(blhaz, aes(x = c(0, cuts[-138]), y = Estimate, xend = cuts, yend = Estimate)) geom_segment() scale_y_continuous(trans = "log", limits = c(.05, 5), breaks = pretty_breaks()) theme_classic() labs(x = "Time (years)", y = "Baseline hazard")
更好的方法是通过使用例如具有节点（k ）的样条来灵活地模拟基线风险。
代码语言：javascript复制 exp(Est.) 2.5% 97.5%(Intercept) 0.074 0.040 0.135ns(time, knots = k)1 0.402 0.177 0.912ns(time, knots = k)2 1.280 0.477 3.432ns(time, knots = k)3 0.576 0.220 1.509ns(time, knots = k)4 1.038 0.321 3.358ns(time, knots = k)5 4.076 0.854 19.452ns(time, knots = k)6 1.040 0.171 6.314sexMale 1.325 0.975 1.801I((age - 65)/10) 1.469 1.300 1.659st3III 1.360 0.952 1.942st3IV 2.361 1.665 3.347
比较不同的策略
我们可以根据特定协变量模式的预测生存曲线比较之前的策略，如65岁的女性患有肿瘤I期或II期。
代码语言：javascript复制newd <- data.frame(sex = "Female", age = 65, st3 = "I II")生存函数的图形表示便于比较。
其他分析
非线性
我们假设年龄对（log）死亡率的影响是线性的。放宽这一假设的可能策略是拟合Cox模型，其中年龄用二次效应建模。
代码语言：javascript复制Call:coxph(formula = Surv(time, all) ~ sex I(age - 65) I((age - 65)^2) st3, data = orca2) n= 267, number of events= 184 coef exp(coef) se(coef) z Pr(>|z|)sexMale 2.903e-01 1.337e 00 1.591e-01 1.825 0.0681I(age - 65) 3.868e-02 1.039e 00 6.554e-03 5.902 3.59e-09I((age - 65)^2) 9.443e-05 1.000e 00 3.576e-04 0.264 0.7917st3III 3.168e-01 1.373e 00 1.838e-01 1.724 0.0847st3IV 8.691e-01 2.385e 00 1.787e-01 4.863 1.16e-06 exp(coef) exp(-coef) lower .95 upper .95sexMale 1.337 0.7481 0.9787 1.826I(age - 65) 1.039 0.9621 1.0262 1.053I((age - 65)^2) 1.000 0.9999 0.9994 1.001st3III 1.373 0.7284 0.9576 1.968st3IV 2.385 0.4193 1.6801 3.385Concordance= 0.674 (se = 0.022 )Rsquare= 0.216 (max possible= 0.999 )Likelihood ratio test= 64.89 on 5 df, p=1e-12Wald test = 63.11 on 5 df, p=3e-12Score (logrank) test = 67.64 on 5 df, p=3e-13非线性（即二次项）的值很高，因此没有证据可以拒绝零假设（即线性假设是合适的）。
如果关系是非线性的，则年龄系数不再可以直接解释。我们可以将HR作为年龄的函数以图形方式呈现。我们需要指定一个指示值; 我们选择65岁的中位年龄值。
代码语言：javascript复制age <- seq(20, 80, 1) - 65 geom_vline(xintercept = 65, lty = 2) geom_hline(yintercept = 1, lty = 2)时间依赖系数
该cox.zph()函数可用于绘制个体预测因子随时间的影响，因此可用于诊断和理解非比例风险。
我们可以通过拟合的阶梯函数来放宽比例风险假设，这意味着在不同的时间间隔内有不同的
包中的survSplit()函数survival将数据集划分。
id sex age stage event st3 tstart time all tgroup1 2 Female 83.08783 III Oral ca. death III 0 0.419 1 12 3 Male 52.59008 II Other death I II 0 5.000 0 13 3 Male 52.59008 II Other death I II 5 7.915 1 24 4 Male 77.08630 I Other death I II 0 2.480 1 15 5 Male 80.33622 IV Oral ca. death IV 0 2.500 1 16 6 Female 82.58132 IV Other death IV 0 0.167 1 1 I((age - 65)/10) st3, data = orca3) coef exp(coef) se(coef) z pI((age - 65)/10) 0.38184 1.46498 0.06255 6.104 1.03e-09st3III 0.28857 1.33451 0.18393 1.569 0.1167st3IV 0.87579 2.40076 0.17963 4.876 1.09e-06relevel(sex, 2)Male:strata(tgroup)tgroup=1 0.42076 1.52312 0.19052 2.209 0.0272relevel(sex, 2)Female:strata(tgroup)tgroup=1 NA NA 0.00000 NA NArelevel(sex, 2)Male:strata(tgroup)tgroup=2 -0.10270 0.90240 0.28120 -0.365 0.7149relevel(sex, 2)Female:strata(tgroup)tgroup=2 NA NA 0.00000 NA NArelevel(sex, 2)Male:strata(tgroup)tgroup=3 1.13186 3.10142 1.09435 1.034 0.3010relevel(sex, 2)Female:strata(tgroup)tgroup=3 NA NA 0.00000 NA NALikelihood ratio test=68.06 on 6 df, p=1.023e-12n= 416, number of events= 184 虽然不显着，但男女比较的风险比在第二时期（5至15年）低于1，而在其他两个时期高于1。
模拟生存百分位数
一个不同但有趣的方法包括模拟生存时间的百分位数。
代码语言：javascript复制Call:ctqr(formula = Surv(time, all) ~ st3, data = orca2, p = p)Coefficients: p = 0.25(Intercept) 2.665 st3III -1.369 st3IV -1.877 Degrees of freedom: 267 total; 225 residualsβ0= 2.665 是参考组中死亡概率等于0.25的时间。另一个被解释为相对度量。
该信息可以直观地比较在肿瘤阶段的水平上分别估计的生存曲线。
代码语言：javascript复制 p = c(p, p - .005, p .005))[-1, ] = 1 - p, xend = time_ref, yend = 1 - p))
对Cox模型中评估生存时间百分位数的可能差异，作为诊断性别和肿瘤阶段年龄的函数。
代码语言：javascript复制ctqr(formula = Surv(time, all) ~ sex I((age - 65)/10) st3, data = orca2, p = seq(0.1, 0.7, 0.1))Coefficients: p = 0.1 p = 0.2 p = 0.3 p = 0.4 p = 0.5 p = 0.6 p = 0.7 (Intercept) 1.44467 2.44379 4.65302 7.73909 10.81386 12.18348 15.19359sexMale -0.09218 -0.27385 -0.85720 -2.49580 -3.27962 -2.81428 -4.01656I((age - 65)/10) -0.19026 -0.39819 -1.20278 -1.93144 -2.39229 -3.03915 -3.52711st3III -0.60994 -1.08534 -1.89357 -2.23741 -3.10478 -2.00037 -1.59213st3IV -1.07679 -1.59566 -2.92700 -3.16652 -4.74759 -4.80838 -5.25810Degrees of freedom: 267 total; 220 residuals结果包括不同百分位数下每种协变量的生存时间差异。
代码语言：javascript复制coef_q <- data.frame(coef(fit_q)) %>% .96 * se )
或者，可以针对一组特定的协方差模式预测生存时间的百分位数。
CIF累积发生率函数
在竞争风险情景中，Kaplan-Meier对特定原因生存的估计通常是不合适的。 我们将考虑事件的累积发生率函数（CIF）
CIF
mstate计算竞争事件的非参数CIF（也称为Aalen-Johansen估计）和相关的标准误差。
head(cif) time Surv CI.Oral ca. death CI.Other death seSurv seCI.Oral ca. death1 0.085 0.9925094 0.007490637 0.000000000 0.005276805 0.0052768052 0.162 0.9887640 0.011235955 0.000000000 0.006450534 0.0064505343 0.167 0.9812734 0.011235955 0.007490637 0.008296000 0.0064505344 0.170 0.9775281 0.011235955 0.011235955 0.009070453 0.0064505345 0.249 0.9737828 0.011235955 0.014981273 0.009778423 0.0064505346 0.252 0.9662921 0.014981273 0.018726592 0.011044962 0.007434315 seCI.Other death1 0.0000000002 0.0000000003 0.0052768054 0.0064505345 0.0074343156 0.008296000我们可以绘制CIF以及生存函数。
通过因子变量的水平来估计累积发生率函数 。
代码语言：javascript复制grid.arrange( ncol = 2)我们可以看到，IV期口腔癌死亡的CIF高于III，甚至更高于I II。相反，对于其他原因死亡率，曲线似乎不随肿瘤阶段而变化。
当我们想要在竞争风险设置中对生存数据进行建模时，有两种常见的策略可以解决不同的问题：
- 针对事件特定风险的Cox模型，例如，兴趣在于预测因素对死亡率的生物效应非常疾病。
- 当我们想要评估因子对事件总体累积发生率的影响时。
CIF Cox模型
代码语言：javascript复制 round(ci.exp(m2haz2), 4) exp(Est.) 2.5% 97.5%sexMale 1.8103 1.1528 2.8431I((age - 65)/10) 1.4876 1.2491 1.7715st3III 1.2300 0.7488 2.0206st3IV 1.6407 0.9522 2.8270原因特异性Cox模型的结果与原因特异性CIF的图形表示一致，即肿瘤IV期仅是口腔癌死亡率的重要风险因素。年龄增加与两种原因的死亡率增加相关（口腔癌死亡率HR = 1.42，其他原因死亡率HR = 1.48）。仅根据其他原因死亡率观察到性别差异（HR = 1.8）。
CRR模型
crr()在cmprsk竞争风险的情况下，包中的函数可用于子分布函数的回归建模。
Call:crr(ftime = time, fstatus = event, cov1 = model.matrix(m2), failcode = "Oral ca. death") coef exp(coef) se(coef) z p-valuesexMale -0.0953 0.909 0.213 -0.447 6.5e-01I((age - 65)/10) 0.2814 1.325 0.093 3.024 2.5e-03st3III 0.3924 1.481 0.258 1.519 1.3e-01st3IV 1.0208 2.775 0.233 4.374 1.2e-05 exp(coef) exp(-coef) 2.5% 97.5%sexMale 0.909 1.100 0.599 1.38I((age - 65)/10) 1.325 0.755 1.104 1.59st3III 1.481 0.675 0.892 2.46st3IV 2.775 0.360 1.757 4.39Num. cases = 267Pseudo Log-likelihood = -501 Pseudo likelihood ratio test = 31.4 on 4 df,m2fg2 <- with(orca2, crr(time, event, cov1 = model.matrix(m2), failcode = "Other death"))summary(m2fg2, Exp = T)Competing Risks RegressionCall:crr(ftime = time, fstatus = event, cov1 = model.matrix(m2), failcode = "Other death") coef exp(coef) se(coef) z p-valuesexMale 0.544 1.723 0.2342 2.324 0.020I((age - 65)/10) 0.197 1.218 0.0807 2.444 0.015st3III 0.130 1.139 0.2502 0.521 0.600st3IV -0.212 0.809 0.2839 -0.748 0.450 exp(coef) exp(-coef) 2.5% 97.5%sexMale 1.723 0.580 1.089 2.73I((age - 65)/10) 1.218 0.821 1.040 1.43st3III 1.139 0.878 0.698 1.86st3IV 0.809 1.237 0.464 1.41Num. cases = 267Pseudo Log-likelihood = -471 Pseudo likelihood ratio test = 9.43 on 4 df,还有问题，联系我们！
