在本文中，我們將描述Cox回歸模型，並提供使用R軟件的實例。
需要進行多元統計建模
在臨床研究中，有許多情況下，幾個已知量（稱為協變量）可能影響患者。
統計模型是一個經常使用的工具，可以同時分析多個因素的生存情況。另外，統計模型提供了每個因素的效應大小。
Cox比例風險模型的基礎
該模型的目的是同時評估幾個因素對生存的影響。換句話說，它使我們能夠檢查特定因素在特定時間點如何影響特定事件（例如，感染，死亡）的發生率。這個速度通常被稱為危險率。預測變量（或因子）在生存分析文獻中通常被稱為協變量。
Cox模型由h（t）表示的危險函數表示。簡而言之，危險函數可以解釋為在時間t死亡的風險。可以估計如下：
h（t）=h0（t）×e×p（b1X1+b2X2+。。。+bpXp）H（Ť）=H0（Ť）×ËXp（b1X1+b2X2+。。。+bpXp）
換言之，高於1的風險比指示與事件概率正相關的協變量，並因此與生存期的長短負相關。
綜上所述，
HR = 1：沒有效果
HR<1：減少危險
HR> 1：危害增加
請注意，在癌症研究中：
危險比> 1（即：b> 0）的協變量被稱為不良預后因素
危險比<1（即：b <0）的協變量被稱為良好的預后因子
Cox模型的一個關鍵假設是，觀察組（或患者）的危險曲線應該是成比例的並且不能交叉。
因此，考克斯模型是一個比例 - 危險模型：任何一組中事件的危險性是其他危險的常數倍數。這一假設意味着，如上所述，各組的危險曲線應該是成比例的，不能交叉。
換句話說，如果一個人在某個初始時間點有死亡風險，是另一個人的兩倍，那么在所有的晚些時候，死亡風險仍然是兩倍。
計算R中的Cox模型
安裝並加載所需的R包
我們將使用兩個R包：
生存計算生存分析
幸存者可視化生存分析結果
安裝軟件包
install.packages(c("survival","survminer"))
加載包
library("survival")library("survminer")
R函數來計算Cox模型：coxph（）
函數coxph（）[在生存包中]可用於計算R中的Cox比例風險回歸模型。
簡化格式如下：
coxph(formula,data,method)
公式：以生存對象為響應變量的線性模型。Survival對象是使用Surv（）函數創建的，如下所示：Surv（time，event）。
數據：包含變量的數據框
方法：用於指定如何處理關系。默認是'efron'。其他選項是“breslow”和“確切”。默認的“efron”通常比一度流行的“breslow”方法更受歡迎。“確切”的方法計算密集得多。
示例數據集
我們將在生存R軟件包中使用肺癌數據。
data("lung")head(lung)
inst：機構代碼
時間：以天為單位的生存時間
狀態：審查狀態1 =審查，2 =死亡
年齡：年齡
性別：男= 1女= 2
ph.ecog：ECOG表現評分（0 =好5 =死）
ph.karno：Karnofsky表現評分（bad = 0-好= 100）由醫師評定
pat.karno：Karnofsky表現評分由患者評估
膳食：餐時消耗的卡路里
wt.loss：過去六個月的體重下降
計算Cox模型
我們將使用以下協變量進行Cox回歸：年齡，性別，ph.ecog和wt.loss。
我們首先計算所有這些變量的單變量Cox分析;那么我們將使用兩個變量來擬合多變量cox分析來描述這些因素如何共同影響生存。
單變量Cox回歸
單變量Cox分析可以計算如下：
res.cox<-coxph(Surv(time,status)~sex,data=lung)res.cox
Call:coxph(formula = Surv(time, status) ~ sex, data = lung)coef exp(coef) se(coef) z psex -0.531 0.588 0.167 -3.18 0.0015Likelihood ratio test=10.6 on 1 df, p=0.00111n= 228, number of events= 165
Cox模型的函數summary（）產生更完整的報告：
summary(res.cox)
Call:coxph(formula = Surv(time, status) ~ sex, data = lung)n= 228, number of events= 165coef exp(coef) se(coef) z Pr(>|z|)sex -0.5310 0.5880 0.1672 -3.176 0.00149 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1exp(coef) exp(-coef) lower .95 upper .95sex 0.588 1.701 0.4237 0.816Concordance= 0.579 (se = 0.022 )Rsquare= 0.046 (max possible= 0.999 )Likelihood ratio test= 10.63 on 1 df, p=0.001111Wald test = 10.09 on 1 df, p=0.001491Score (logrank) test = 10.33 on 1 df, p=0.001312
Cox回歸結果可以解釋如下：
統計顯着性。標記為“z”的列給出了Wald統計值。它對應於每個回歸系數與其標准誤差的比率（z = coef / se（coef））。wald統計評估是否beta（ββ）系數在統計上顯着不同於0。從上面的輸出，我們可以得出結論，變量性別具有高度統計學意義的系數。
回歸系數。Cox模型結果中要注意的第二個特征是回歸系數（coef）的符號。一個積極的信號意味着危險（死亡風險）較高，因此對於那些變量值較高的受試者，預后更差。變量性被編碼為數字向量。1：男，2：女。Cox模型的R總結給出了第二組相對於第一組，即女性與男性的風險比（HR）。性別的β系數= -0.53表明在這些數據中，女性的死亡風險（低存活率）低於男性。
危害比例。指數系數（exp（coef）= exp（-0.53）= 0.59）也稱為風險比，給出協變量的效應大小。例如，女性（性別= 2）將危害降低了0.59倍，即41％。女性與預后良好相關。
風險比的置信區間。總結結果還給出了風險比（exp（coef））的95％置信區間的上限和下限，下限95％界限= 0.4237，上限95％界限= 0.816。
全球統計學意義的模型。最后，輸出為模型的總體顯着性提供了三個替代測試的p值：可能性比率測試，Wald測試和得分logrank統計。這三種方法是漸近等價的。對於足夠大的N，他們會得到相似的結果。對於小N來說，它們可能有所不同。似然比檢驗對於小樣本量具有更好的表現，所以通常是優選的。
要一次性將單變量coxph函數應用於多個協變量，請輸入：
covariates<-c("age","sex","ph.karno","ph.ecog","wt.loss")univ_formulas<-sapply(covariates,function(x)as.formula(paste('Surv(time, status)~',x)))univ_models<-lapply(univ_formulas,function(x){coxph(x,data=lung)})# Extract datauniv_results<-lapply(univ_models,function(x){x<-summary(x)p.value<-signif(x$wald["pvalue"],digits=2)wald.test<-signif(x$wald["test"],digits=2)beta<-signif(x$coef[1],digits=2);#coeficient betaHR<-signif(x$coef[2],digits=2);#exp(beta)HR.confint.lower<-signif(x$conf.int[,"lower .95"],2)HR.confint.upper<-signif(x$conf.int[,"upper .95"],2)HR<-paste0(HR," (",HR.confint.lower,"-",HR.confint.upper,")")res<-c(beta,HR,wald.test,p.value)names(res)<-c("beta","HR (95% CI for HR)","wald.test","p.value")return(res)#return(exp(cbind(coef(x),confint(x))))})res<-t(as.data.frame(univ_results,check.names=FALSE))as.data.frame(res)
beta HR (95% CI for HR) wald.test p.valueage 0.019 1 (1-1) 4.1 0.042sex -0.53 0.59 (0.42-0.82) 10 0.0015ph.karno -0.016 0.98 (0.97-1) 7.9 0.005ph.ecog 0.48 1.6 (1.3-2) 18 2.7e-05wt.loss 0.0013 1 (0.99-1) 0.05 0.83
上面的輸出顯示了與總體生存相關的每個變量的回歸β系數，效應大小（作為危害比給出）和統計顯着性。每個因素通過單獨的單變量Cox回歸評估。
從上面的輸出中，
變量性別，年齡和ph.ecog具有高度統計學意義的系數，而ph.karno的系數不顯着。
年齡和ph.ecog有正的beta系數，而性別有負系數。因此，年齡越大和ph.ecog越高，生存率越差，而女性（性別= 2）與更好的生存相關。
現在我們要描述這些因素如何共同影響生存。為了回答這個問題，我們將執行多元Cox回歸分析。由於變量ph.karno在單變量Cox分析中不顯着，我們將在多變量分析中跳過它。我們將3個因素（性別，年齡和ph.ecog）納入多變量模型。
多變量Cox回歸分析
時間常數協變量的死亡時間的Cox回歸如下所示：
res.cox<-coxph(Surv(time,status)~age+sex+ph.ecog,data=lung)summary(res.cox)
Call:coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)n= 227, number of events= 164(1 observation deleted due to missingness)coef exp(coef) se(coef) z Pr(>|z|)age 0.011067 1.011128 0.009267 1.194 0.232416sex -0.552612 0.575445 0.167739 -3.294 0.000986 ***ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1exp(coef) exp(-coef) lower .95 upper .95age 1.0111 0.9890 0.9929 1.0297sex 0.5754 1.7378 0.4142 0.7994ph.ecog 1.5900 0.6289 1.2727 1.9864Concordance= 0.637 (se = 0.026 )Rsquare= 0.126 (max possible= 0.999 )Likelihood ratio test= 30.5 on 3 df, p=1.083e-06Wald test = 29.93 on 3 df, p=1.428e-06Score (logrank) test = 30.5 on 3 df, p=1.083e-06
所有三個整體測試（可能性，Wald和得分）的p值都是顯着的，表明該模型是顯着的。這些測試評估了所有的beta（ββ）為0.在上面的例子中，檢驗統計是完全一致的，綜合零假設被完全拒絕。
在多變量Cox分析中，協變量性別和ph.ecog仍然顯着（p<0.05）。然而，協變量年齡並不顯着（P = 0.23，這比0.05）。
性別p值為0.000986，危險比HR = exp（coef）= 0.58，表明患者性別和死亡風險降低之間有很強的關系。協變量的風險比可以解釋為對風險的倍增效應。例如，保持其他協變量不變，女性（性別= 2）將危害降低0.58倍，即42％。我們的結論是，女性與良好的預后相關。
類似地，ph.ecog的p值是4.45e-05，危險比HR = 1.59，表明ph.ecog值與死亡風險增加之間的強關系。保持其他協變量不變，ph.ecog值越高，生存率越差。
相比之下，年齡的p值現在是p = 0.23。風險比HR = exp（coef）= 1.01，95％置信區間為0.99至1.03。由於HR的置信區間為1，這些結果表明，年齡在調整了ph值和患者性別后對HR的差異的貢獻較小，並且僅趨向顯着性。例如，保持其他協變量不變，額外的年齡會導致每日死亡危險因素為exp（beta）= 1.01或1％，這不是一個重要的貢獻。
可視化估計的生存時間分布
已經將Cox模型擬合到數據中，可以在特定風險組的任何給定時間點可視化預測的存活比例。函數survfit（）估計生存比例，默認情況下為協變量的平均值。
# Plot the baseline survival functionggsurvplot(survfit(res.cox),color="#2E9FDF",ggtheme=theme_minimal())
我們可能希望展示估計的生存如何取決於感興趣的協變量的值。
考慮到，我們要評估性別對估計生存概率的影響。在這種情況下，我們構建一個兩行的新數據框，每個性別值一個;其他協變量被固定為它們的平均值（如果它們是連續變量的話）或者它們的最低水平（如果它們是離散變量的話）。對於虛擬協變量，平均值是數據集中編碼1的比例。這個數據幀通過newdata參數傳遞給survfit（）：
# Create the new datasex_df<-with(lung,data.frame(sex=c(1,2),age=rep(mean(age,na.rm=TRUE),2),ph.ecog=c(1,1)))sex_df
sex age ph.ecog1 1 62.44737 12 2 62.44737 1
# Survival curvesfit<-survfit(res.cox,newdata=sex_df)ggsurvplot(fit,conf.int=TRUE,legend.labs=c("Sex=1","Sex=2"),ggtheme=theme_minimal())
轉自： https://www.sohu.com/a/205899780_826434