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

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

R语言生存分析数据分析可视化案例(上):https://developer.aliyun.com/article/1493709


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

我们可以检查数据是否与每个变量的比例风险假设分别和全局一致。

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模型可以重新参数化为比例风险模型。

显示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比例模型的系数的等效解释。

通过将参数提供fnsummaryplot方法,可以汇总或绘制拟合模型的参数的任何函数。例如,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)

拟合条件泊松回归,其中时间的影响(作为因子变量)可以被边缘化(不估计来提高计算效率)。

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

如果我们想要估计基线风险,我们还需要估计泊松模型中时间的影响。

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。


模拟生存百分位数


一个不同但有趣的方法包括模拟生存时间的百分位数。

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的时间。另一个被解释为相对度量。

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

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

对Cox模型中评估生存时间百分位数的可能差异,作为诊断性别和肿瘤阶段年龄的函数。

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

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模型,例如,兴趣在于预测因素对死亡率的生物效应非常疾病。
  • 当我们想要评估因子对事件总体累积发生率的影响时。

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,


相关文章
|
5天前
|
JSON 数据挖掘 API
在会议系统工程中,Python可以用于多种任务,如网络请求(用于视频会议的连接和会议数据的传输)、数据分析(用于分析会议参与者的行为或会议效果)等。
在会议系统工程中,Python可以用于多种任务,如网络请求(用于视频会议的连接和会议数据的传输)、数据分析(用于分析会议参与者的行为或会议效果)等。
|
15天前
|
JavaScript Java 测试技术
基于springboot+vue.js+uniapp小程序的数据分析岗位招聘信息与分析附带文章源码部署视频讲解等
基于springboot+vue.js+uniapp小程序的数据分析岗位招聘信息与分析附带文章源码部署视频讲解等
7 0
|
1月前
|
JSON 数据挖掘 API
数据分析实战丨基于pygal与requests分析GitHub最受欢迎的Python库
数据分析实战丨基于pygal与requests分析GitHub最受欢迎的Python库
29 2
|
25天前
|
数据采集 Web App开发 数据可视化
程序员必知:对厦门二手房的数据分析与可视化分析
程序员必知:对厦门二手房的数据分析与可视化分析
26 0
|
2月前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
2月前
|
机器学习/深度学习 数据可视化
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为
|
2月前
|
数据采集 数据可视化
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
|
2月前
利用R语言进行典型相关分析实战
利用R语言进行典型相关分析实战
|
2月前
|
Web App开发 数据可视化 数据挖掘
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)

热门文章

最新文章