之前分別介紹了生存分析中的壽命表法,、K-M曲線,、logrank檢驗(yàn)：R語言生存分析的實(shí)現(xiàn)
以及Cox回歸的構(gòu)建、可視化以及比例風(fēng)險(xiǎn)檢驗(yàn)的內(nèi)容：R語言生存分析：Cox回歸
本次主要介紹如果數(shù)據(jù)不符合PH假設(shè)時(shí)采取的方法,。
時(shí)間依存協(xié)變量的Cox回歸和時(shí)間依存系數(shù)Cox回歸
關(guān)于時(shí)依協(xié)變量,、時(shí)依系數(shù)的基礎(chǔ)知識，大家可以參考這幾篇文章：
- survival包的案例介紹：Using Time Dependent Covariates and Time Dependent
Coefcients in the Cox Model[1]
- 醫(yī)咖會：一文詳解時(shí)依協(xié)變量[2]
- 7code：含時(shí)依協(xié)變量的Cox回歸[3]
如果不能滿足PH假設(shè),，可以考慮使用時(shí)依協(xié)變量或者時(shí)依系數(shù)Cox回歸,，時(shí)依協(xié)變量和時(shí)依系數(shù)是兩個(gè)概念，簡單來說就是如果一個(gè)協(xié)變量本身會隨著時(shí)間而改變,，這種叫時(shí)依協(xié)變量,，如果是協(xié)變量的系數(shù)隨著時(shí)間改變，這種叫時(shí)依系數(shù),。
這里以survival包的veteran數(shù)據(jù)集為例，演示如何處理此類不符合PH檢驗(yàn)的情況,。
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 ...
這個(gè)數(shù)據(jù)集中的變量解釋如下圖：
veteran首先構(gòu)建普通的Cox回歸,，進(jìn)行等比例風(fēng)險(xiǎn)假設(shè)，這里只選擇了trt/prior/karno3個(gè)變量,，而且trt/prior作為分類變量并沒有轉(zhuǎn)換為因子型,，因?yàn)槎诸愖兞繑?shù)值型和因子型的結(jié)果是一樣的，轉(zhuǎn)不轉(zhuǎn)換沒啥影響,！如果你還不懂分類變量在r語言中的編碼方案,，一定要看這篇：分類變量進(jìn)行回歸分析時(shí)的編碼方案
fit <- coxph(Surv(time, status) ~ trt + prior + karno, data = veteran)
# 進(jìn)行PH檢驗(yàn)
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假設(shè)的,。
通過圖形化方法查看PH檢驗(yàn)的結(jié)果：
#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)
黑色實(shí)線以及兩側(cè)的虛線是karno的系數(shù)隨著時(shí)間變化的曲線,，綠色虛線是假設(shè)karno符合PH檢驗(yàn)時(shí)的總體估計(jì)線，紅色實(shí)線是參考線,。
這張圖反映了karno變量的系數(shù)隨著時(shí)間的改變,，karno偏離的比較厲害（上面注釋掉的代碼可以都運(yùn)行看看其他變量的情況），系數(shù)最開始接近-0.05,，然后逐漸趨于0,，最后又開始趨向-0.05,，所以它的系數(shù)是一直在隨著時(shí)間改變的，不符合比例風(fēng)險(xiǎn)假設(shè),。
對時(shí)間分層
這種情況下一個(gè)比較簡單的解決方式是對時(shí)間使用分層函數(shù),。根據(jù)上面的圖示我們知道karno的系數(shù)大概分為3層（3段），可以根據(jù)兩個(gè)拐點(diǎn)進(jìn)行分層,，通過survival中的survSplit()實(shí)現(xiàn),。
vet2 <- survSplit(Surv(time, status) ~ ., data= veteran,
cut=c(90, 180), # 兩個(gè)拐點(diǎn)把時(shí)間分為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
結(jié)果多了兩列：tstart/tgroup。
受試者1（id編號為1）在第72天的時(shí)候死了,，所以數(shù)據(jù)和之前一樣,。受試者2和3（id為2和3）雖然時(shí)間在變，但是直到第3層才死去,，karno的值沒有變化,。
重新擬合Cox模型，此時(shí)tgroup是分好的層,，所以要用strata(),，另外karno會隨著時(shí)間變化，和時(shí)間有交互,，所以用karno:strata(tgroup),。
# 注意此時(shí)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
結(jié)果表明karno這個(gè)變量只有在tgroup=1（第1層,，前3個(gè)月）才有意義,，后面兩層是沒有意義的。
再次進(jìn)行PH檢驗(yàn)：
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
這時(shí)karno:strata(tgroup)就滿足了等比例風(fēng)險(xiǎn)假設(shè),。
連續(xù)性時(shí)依系數(shù)變換
除了對時(shí)間進(jìn)行分層外,，還有一種解決方法。
上面的圖中我們可以看出karno系數(shù)隨時(shí)間變化的曲線明顯不是線性的,，我們可以通過數(shù)據(jù)變換把它變成類似線性的,，比如取log，這種變換通過tt(time transform)函數(shù)實(shí)現(xiàn),。
這種方法實(shí)際上是通過tt()函數(shù)構(gòu)建了一個(gè)時(shí)依協(xié)變量,，但是這樣做是為了解決系數(shù)隨著時(shí)間改變的問題（也就是為了解決時(shí)依系數(shù)的問題）。
fit3 <- coxph(Surv(time, status) ~ trt + prior + karno + tt(karno), # 對karno進(jìn)行變換
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
此時(shí)karno的時(shí)依系數(shù)估計(jì)為：-0.124662 * log(t + 20),。
在構(gòu)建時(shí)依協(xié)變量時(shí),，可以選擇x * t、x * log(t),、x * log(t + 20),、x * log(t + 200)等等，沒有明確的規(guī)定，要結(jié)合結(jié)果和圖示進(jìn)行選擇,，可以參考馮國雙老師的文章：一文詳解時(shí)依協(xié)變量[4],。
我們可以把現(xiàn)在的時(shí)依系數(shù)估計(jì)和經(jīng)過變換后的的PH檢驗(yàn)畫在一起，看看變換后的效果：
# 變換后的PH檢驗(yàn)
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) # 整體估計(jì)
abline(coef(fit3)[3:4],lwd=2,lty=3,col="blue") # 現(xiàn)在的估計(jì)
可以看到變換后結(jié)果好多了（藍(lán)色虛線,，和黑色曲線相比較）,，雖然還是有一點(diǎn)傾斜。
以上是兩種處理不滿足PH假設(shè)的方法,，實(shí)際還有很多種方法,，比較常用的是對時(shí)間進(jìn)行分層，其他方法有機(jī)會繼續(xù)介紹,！
參考資料
[1]survival包: https://cran./web/packages/survival/vignettes/timedep.pdf
[2]醫(yī)咖會時(shí)依協(xié)變量: https://mp.weixin.qq.com/s/CSwO42ucfeKipHl32NoOCA
[3]7code時(shí)依協(xié)變量: https://mp.weixin.qq.com/s/sBJ8Xg_e3uBfUnDOIPoxug
[4]醫(yī)咖會時(shí)依協(xié)變量: https://mp.weixin.qq.com/s/CSwO42ucfeKipHl32NoOCA
