library(survival)
library(survminer)
## Warning: package 'ggplot2' was built under R version 4.5.324 生存分析
生存分析是临床常用统计方法，一旦和时间扯上关系，分析就变得复杂多了，此时不再是单一的因变量，还需要考虑时间给因变量和自变量带来的各种影响。
本次主要演示R语言做生存分析的一些方法。比如寿命表、K-M曲线、logrank检验。
本推文不涉及理论，只有实操，想要了解生存分析的理论的请自行学习。不涉及理论，并不代表理论不重要，在以后的机器学习和临床预测模型的相关推文中，会经常用到这些理论，建议大家学习一下。
24.1 生存过程的描述
使用survival包中的lung数据集用于演示，这是一份关于肺癌患者的生存数据。time是生存时间，以天为单位，status是生存状态，1代表删失，2代表死亡。但是一般在生存分析中我们喜欢用1代表死亡，用0代表删失，所以我们更改一下（其实不改也可以，你记住就行）。
df <- lung
df$status <- ifelse(df$status == 2,1,0)
str(df)
## 'data.frame': 228 obs. of 10 variables:
## $ inst : num 3 3 3 5 1 12 7 11 1 7 ...
## $ time : num 306 455 1010 210 883 ...
## $ status : num 1 1 0 1 1 0 1 1 1 1 ...
## $ age : num 74 68 56 57 60 74 68 71 53 61 ...
## $ sex : num 1 1 1 1 1 1 2 2 1 1 ...
## $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
## $ meal.cal : num 1175 1225 NA 1150 NA ...
## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...首先把生存时间和生存状态用Surv()放到一起，可以看到有+的就是截尾数据。
Surv(time = lung$time, event = lung$status)
## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166 170 654
## [13] 728 71 567 144 613 707 61 88 301 81 624 371
## [25] 394 520 574 118 390 12 473 26 533 107 53 122
## [37] 814 965+ 93 731 460 153 433 145 583 95 303 519
## [49] 643 765 735 189 53 246 689 65 5 132 687 345
## [61] 444 223 175 60 163 65 208 821+ 428 230 840+ 305
## [73] 11 132 226 426 705 363 11 176 791 95 196+ 167
## [85] 806+ 284 641 147 740+ 163 655 239 88 245 588+ 30
## [97] 179 310 477 166 559+ 450 364 107 177 156 529+ 11
## [109] 429 351 15 181 283 201 524 13 212 524 288 363
## [121] 442 199 550 54 558 207 92 60 551+ 543+ 293 202
## [133] 353 511+ 267 511+ 371 387 457 337 201 404+ 222 62
## [145] 458+ 356+ 353 163 31 340 229 444+ 315+ 182 156 329
## [157] 364+ 291 179 376+ 384+ 268 292+ 142 413+ 266+ 194 320
## [169] 181 285 301+ 348 197 382+ 303+ 296+ 180 186 145 269+
## [181] 300+ 284+ 350 272+ 292+ 332+ 285 259+ 110 286 270 81
## [193] 131 225+ 269 225+ 243+ 279+ 276+ 135 79 59 240+ 202+
## [205] 235+ 105 224+ 239 237+ 173+ 252+ 221+ 185+ 92+ 13 222+
## [217] 192+ 183 211+ 175+ 197+ 203+ 116 188+ 191+ 105+ 174+ 177+如果只是想要描述一下这份生存数据，可以使用寿命表法或者K-M曲线，在R中可以通过survfit()实现。
# 构建生存曲线
fit <- survfit(Surv(time, status) ~ 1, data = df)
# 寿命表，surv_summary比默认的summary()更好
ss <- surv_summary(fit)
head(ss)
## time n.risk n.event n.censor surv std.err upper lower
## 1 5 228 1 0 0.9956140 0.004395615 1.0000000 0.9870734
## 2 11 227 3 0 0.9824561 0.008849904 0.9996460 0.9655619
## 3 12 224 1 0 0.9780702 0.009916654 0.9972662 0.9592437
## 4 13 223 2 0 0.9692982 0.011786516 0.9919508 0.9471630
## 5 15 221 1 0 0.9649123 0.012628921 0.9890941 0.9413217
## 6 26 220 1 0 0.9605263 0.013425540 0.9861367 0.9355811画出生存曲线，横坐标是生存时间，纵坐标是生存率。
ggsurvplot(fit,
conf.int = TRUE, # 可信区间
palette= 'blue', # 更改配色
surv.median.line = "hv", # 中位生存时间
ggtheme = theme_bw() # 更改主题
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
## Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
24.2 生存过程的比较
如果通过某个变量把数据分为多组，然后检验不同组别之间的生存时间（生存曲线）有无差别，则可以通过logrank检验或者breslow检验。
在R语言中通过survdiff()实现logrank检验。
fit <- survdiff(Surv(time, status) ~ sex, data = df)
fit
## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = df)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 112 91.6 4.55 10.3
## sex=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001可以用神包broom提取数据：
broom::tidy(fit)
## # A tibble: 2 × 4
## sex N obs exp
## <chr> <dbl> <dbl> <dbl>
## 1 1 138 112 91.6
## 2 2 90 53 73.4
broom::glance(fit)
## # A tibble: 1 × 3
## statistic df p.value
## <dbl> <dbl> <dbl>
## 1 10.3 1 0.00131对于不同组别之间生存曲线的检验，也可以通过K-M图示的方法：
fit.logrank <- survfit(Surv(time, status) ~ sex, data = df)
# 这一步输出太多，我注释掉了，可以自己运行看看
# surv_summary(fit.logrank) # 可以查看寿命表通过ggsurvplot()进行可视化，非常多的细节可以修改，超级详细的教程可以参考我的另一篇推文：R语言生存曲线的可视化(超详细)
ggsurvplot(fit.logrank,
data = df,
surv.median.line = "hv", # Add medians survival
# Change legends: title & labels
legend.title = "Sex",
legend.labs = c("Male", "Female"),
# Add p-value and tervals
pval = TRUE, # 这里P值直接写数字也行
conf.int = TRUE,
# Add risk table
risk.table = TRUE,
tables.height = 0.2,
tables.theme = theme_cleantable(),
ncensor.plot = TRUE,
# Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
# or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
palette = c("#E7B800", "#2E9FDF"),
ggtheme = theme_bw(), # Change ggplot2 theme
# Change font size, style and color
main = "Survival curve",
font.main = c(16, "bold", "darkblue"),
font.x = c(14, "bold.italic", "red"),
font.y = c(14, "bold.italic", "darkred"),
font.tickslab = c(12, "plain", "darkgreen")
)
## Warning in geom_segment(aes(x = 0, y = max(y2), xend = max(x1), yend = max(y2)), : All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
自带的surv_cutpoint()可用于寻找最佳切点，但是只能用于连续性数据。
使用myeloma数据进行演示。
rm(list = ls())
# 0. Load some data
data(myeloma)
head(myeloma)
## molecular_group chr1q21_status treatment event time CCND1 CRIM1
## GSM50986 Cyclin D-1 3 copies TT2 0 69.24 9908.4 420.9
## GSM50988 Cyclin D-2 2 copies TT2 0 66.43 16698.8 52.0
## GSM50989 MMSET 2 copies TT2 0 66.50 294.5 617.9
## GSM50990 MMSET 3 copies TT2 1 42.67 241.9 11.9
## GSM50991 MAF <NA> TT2 0 65.00 472.6 38.8
## GSM50992 Hyperdiploid 2 copies TT2 0 65.20 664.1 16.9
## DEPDC1 IRF4 TP53 WHSC1
## GSM50986 523.5 16156.5 10.0 261.9
## GSM50988 21.1 16946.2 1056.9 363.8
## GSM50989 192.9 8903.9 1762.8 10042.9
## GSM50990 184.7 11894.7 946.8 4931.0
## GSM50991 212.0 7563.1 361.4 165.0
## GSM50992 341.6 16023.4 2096.3 569.2寻找最佳切点：
# 1. Determine the optimal cutpoint of variables
res.cut <- surv_cutpoint(myeloma, time = "time", event = "event",
variables = c("DEPDC1", "WHSC1", "CRIM1") # 找这3个变量的最佳切点
)
summary(res.cut)
## cutpoint statistic
## DEPDC1 279.8 4.275452
## WHSC1 3205.6 3.361330
## CRIM1 82.3 1.968317查看根据最佳切点进行分组后的数据分布情况：
# 2. Plot cutpoint for DEPDC1
plot(res.cut, "DEPDC1", palette = "npg")
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
## Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
## Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## $DEPDC1
根据最佳切点重新划分数据，这样数据就根据最佳切点变成了高表达/低表达组。
# 3. Categorize variables
res.cat <- surv_categorize(res.cut)
head(res.cat)
## time event DEPDC1 WHSC1 CRIM1
## GSM50986 69.24 0 high low high
## GSM50988 66.43 0 low low low
## GSM50989 66.50 0 low high high
## GSM50990 42.67 1 low high low
## GSM50991 65.00 0 low low low
## GSM50992 65.20 0 high low low根据最佳切点绘制生存曲线：
# 4. Fit survival curves and visualize
library("survival")
fit <- survfit(Surv(time, event) ~DEPDC1, data = res.cat)
ggsurvplot(fit, data = res.cat, risk.table = TRUE, conf.int = TRUE)
## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## Warning: tables.theme is not a valid theme.
## Please use `theme()` to construct themes.
确定最佳切点的R包还有非常多，公众号医学和生信笔记后台回复最佳截点即可获取全集。
24.3 Cox回归
本文不涉及理论，只有实操，想要了解生存分析的理论的请自行学习。
使用survival包中的lung数据集用于演示，这是一份关于肺癌患者的生存数据。time是生存时间，以天为单位，status是生存状态，1代表删失，2代表死亡。
rm(list = ls())
library(survival)
library(survminer)
str(lung)
## 'data.frame': 228 obs. of 10 variables:
## $ inst : num 3 3 3 5 1 12 7 11 1 7 ...
## $ time : num 306 455 1010 210 883 ...
## $ status : num 2 2 1 2 2 1 2 2 2 2 ...
## $ age : num 74 68 56 57 60 74 68 71 53 61 ...
## $ sex : num 1 1 1 1 1 1 2 2 1 1 ...
## $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
## $ meal.cal : num 1175 1225 NA 1150 NA ...
## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...可以使用cox回归探索危险因素。分类变量需要变为因子型，这样在进行回归时会自动进行哑变量设置。
lung$sex <- factor(lung$sex, labels = c("female","male"))
lung$ph.ecog <- factor(lung$ph.ecog, labels = c("asymptomatic", "symptomatic",
'in bed <50%','in bed >50%'))
str(lung)
## 'data.frame': 228 obs. of 10 variables:
## $ inst : num 3 3 3 5 1 12 7 11 1 7 ...
## $ time : num 306 455 1010 210 883 ...
## $ status : num 2 2 1 2 2 1 2 2 2 2 ...
## $ age : num 74 68 56 57 60 74 68 71 53 61 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 2 2 1 1 ...
## $ ph.ecog : Factor w/ 4 levels "asymptomatic",..: 2 1 1 2 1 2 3 3 2 3 ...
## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
## $ meal.cal : num 1175 1225 NA 1150 NA ...
## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...拟合多因素Cox回归模型，这里我们只用sex/age两个变量做演示：
fit.cox <- coxph(Surv(time, status) ~ sex + age + ph.karno, data = lung)
# 查看结果
summary(fit.cox)
## Call:
## coxph(formula = Surv(time, status) ~ sex + age + ph.karno, data = lung)
##
## n= 227, number of events= 164
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sexmale -0.497170 0.608249 0.167713 -2.964 0.00303 **
## age 0.012375 1.012452 0.009405 1.316 0.18821
## ph.karno -0.013322 0.986767 0.005880 -2.266 0.02348 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sexmale 0.6082 1.6441 0.4378 0.8450
## age 1.0125 0.9877 0.9940 1.0313
## ph.karno 0.9868 1.0134 0.9755 0.9982
##
## Concordance= 0.637 (se = 0.025 )
## Likelihood ratio test= 18.81 on 3 df, p=3e-04
## Wald test = 18.73 on 3 df, p=3e-04
## Score (logrank) test = 19.05 on 3 df, p=3e-04结果解读和logistic回归的结果解读类似：R语言logistic回归的细节解读
coef是回归系数，exp(coef)是HR值，se(coef)是回归系数的标准误，z是Wald检验的z值，Pr(>|z|)是回归系数的P值，lower .95/upper .95是HR值的95%可信区间。
Concordance= 0.645是Cox回归的C-index，最后给出了Likelihood ratio test似然比检验的统计量、自由度、P值；Wald test的统计量、自由度、P值；Score (logrank) test的统计量、自由度、P值。
想获得整洁的结果不需要自己提取，只要用神包broom即可：
broom::tidy(fit.cox, exponentiate = T, conf.int = T)
## # A tibble: 3 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 sexmale 0.608 0.168 -2.96 0.00303 0.438 0.845
## 2 age 1.01 0.00940 1.32 0.188 0.994 1.03
## 3 ph.karno 0.987 0.00588 -2.27 0.0235 0.975 0.998estimate：HR值（exp(coef)）std.error：回归系数的标准误（se(coef)）statistic：Wald检验的z值p.value：回归系数的P值conf.low/conf.high：HR的95%的可信区间
构建好Cox回归后，也可以用函数单独提取想要的结果，以下图片展示了可用于提取模型信息的函数，和logistic回归差不多：
进行Cox回归必须要符合等比例风险假设，关于什么是等比例风险假设，可以参考郑老师的这篇文章：生存分析COX回归，小心你的数据不符合应用条件
等比例风险的检验可以通过很多方法进行，比如K-M曲线，一般如果有交叉，那么可能不符合等比例风险假设，还可以通过各种残差分布来检验。
下面是Cox回归的等比例风险假设检验，检验方法是基于Schoenfeld残差：
ftest <- cox.zph(fit.cox)
ftest
## chisq df p
## sex 3.085 1 0.0790
## age 0.478 1 0.4892
## ph.karno 8.017 1 0.0046
## GLOBAL 10.359 3 0.0157可以看到ph.karno的P值是小于0.05的，其实是不满足等比例风险假设的，下一篇推文会说到不符合等比例风险假设时该怎么办。
这种方法是基于Schoenfeld残差，检验结果可以通过图示画出来：
library(survminer)
ggcoxzph(ftest)
## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
可以看到sex和age的回归系数随着时间变化基本没啥变化，稳定在0水平线上，和上面的检验结果一样。
还可以通过以下方式查看残差的变化：
ggcoxdiagnostics(fit.cox, type = "schoenfeld")
这张图反映的也是回归系数随时间的变化趋势，和上面的图意思一样，如果符合比例风险假设，那么结果应该是一条水平线，从图示来看，这3个变量都是有点问题的，但是真实数据往往不可能是完美的，很少有完全符合要求的数据。
除了Schoenfeld残差外，ggcoxdiagnostics()还支持其他类型，比如：“martingale”, “deviance”, “score”,“dfbetas” and “scaledsch”在，只需要在type参数中提供合适的类型即可。
cox回归也是回归分析的一种，可以计算出回归系数和95%的可信区间，因此结果可以通过森林图展示：
# 为了森林图好看点，多选几个变量
fit.cox <- coxph(Surv(time, status) ~ . , data = lung)
ggforest(fit.cox, data = lung,
main = "Hazard ratio",
cpositions = c(0.01, 0.15, 0.35), # 更改前三列的相对位置
fontsize = 0.7,
refLabel = "reference",
noDigits = 2
)
这个结果如果你觉得不好看，或者你还有其他的森林图想做到统一的样式，可以考虑我公众号中介绍的画森林图的方法进行个性化定制:
以上是Cox回归的主要内容，大家有问题可以加群或者评论区留言。
24.4 时间依存协变量的Cox回归和时间依存系数Cox回归
之前分别介绍了生存分析中的寿命表法、K-M曲线、logrank检验，以及Cox回归的构建、可视化以及比例风险检验的内容。
本次主要介绍如果数据不符合PH假设时采取的方法。
关于时依协变量、时依系数的基础知识，大家可以参考这几篇文章：
- survival包的案例介绍：Using Time Dependent Covariates and Time Dependent Coefcients in the Cox Model
- 医咖会：一文详解时依协变量
- 7code：含时依协变量的Cox回归
如果不能满足PH假设，可以考虑使用时依协变量或者时依系数Cox回归，时依协变量和时依系数是两个概念，简单来说就是如果一个协变量本身会随着时间而改变，这种叫时依协变量，如果是协变量的系数随着时间改变，这种叫时依系数。
这里以survival包的veteran数据集为例，演示如何处理此类不符合PH检验的情况。
rm(list = ls())
library(survival)
str(veteran)
## 'data.frame': 137 obs. of 8 variables:
## $ trt : num 1 1 1 1 1 1 1 1 1 1 ...
## $ celltype: Factor w/ 4 levels "squamous","smallcell",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ time : num 72 411 228 126 118 10 82 110 314 100 ...
## $ status : num 1 1 1 1 1 1 1 1 1 0 ...
## $ karno : num 60 70 60 60 70 20 40 80 50 70 ...
## $ diagtime: num 7 5 3 9 11 5 10 29 18 6 ...
## $ age : num 69 64 38 63 65 49 69 68 43 70 ...
## $ prior : num 0 10 0 10 10 0 10 0 0 0 ...这个数据集中的变量解释如下图：
首先构建普通的Cox回归，进行等比例风险假设，这里只选择了trt/prior/karno3个变量，而且trt/prior作为分类变量并没有转换为因子型，因为二分类变量数值型和因子型的结果是一样的，转不转换没啥影响！
fit <- coxph(Surv(time, status) ~ trt + prior + karno, data = veteran)
# 进行PH检验
zp <- cox.zph(fit)
zp
## chisq df p
## trt 0.288 1 0.59125
## prior 2.168 1 0.14087
## karno 12.138 1 0.00049
## GLOBAL 18.073 3 0.00042可以看到变量karno的P值小于0.05，是不满足PH假设的。
通过图形化方法查看PH检验的结果：
#op <- par(mfrow=c(1,3))
#plot(zp)
#par(op)
#ggcoxdiagnostics(fit, type = "schoenfeld")
plot(zp[3])
abline(0,0, col="red") # 0水平线
abline(h=fit$coef[3], col="green", lwd=2, lty=2) 
黑色实线以及两侧的虚线是karno的系数随着时间变化的曲线，绿色虚线是假设karno符合PH检验时的总体估计线，红色实线是参考线。
这张图反映了karno变量的系数随着时间的改变，karno偏离的比较厉害（上面注释掉的代码可以都运行看看其他变量的情况），系数最开始接近-0.05，然后逐渐趋于0，最后又开始趋向-0.05，所以它的系数是一致在随着时间改变的，不符合比例风险假设。
24.4.1 对时间分层
这种情况下一个比较简单的解决方式是对时间使用分层函数。根据上面的图示我们知道karno的系数大概分为3层（3段），可以根据两个拐点进行分层，通过survival中的survSplit()实现。
vet2 <- survSplit(Surv(time, status) ~ ., data= veteran,
cut=c(90, 180), # 两个拐点把时间分为3层（3段）
episode= "tgroup",
id="id")
vet2[1:7, c("id", "tstart", "time", "status", "tgroup", "age", "karno")]
## id tstart time status tgroup age karno
## 1 1 0 72 1 1 69 60
## 2 2 0 90 0 1 64 70
## 3 2 90 180 0 2 64 70
## 4 2 180 411 1 3 64 70
## 5 3 0 90 0 1 38 60
## 6 3 90 180 0 2 38 60
## 7 3 180 228 1 3 38 60结果多了两列：tstart/tgroup。
受试者1（id编号为1）在第72天的时候死了，所以数据和之前一样。受试者2和3（id为2和3）虽然时间在变，但是直到第3层才死去，karno的值没有变化。
重新拟合Cox模型，此时tgroup是分好的层，所以要用strata()，另外karno会随着时间变化，和时间有交互，所以用karno:strata(tgroup)。
# 注意此时Surv()的用法！
fit2 <- coxph(Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup), data = vet2)
fit2
## Call:
## coxph(formula = Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup),
## data = vet2)
##
## coef exp(coef) se(coef) z p
## trt -0.011025 0.989035 0.189062 -0.058 0.953
## prior -0.006107 0.993912 0.020355 -0.300 0.764
## karno:strata(tgroup)tgroup=1 -0.048755 0.952414 0.006222 -7.836 4.64e-15
## karno:strata(tgroup)tgroup=2 0.008050 1.008083 0.012823 0.628 0.530
## karno:strata(tgroup)tgroup=3 -0.008349 0.991686 0.014620 -0.571 0.568
##
## Likelihood ratio test=63.04 on 5 df, p=2.857e-12
## n= 225, number of events= 128结果表明karno这个变量只有在tgroup=1（第1层，前3个月）才有意义，后面两层是没有意义的。
再次进行PH检验：
cox.zph(fit2)
## chisq df p
## trt 1.72 1 0.189
## prior 3.81 1 0.051
## karno:strata(tgroup) 3.04 3 0.385
## GLOBAL 8.03 5 0.154这时karno:strata(tgroup)就满足了等比例风险假设。
24.4.2 连续性时依系数变换
除了对时间进行分层外，还有一种解决方法。
上面的图中我们可以看出karno系数随时间变化的曲线明显不是线性的，我们可以通过数据变换把它变成类似线性的，比如取log，这种变换通过tt(time transform)函数实现。
这种方法实际上是通过tt()函数构建了一个时依协变量，但是这样做是为了解决系数随着时间改变的问题（也就是为了解决时依系数的问题）。
fit3 <- coxph(Surv(time, status) ~ trt + prior + karno + tt(karno), # 对karno进行变换
data = veteran,
tt = function(x, t, ...) x * log(t+20) # 具体变换方式
)
fit3
## Call:
## coxph(formula = Surv(time, status) ~ trt + prior + karno + tt(karno),
## data = veteran, tt = function(x, t, ...) x * log(t + 20))
##
## coef exp(coef) se(coef) z p
## trt 0.016478 1.016614 0.190707 0.086 0.93115
## prior -0.009317 0.990726 0.020296 -0.459 0.64619
## karno -0.124662 0.882795 0.028785 -4.331 1.49e-05
## tt(karno) 0.021310 1.021538 0.006607 3.225 0.00126
##
## Likelihood ratio test=53.84 on 4 df, p=5.698e-11
## n= 137, number of events= 128此时karno的时依系数估计为：-0.124662 * log(t + 20)。
在构建时依协变量时，可以选择x * t、x * log(t)、x * log(t + 20)、x * log(t + 200)等等，没有明确的规定，要结合结果和图示进行选择，可以参考冯国双老师的文章：一文详解时依协变量。
我们可以把现在的时依系数估计和经过变换后的的PH检验画在一起，看看变换后的效果：
# 变换后的PH检验
zp <- cox.zph(fit, transform = function(time) log(time + 20))
# 画图
plot(zp[3])
abline(0,0, col="red") # 0水平线
abline(h=fit$coef[3], col="green", lwd=2, lty=2) # 整体估计
abline(coef(fit3)[3:4],lwd=2,lty=3,col="blue") # 现在的估计
可以看到变换后结果好多了（蓝色虚线，和黑色曲线相比较），虽然还是有一点倾斜。
以上是两种处理不满足PH假设的方法，实际还有很多种方法，比较常用的是对时间进行分层，其他方法有机会继续介绍。
24.5 参考资料
- http://www.sthda.com/english/wiki/survival-analysis-basics
- https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html
- survival包帮助文档
- https://mp.weixin.qq.com/s/2rwxeaF_M0UnqPi2F9JNxA