R語言生存分析數據分析可視化案例

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

目標

本文的目的是對如何在R中進行生存分析進行簡短而全面的評估。關於該主題的文獻很普遍，僅涉及有限數量的（常見）問題/特徵。
可用的R包數量反映了對該主題的研究範圍。 

---

R包

可使用各類R包來解決特定問題，而且還有替代功能來解決常見問題。如下是本次審查中用於讀取，管理，分析和顯示數據的軟件包。
運行如下行以安裝和加載所需的包。

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ŤT

爲了定義失效時間隨機變量，咱們須要：
1。時間起源（診斷OSCC），
2。時間尺度（診斷後的年數，年齡），
3。事件的定義。咱們將首先考慮總死亡率（或全因死亡率） 。

圖1：轉換的框圖。

Alive Oral ca. death    Other death            109            122            107

FALSE  TRUE   109   229

以圖形方式顯示觀察到的隨訪時間對於生存數據的分析很是有幫助。 

 

OSCC死亡更有可能在診斷後早期發生，而不是其餘緣由引發的死亡。審查的類型怎麼樣？

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

而後將建立的生存對象用做生存分析的其餘特定函數中的響應變量。

---

估計生存功能

非參數估計

咱們將首先介紹一類非參數估計（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,還有問題，聯繫咱們！