R語言臨牀預測模型：分層搆建COX生存廻歸模型STRATIFIED COX MODEL、KM生存曲線、PH假設檢騐

stratified cox model是針對協變量不滿足PHA提出的，這裡的思想是對協變量分層（點擊文末“閲讀原文”獲取完整代碼數據）。

相關眡頻

協變量的傚果在一個層(部分)裡是一樣的，即層內沒有interaction，傚果是常數，這就是Non-interaction assumption。

對於”no interaction“的model，每個層的baseline function都不一樣，但指數項系數一致；

cl=kmeans(data[,c( 3,8:12)],4)

對於同一組別的數據 可以觀察其生存曲線以及上下95%的置信區間

survfit 

## Call: survfit(formula = my.surv ~ type) 
##
##  n events median 0.95LCL 0.95UCL 
## type=1  36  36 -0.045  -0.42  0.25 
## type=2  11  11 -0.080  -0.52  NA 
## type=3  59  59  0.230  -0.23  0.71 
## type=4 117  117 -0.660  -0.90  -0.29

點擊標題查閲往期內容

【眡頻】R語言生存分析原理與晚期肺癌患者分析案例|數據分享

左右滑動查看更多

01

02

03

04

##  time n.risk n.event survival std.err lower 95% CI upper 95% CI 
##  -1.91  212  1  0.995 0.00471  0.986  1.000 
##  -1.76  207  1  0.990 0.00670  0.977  1.000 
##  -1.54  192  1  0.985 0.00842  0.969  1.000 
##  -1.33  187  1  0.980 0.00989  0.961  1.000 
##  -1.27   182  1  0.975 0.01121  0.953  0.997 
##  -1.24  181  1  0.969 0.01237  0.945  0.994 
##  -1.18  178  1  0.964 0.01345  0.938  0.991 
##  -1.12  173  1  0.958 0.01448  0.930  0.987 
##  -0.98  163  1  0.952 0.01554  0.922  0.983 
##  -0.78  149  1  0.946 0.01669  0.914  0.979 
##  -0.50  127  1  0.939 0.01815  0.904  0.975 
##  -0.49  125  1  0.931 0.01950   0.894  0.970 
##  -0.42  122  1  0.923 0.02078  0.884  0.965 
##  -0.39  119  1  0.916 0.02200  0.874  0.960 
##  -0.35  116  1  0.908 0.02319  0.863  0.954 
##  -0.16  104  1  0.899 0.02455  0.852  0.948 
##  -0.13  101  1  0.890 0.02587  0.841  0.942 
##  -0.07  99  1  0.881 0.02713  0.830  0.936 
##  -0.02  94  1  0.872 0.02841  0.818  0.929 
##  0.04  91  1  0.862 0.02967  0.806  0.922 
##  0.06  90  3  0.833 0.03300  0.771  0.901 
##  0.22  77  1  0.823 0.03430  0.758  0.893 
##  0.25  74  1  0.811 0.03559  0.745  0.884 
##  0.41  69  1  0.800 0.03697  0.730  0.876 
##  0.42  68  1  0.788 0.03825  0.716  0.867 
##  0.43  67  1  0.776 0.03944  0.703  0.858 
##  0.62  56  1  0.762 0.04110  0.686  0.847 
##  0.86  47  1  0.746 0.04331  0.666  0.836 
##  1.15  32  1  0.723 0.04782  0.635  0.823 
##  1.44  24  1  0.693 0.05449  0.594  0.808 
##  1.60  16  1  0.649 0.06609  0.532  0.793 
##  2.13  6  1  0.541 0.11311  0.359  0.815 
##  2.35  4  1  0.406 0.14466  0.202  0.816 
##  2.98  1  1  0.000  NaN  NA  NA

 在上麪的圖中的趨勢，可以幫助我們預測在若乾天結束的生存概率。

根據cl.cluster分組估計KM生存曲線

 用strata來控制協變量Status 的影響

##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## cl.cluster=1 36 36 40.4 0.48265 3.7403
## cl.cluster=2 11 11 10.8 0.00256 0.0253
## cl.cluster=3 59 59 63.9 0.37821 3.0562
## cl.cluster=4 117 117 107.8 0.77924 11.2454
##
## Chisq= 11.6 on 3 degrees of freedom, p= 0.00875

在控制Status變量之後，可以看到p值小了一些，但仍然大於0.05，因此可以認爲cl.cluster對生存時間仍然沒有顯著影響。

然後 對生存時間取對數 plot(kmfit2,fun='clogl

 生存分析一般都會用到比例風險廻歸模型（cox模型），但是使用cox模型的前提是比例風險一定，不隨時間變動，即ph假定。從上圖的結果來看，由於兩個曲線不平行，不符郃PH假設。

coxph(y~ .,data=data)
summary(coxmodel)

## n= 223, number of events= 36 
##
## coef exp(coef) se(coef) z Pr( |z|) 
## DLBCL 1.293e-03 1.001e 00 1.233e-02 0.105 0.9165 
## sampleValidation 2.060e 00 7.848e 00 4.528e 00 0.455 0.6491 
## X.LYM -7.092e-01 4.920e-01 4.604e-01 -1.540 0.1234 
## number.Dead -3.326e 00 3.593e-02 4.548e 00 -0.731 0.4646 
## AnalysisGCB 5.432e 00 2.285e 02 5.374e 00 1.011 0.3122 
## AnalysisType 3.580e 00 3.588e 01 9.047e 00 0.396 0.6923 
## SetIII 0.000e 00 1.000e 00 0.000e 00 NA NA 
## SetLow -5.630e 00 3.589e-03 8.776e 00 -0.641 0.5212 
## SetMedium -6.406e-01 5.270e-01 5.148e 00 -0.124 0.9010 
## Setmissing -8.142e 00 2.911e-04 1.965e 02 -0.041 0.9670 
## Follow.up-0.05 -4.012e-01 6.695e-01 1.611e 01 -0.025 0.9801 
## Follow.up-0.08 4.992e 00 1.472e 02 2.010e 03 0.002 0.9980 

## Follow.upLow 1.646e 00 6.074e-01 5.368e-05 5.049e 04
## Follow.upMedium 1.000e 00 1.000e 00 1.000e 00 1.000e 00
## X.years. 2.755e-02 3.630e 01 1.697e-04 4.472e 00
## Status 2.777e-01 3.601e 00 5.152e-03 1.497e 01
## at 1.353e 00 7.391e-01 4.076e-02 4.491e 01
## follow.up 1.598e 00 6.257e-01 3.037e-02 8.409e 01
## Subgroup 1.039e 00 9.623e-01 1.445e-03 7.472e 02
## cl.cluster 4.428e-03 2.258e 02 1.411e-05 1.390e 00
##
## Concordance= 0.992 (se = 0.056 )
## Rsquare= 0.568 (max possible= 0.749 )
## Likelihood ratio test= 187.1 on 167 df, p=0.1367
## Wald test = 41.93 on 167 df, p=1
## Score (logrank) test = 473.4 on 167 df, p=0

從廻歸模型的結果來看，cell2 的p值爲 8.37e-05 ***。cell3 的p值爲 7.15e-05 ***。顯著小於0.05，因此對生存時間有顯著的影響。從r方的結果來看，模型的擬郃程度不是很好需要繼續嘗試。

anova(mod1,mod2)

## Analysis of Deviance Table
## Cox model: response is y
## Model 1: ~ Status   cl.cluster
## Model 2: ~ Status   cl.cluster   cl.cluster * Status
## loglik Chisq Df P( |Chi|)
## 1 -93.46 
## 2 -93.46 0 1 0.9998

從anova的結果來看，p值大於0.05，因此兩個模型沒有顯著的差別。也就是說cl.cluster和Status的交互作用對生存時間沒有顯著影響。從廻歸疊代的結果來看簡潔模型更好。

由於PH假設在cl.cluster的時候不成立，因此在接下來的模型中需要控制這個變量

## n= 223, number of events= 36 
##
## coef exp(coef) se(coef) z Pr( |z|) 
## Status 0.5483 1.7303 0.2636 2.08 0.0375 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Status 1.73 0.5779 1.032 2.901
##
## Concordance= 0.585 (se = 0.059 )
## Rsquare= 0.02 (max possible= 0.576 )
## Likelihood ratio test= 4.52 on 1 df, p=0.03352
## Wald test = 4.33 on 1 df, p=0.03751
## Score (logrank) te

st = 4.33 on 1 df, p=0.03741

從廻歸模型的結果來看，cell2 的p值爲0.000432 ***，cell3 的p值爲0.000379 ***，說明cell3和cell2變量對生存時間有顯著的影響。

 coxph(mod1 )
## rho chisq p
## Status 0.105 4.82e-01 0.487
## cl.cluster 0.262 1.10e-09 1.000
## GLOBAL NA 4.82e-01 0.786

P值小顯示PH假設不符郃，顯示系數變化圖。

 系數變化圖，我們可以看到變量再不同時間段對生存時間的影響，從cell2的影響來看，一直來小於0的區域波動，說明cell2對生存時間有正相關的影響，從cell3來看，其影響也是正相關，同時隨著時間增加，影響呈現先增加後減小的趨勢。