R語言如何實現Cox回歸模型

這篇“R語言如何實現Cox回歸模型”文章的知識點大部分人都不太理解，所以小編給大家總結了以下內容，內容詳細，步驟清晰，具有一定的借鑒價值，希望大家閱讀完這篇文章能有所收獲，下面我們一起來看看這篇“R語言如何實現Cox回歸模型”文章吧。

R語言中的Cox模型分析

單變量Cox回歸

summary(res.cox)
res.cox <- coxph(Surv(time, status) ~ sex, data =  lung)
summary(res.cox)

summary的結果：

Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)
  n= 228, number of events= 165 
       coef 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 ‘ ’ 1
    exp(coef) exp(-coef) lower .95 upper .95
sex     0.588      1.701    0.4237     0.816
Concordance= 0.579  (se = 0.021 )
Likelihood ratio test= 10.63  on 1 df,   p=0.001
Wald test            = 10.09  on 1 df,   p=0.001
Score (logrank) test = 10.33  on 1 df,   p=0.001

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來說，它們可能有所不同。似然比檢驗對于小樣本量具有更好的表現，所以通常是優選的。

批量單變量cox分析

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 data 
univ_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 beta
                          HR <-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.value
age       0.019            1 (1-1)       4.1   0.042
sex       -0.53   0.59 (0.42-0.82)        10  0.0015
ph.karno -0.016      0.98 (0.97-1)       7.9   0.005
ph.ecog    0.48        1.6 (1.3-2)        18 2.7e-05
wt.loss  0.0013         1 (0.99-1)      0.05    0.83

多變量Cox回歸

要一次性將單變量coxph函數應用于多個協變量，請輸入：

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.232416    
sex     -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 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
age        1.0111     0.9890    0.9929    1.0297
sex        0.5754     1.7378    0.4142    0.7994
ph.ecog    1.5900     0.6289    1.2727    1.9864

Concordance= 0.637  (se = 0.025 )
Likelihood ratio test= 30.5  on 3 df,   p=1e-06
Wald test            = 29.93  on 3 df,   p=1e-06
Score (logrank) test = 30.5  on 3 df,   p=1e-06

所有3個整體測試（可能性，Wald 和 得分）的p值都是顯著的，表明該模型是顯著的。這些測試評估了所有的beta（ββ）為0.在上面的例子中，檢驗統計是完全一致的，綜合零假設被完全拒絕。

在多變量Cox分析中，協變量性別和ph.ecog仍然顯著（p <0.05）。然而，協變量年齡并不顯著（P = 0.17，這比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（）估計生存比例，默認情況下為協變量的平均值。

繪制生存曲線：

kmfit<-survfit(Surv(time, status) ~ sex, data =  lung)
ggsurvplot(kmfit,
              pval = TRUE, conf.int = TRUE,
              risk.table = TRUE, 
              risk.table.col = "strata", 
              linetype = "strata", 
              surv.median.line = "hv", 
              ggtheme = theme_bw(), 
              palette = c("#E7B800", "#2E9FDF"))

以上就是關于“R語言如何實現Cox回歸模型”這篇文章的內容，相信大家都有了一定的了解，希望小編分享的內容對大家有幫助，若想了解更多相關的知識內容，請關注億速云行業資訊頻道。