目標
本文的目的是對如何在R中進行生存分析進行簡短而全面的評估。關於該主題的文獻很廣泛，僅涉及有限數量的（常見）問題/特征。
可用的R包數量反映了對該主題的研究范圍。
R包
可以使用各種R包來解決特定問題，並且還有替代功能來解決常見問題。以下是本次審查中用於讀取，管理，分析和顯示數據的軟件包。
運行以下行以安裝和加載所需的包。
if (!require(pacman)) install.packages("pacman") pacman::p_load(tidyverse, survival, ggfortify, survminer, plotly, gridExtra, Epi, KMsurv, gnm, cmprsk, mstate, flexsurv, splines, epitools, eha, shiny, ctqr, scales)
數據
該評價將基於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 event
1 1 Male 65.42274 unkn 5.081 Alive
2 2 Female 83.08783 III 0.419 Oral ca. death
3 3 Male 52.59008 II 7.915 Other death
4 4 Male 77.08630 I 2.480 Other death
5 5 Male 80.33622 IV 2.500 Oral ca. death
6 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 生存數據分析
生存分析側重於事件數據的時間，通常稱為故障時間，。在我們的例子中，是診斷后的死亡時間。ŤTŤ≥ 0T≥0ŤT
為了定義失效時間隨機變量，我們需要：
1。時間起源（診斷OSCC），
2。時間尺度（診斷后的年數，年齡），
3。事件的定義。我們將首先考慮總死亡率（或全因死亡率） 。
Alive Oral ca. death Other death
109 122 107 orca <- mutate(orca, all = event != "Alive") table(orca$all)
FALSE TRUE
109 229 以圖形方式顯示觀察到的隨訪時間對於生存數據的分析非常有幫助。
OSCC死亡更有可能在診斷后早期發生，而不是其他原因引起的死亡。審查的類型怎么樣？
str(su_obj) '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 lower
1 0.085 338 2 0 0.9940828 0.004196498 1.0000000 0.9859401
2 0.162 336 2 0 0.9881657 0.005952486 0.9997618 0.9767041
3 0.167 334 4 0 0.9763314 0.008468952 0.9926726 0.9602592
4 0.170 330 2 0 0.9704142 0.009497400 0.9886472 0.9525175
5 0.246 328 1 0 0.9674556 0.009976176 0.9865584 0.9487228
6 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.hazard
0-1 338 0 338.0 64 1.0000 0.1893 0.2092 0.0000 0.0213 0.0260
1-2 274 4 272.0 41 0.8107 0.1222 0.1630 0.0213 0.0179 0.0254
2-3 229 9 224.5 21 0.6885 0.0644 0.0981 0.0252 0.0136 0.0214
3-4 199 12 193.0 20 0.6241 0.0647 0.1093 0.0265 0.0140 0.0244
4-5 167 9 162.5 13 0.5594 0.0448 0.0833 0.0274 0.0121 0.0231
5-6 145 14 138.0 13 0.5146 0.0485 0.0989 0.0279 0.0131 0.0274
6-7 118 5 115.5 8 0.4662 0.0323 0.0717 0.0283 0.0112 0.0254
7-8 105 8 101.0 9 0.4339 0.0387 0.0933 0.0286 0.0126 0.0311
8-9 88 7 84.5 1 0.3952 0.0047 0.0119 0.0288 0.0047 0.0119
9-10 80 4 78.0 8 0.3905 0.0401 0.1081 0.0288 0.0137 0.0382
10-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.upper
25 0.25 1.333 1.084 1.834 1.333 1.084 1.747
50 0.50 5.418 4.331 6.916 5.418 4.244 6.913
75 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.00791
N = 338, Events: 229, Censored: 109
Total time at risk: 1913.673
Log-likelihood = -715.1802, df = 1
AIC = 1432.36同樣，可以用非參數估計器圖形地比較不同的方法
生存曲線的比較
例如，腫瘤階段是癌症存活研究中的重要預后因素。我們可以估計和繪制不同顏色的不同組（階段）的生存曲線。
stage D Y x pt rate lower upper conf.level
1 I 25 336.776 25 336.776 0.07423332 0.04513439 0.1033322 0.95
2 II 51 556.700 51 556.700 0.09161128 0.06646858 0.1167540 0.95
3 III 51 464.836 51 464.836 0.10971611 0.07960454 0.1398277 0.95
4 IV 57 262.552 57 262.552 0.21709985 0.16073995 0.2734597 0.95
5 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.95UCL
stage=I 50 25 10.56 6.17 NA
stage=II 77 51 7.92 4.92 13.34
stage=III 72 51 7.41 3.92 9.90
stage=IV 68 57 2.00 1.08 4.82
stage=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/V
stage=I 50 25 39.9 5.573 6.813
stage=II 77 51 63.9 2.606 3.662
stage=III 72 51 54.1 0.174 0.231
stage=IV 68 57 33.2 16.966 20.103
stage=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/V
stage=I 50 14.5 25.2 4.500 7.653
stage=II 77 29.3 39.3 2.549 4.954
stage=III 72 30.7 33.8 0.284 0.521
stage=IV 68 40.3 22.7 13.738 21.887
stage=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.012947
I((age - 65)/10) 0.41603 1.51593 0.05641 7.375 1.65e-13
stageII 0.03492 1.03554 0.24667 0.142 0.887421
stageIII 0.34545 1.41262 0.24568 1.406 0.159708
stageIV 0.88542 2.42399 0.24273 3.648 0.000265
stageunkn 0.58441 1.79393 0.25125 2.326 0.020016
exp(coef) exp(-coef) lower .95 upper .95
sexMale 1.421 0.7037 1.0771 1.875
I((age - 65)/10) 1.516 0.6597 1.3573 1.693
stageII 1.036 0.9657 0.6386 1.679
stageIII 1.413 0.7079 0.8728 2.286
stageIV 2.424 0.4125 1.5063 3.901
stageunkn 1.794 0.5574 1.0963 2.935
Concordance= 0.674 (se = 0.02 )
Rsquare= 0.226 (max possible= 0.999 )
Likelihood ratio test= 86.76 on 6 df, p=<2e-16
Wald test = 80.5 on 6 df, p=3e-15
Score (logrank) test = 82.86 on 6 df, p=9e-16我們可以使用函數檢查數據是否與每個變量的比例風險假設分別和全局一致cox.zph()。
rho chisq p
sexMale -0.00137 0.000439 0.983
I((age - 65)/10) 0.07539 1.393597 0.238
stageII -0.04208 0.411652 0.521
stageIII -0.06915 1.083755 0.298
stageIV -0.10044 2.301780 0.129
stageunkn -0.09663 2.082042 0.149
GLOBAL 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.8074
I((age - 65)/10) 1.4624 1.2947 1.6519
st3III 1.3620 0.9521 1.9482
st3IV 2.3828 1.6789 3.3818顯示和圖形化比較多變量Cox模型的結果的便捷方式是通過森林圖。
讓我們逐步繪制預測的生存曲線，根據擬合的模型確定性別和年齡的值
newd sex age st3 id
1 Male 40 I+II 1
2 Female 40 I+II 2
3 Male 80 I+II 3
4 Female 80 I+II 4
5 Male 40 III 5
6 Female 40 III 6
7 Male 80 III 7
8 Female 80 III 8
9 Male 40 IV 9
10 Female 40 IV 10
11 Male 80 IV 11
12 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 NA
scale NA 13.53151 9.97582 18.35456 2.10472 NA NA
sexMale 0.53184 -0.33905 -0.66858 -0.00951 0.16813 0.71245 0.51243
I((age - 65)/10) -0.15979 -0.41836 -0.54898 -0.28773 0.06665 0.65813 0.57754
st3III 0.26966 -0.32567 -0.70973 0.05839 0.19595 0.72204 0.49178
st3IV 0.25468 -0.95656 -1.33281 -0.58030 0.19197 0.38421 0.26374
U95%
shape NA
scale NA
sexMale 0.99053
I((age - 65)/10) 0.74996
st3III 1.06012
st3IV 0.55973
N = 267, Events: 184, Censored: 83
Total time at risk: 1620.864
Log-likelihood = -545.858, df = 6
AIC = 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 p
sex
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 ucl
1 1 6.507834 4.898889 8.631952
sex=Female, I((age - 65)/10)=0, st3=I+II
time est lcl ucl
1 1 9.134466 6.801322 12.33771將結果與Cox模型的結果進行比較。
survfit(m2, newdata = newd)Call: survfit(formula = m2, newdata = newd)
n events median 0.95LCL 0.95UCL
1 267 184 7.00 5.25 10.6
2 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.8074
I((age - 65)/10) 1.4624 1.2947 1.6519 1.4624 1.2947 1.6519
st3III 1.3620 0.9521 1.9482 1.3620 0.9521 1.9482
st3IV 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.135
ns(time, knots = k)1 0.402 0.177 0.912
ns(time, knots = k)2 1.280 0.477 3.432
ns(time, knots = k)3 0.576 0.220 1.509
ns(time, knots = k)4 1.038 0.321 3.358
ns(time, knots = k)5 4.076 0.854 19.452
ns(time, knots = k)6 1.040 0.171 6.314
sexMale 1.325 0.975 1.801
I((age - 65)/10) 1.469 1.300 1.659
st3III 1.360 0.952 1.942
st3IV 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.0681
I(age - 65) 3.868e-02 1.039e+00 6.554e-03 5.902 3.59e-09
I((age - 65)^2) 9.443e-05 1.000e+00 3.576e-04 0.264 0.7917
st3III 3.168e-01 1.373e+00 1.838e-01 1.724 0.0847
st3IV 8.691e-01 2.385e+00 1.787e-01 4.863 1.16e-06
exp(coef) exp(-coef) lower .95 upper .95
sexMale 1.337 0.7481 0.9787 1.826
I(age - 65) 1.039 0.9621 1.0262 1.053
I((age - 65)^2) 1.000 0.9999 0.9994 1.001
st3III 1.373 0.7284 0.9576 1.968
st3IV 2.385 0.4193 1.6801 3.385
Concordance= 0.674 (se = 0.022 )
Rsquare= 0.216 (max possible= 0.999 )
Likelihood ratio test= 64.89 on 5 df, p=1e-12
Wald test = 63.11 on 5 df, p=3e-12
Score (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 tgroup
1 2 Female 83.08783 III Oral ca. death III 0 0.419 1 1
2 3 Male 52.59008 II Other death I+II 0 5.000 0 1
3 3 Male 52.59008 II Other death I+II 5 7.915 1 2
4 4 Male 77.08630 I Other death I+II 0 2.480 1 1
5 5 Male 80.33622 IV Oral ca. death IV 0 2.500 1 1
6 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 p
I((age - 65)/10) 0.38184 1.46498 0.06255 6.104 1.03e-09
st3III 0.28857 1.33451 0.18393 1.569 0.1167
st3IV 0.87579 2.40076 0.17963 4.876 1.09e-06
relevel(sex, 2)Male:strata(tgroup)tgroup=1 0.42076 1.52312 0.19052 2.209 0.0272
relevel(sex, 2)Female:strata(tgroup)tgroup=1 NA NA 0.00000 NA NA
relevel(sex, 2)Male:strata(tgroup)tgroup=2 -0.10270 0.90240 0.28120 -0.365 0.7149
relevel(sex, 2)Female:strata(tgroup)tgroup=2 NA NA 0.00000 NA NA
relevel(sex, 2)Male:strata(tgroup)tgroup=3 1.13186 3.10142 1.09435 1.034 0.3010
relevel(sex, 2)Female:strata(tgroup)tgroup=3 NA NA 0.00000 NA NA
Likelihood ratio test=68.06 on 6 df, p=1.023e-12
n= 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.19359
sexMale -0.09218 -0.27385 -0.85720 -2.49580 -3.27962 -2.81428 -4.01656
I((age - 65)/10) -0.19026 -0.39819 -1.20278 -1.93144 -2.39229 -3.03915 -3.52711
st3III -0.60994 -1.08534 -1.89357 -2.23741 -3.10478 -2.00037 -1.59213
st3IV -1.07679 -1.59566 -2.92700 -3.16652 -4.74759 -4.80838 -5.25810
Degrees 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. death
1 0.085 0.9925094 0.007490637 0.000000000 0.005276805 0.005276805
2 0.162 0.9887640 0.011235955 0.000000000 0.006450534 0.006450534
3 0.167 0.9812734 0.011235955 0.007490637 0.008296000 0.006450534
4 0.170 0.9775281 0.011235955 0.011235955 0.009070453 0.006450534
5 0.249 0.9737828 0.011235955 0.014981273 0.009778423 0.006450534
6 0.252 0.9662921 0.014981273 0.018726592 0.011044962 0.007434315
seCI.Other death
1 0.000000000
2 0.000000000
3 0.005276805
4 0.006450534
5 0.007434315
6 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.8431
I((age - 65)/10) 1.4876 1.2491 1.7715
st3III 1.2300 0.7488 2.0206
st3IV 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-value
sexMale -0.0953 0.909 0.213 -0.447 6.5e-01
I((age - 65)/10) 0.2814 1.325 0.093 3.024 2.5e-03
st3III 0.3924 1.481 0.258 1.519 1.3e-01
st3IV 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.38
I((age - 65)/10) 1.325 0.755 1.104 1.59
st3III 1.481 0.675 0.892 2.46
st3IV 2.775 0.360 1.757 4.39
Num. cases = 267
Pseudo 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 Regression
Call:
crr(ftime = time, fstatus = event, cov1 = model.matrix(m2), failcode = "Other death")
coef exp(coef) se(coef) z p-value
sexMale 0.544 1.723 0.2342 2.324 0.020
I((age - 65)/10) 0.197 1.218 0.0807 2.444 0.015
st3III 0.130 1.139 0.2502 0.521 0.600
st3IV -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.73
I((age - 65)/10) 1.218 0.821 1.040 1.43
st3III 1.139 0.878 0.698 1.86
st3IV 0.809 1.237 0.464 1.41
Num. cases = 267
Pseudo Log-likelihood = -471
Pseudo likelihood ratio test = 9.43 on 4 df,