R语言生存分析数据分析可视化案例

时间 2020-07-20

标签语言生存分析数据可视化案例繁體版

原文原文链接

原文连接：http://tecdat.cn/?p=2858

目标

本文的目的是对如何在R中进行生存分析进行简短而全面的评估。关于该主题的文献很普遍，仅涉及有限数量的（常见）问题/特征。
可用的R包数量反映了对该主题的研究范围。并发

R包

可使用各类R包来解决特定问题，而且还有替代功能来解决常见问题。如下是本次审查中用于读取，管理，分析和显示数据的软件包。
运行如下行以安装和加载所需的包。dom

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中。测试

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 death

summary(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

生存数据分析

生存分析侧重于事件数据的时间，一般称为故障时间，。在咱们的例子中，是诊断后的死亡时间。ŤTŤ≥ 0T≥0ŤTflex

为了定义失效时间随机变量，咱们须要：
1。时间起源（诊断OSCC），
2。时间尺度（诊断后的年数，年龄），
3。事件的定义。咱们将首先考虑总死亡率（或全因死亡率）。ui

图1：转换的框图。spa

Alive Oral ca. death    Other death            109            122            107

FALSE  TRUE   109   229

以图形方式显示观察到的随访时间对于生存数据的分析很是有帮助。设计

OSCC死亡更有可能在诊断后早期发生，而不是其余缘由引发的死亡。审查的类型怎么样？3d

'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"

而后将建立的生存对象用做生存分析的其余特定函数中的响应变量。rest

估计生存功能

非参数估计

咱们将首先介绍一类非参数估计（a。）。

Kaplan–Meier

生存曲线基于每一个独特死亡时间的风险数量和事件数量的列表。包的survfit()功能survival使用不一样的方法建立（估计）生存曲线。

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

该print()函数仅返回估计的生存曲线的摘要。

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.9449699

该包中ggsurvplot()的专用功能survminer提供了估计的生存曲线的信息性说明。有关?ggsurvplot不一样可能性（参数）的说明。

默认的KM图表显示了生存函数。有几种替代方案/功能可供使用

可升降的或精算的估算器

生命方法在精算师和人口统计学中很是广泛。它特别适用于分组数据。

为了在实际示例中显示此方法，咱们首先须要建立聚合数据，即将后续分组并在每一个层中计算风险，事件和审查的人数。

基于分组的数据，咱们估计会用存活曲线lifetab()的KMsurv包。

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.0352

Nelson-Aalen估计

图形比较

能够绘制不一样的生存函数估计值来评估潜在的差别。

能够从估计的存活曲线导出诸如分位数的集中趋势的度量。

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模型。

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

一样，能够用非参数估计器图形地比较不一样的方法

生存曲线的比较

例如，肿瘤阶段是癌症存活研究中的重要预后因素。咱们能够估计和绘制不一样颜色的不一样组（阶段）的生存曲线。

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行。

# 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测试

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

咱们可使用函数检查数据是否与每一个变量的比例风险假设分别和全局一致cox.zph()。

rho    chisq     p

sexMale          -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。此外，估计数中第一阶段和第二阶段之间未发现任何差别。另外一方面，阶段未知的群体是来自不一样真实阶段的患者的复杂混合物。所以，谨慎的作法是将这些主题从数据中排除，并将前两个阶段组合为一个。

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模型的结果的便捷方式是经过森林图。

让咱们逐步绘制预测的生存曲线，根据拟合的模型肯定性别和年龄的值

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 12

AFT模型

参数模型假设存活时间的分布。

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模型能够从新参数化为比例风险模型（具备来自指数/威布尔分布族的基线危险）。
这可使用包中的weibreg()功能显示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模型的结果进行比较。

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)

包中的gnm()函数gnm适合分裂数据上的条件泊松，其中时间的影响（做为因子变量）能够被边缘化（不估计以提升计算效率）。

mod_poi <- gnm(all ~ sex + I((age-65)/10) + st3, data = orca_splitted,                family = poisson, eliminate = factor(time))summary(mod_poi)

将从条件Poisson得到的估计值与cox比例风险模型进行比较。

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

若是咱们想要估计基线危险，咱们还须要估计泊松模型中时间的影响（OBS咱们还须要包括时间间隔的（对数）长度做为偏移）。

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))

基线危险包括阶梯函数，其中速率在每一个时间间隔内是恒定的。

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 \）的样条来灵活地模拟基线危险。

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期。

newd <- data.frame(sex = "Female", age = 65, st3 = "I+II")

生存函数的图形表示便于比较。

其余分析

非线性

咱们假设年龄对（log）死亡率的影响是线性的。放宽这一假设的可能策略是适合Cox模型，其中年龄用二次效应建模。

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岁的中位年龄值。

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。该cox.zph()函数可用于检查在分割分析中是否仍然偏离比例假设。

模拟生存百分位数

一个不一样但有趣的视角包括模拟生存时间的百分位数。

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=2.665是参考组中死亡几率等于0.25的时间。另外一个被解释为相对度量。

该信息能够直观地比较在肿瘤阶段的水平上分别估计的存活曲线。

p = c(p, p - .005, p + .005))[-1, ]  = 1 - p, xend = time_ref,                                  yend = 1 - p))

对Cox模型中表示的那个进行补充分析_m2_评估生存时间百分位数的可能差别，做为诊断，性别和肿瘤阶段年龄的函数

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

结果包括不一样百分位数下每种协变量的存活时间差别。图形表示能够促进结果的解释。

coef_q <- data.frame(coef(fit_q)) %>%    .96 * se    )

或者，能够针对一组特定的协方差模式预测存活时间的百分位数。

CIF累积发生率函数

对于因某一特定缘由而死亡的风险或危险，一般是主要的兴趣。因为竞争缘由阻止受试者发展事件，可能没法观察到缘由特定事件。竞争事件不只发生在特定缘由的死亡率上，并且每次事件都会阻止并发事件发生时更多。
在咱们的例子中，咱们感兴趣的是模拟口腔癌死亡风险，而且死于其余缘由将被视为竞争事件。

在竞争风险情景中，事件特定的生存得到审查其余事件（也就是天真的Kaplan-Meier对特定缘由生存的估计）一般是不合适的。
咱们将考虑事件的累积发生率函数（CIF）

CIF

包中的Cuminc()函数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（一个堆叠的另外一个）以及派生的无事件生存函数。

已经实施了扩展以经过因子变量的水平来估计累积发生率函数。

grid.arrange(

ncol = 2)

咱们能够看到，IV期口腔癌死亡的CIF高于III，甚至更高于I + II。相反，对于其余缘由死亡率，曲线彷佛不随肿瘤阶段而变化。

当咱们想要在竞争风险设置中对生存数据进行建模时，有两种常见的策略能够解决不一样的问题：
- 针对事件特定危害的Cox模型，例如，兴趣在于预测因素对死亡率的生物效应很是疾病，每每致使相关的结果。
- 当咱们想要评估因子对事件整体累积发生率的影响时，细分模型的细分模型为。Cc

CIF Cox模型

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,还有问题，联系咱们！