Cox模型變量選擇的R軟件實現

韋新星+趙麗棉

【摘要】本文基于R軟件，對Cox模型的變量選擇進行案例分析。結果表明，R軟件不僅能很好地對Cox模型進行變量選擇，操作簡便，準確度高，而且還能繪制特定變量的生存函數圖，供進一步對比分析使用。

【關鍵詞】R Cox模型 變量選擇 生存函數圖

【中圖分類號】G642.1 【文獻標識碼】A 【文章編號】2095-3089（2017）39-0217-02

一、引言

Cox模型是英國統計學家Cox于1972年提出的[1]。Cox模型一經提出，便被廣泛應用于生存數據的處理和分析中。運用Cox模型進行變量選擇，人們不僅可以處理內含刪失數據的生存數據，而且還可以處理不服從正態分布的生存數據，這為生存數據的處理和分析提供了有效的途徑。

然而，在建立Cox模型進行變量選擇的過程中，由于Cox模型隸屬半參數模型，模型表達式較為復雜，若按照公式進行處理，則不僅計算量大，而且還容易出錯。為此，本文在前人研究的基礎上，基于R軟件，對Cox模型的變量選擇進行R軟件的實現研究，以解決上述問題帶來的影響。

二、Cox模型概述

Cox模型，又稱Cox比例風險模型，是一種處理多因素生存數據的有效方法。其一般形式為：

其中，為協變量，為回歸系數，為與無關的基本風險函數。之所以又將Cox模型稱為Cox比例風險模型，是因為對于不同的協變量和，該模型都能滿足一個特點，即風險函數之比與無關，類似于沒有的情況[2]。于是通常用風險函數來描述被觀測對象在某一時間點存活的條件下，在未來時間里死亡的概率。另一方面，個體生存時間大于時間的概率被命名為生存函數。

Cox模型中的回歸系數通常用偏似然函數來估計[4]。而生存分析中的很多檢驗也都是基于似然函數或偏似然函數的漸進性質實現的，如似然比檢驗、Wald檢驗和得分檢驗等，若這些檢驗的統計量值大于臨界值，或檢驗統計量的P值小于檢驗水平，則拒絕原假設，認為回歸系數有顯著影響。

三、R軟件實現

本節從典型的生存案例出發，基于R軟件，實現Cox模型變量選擇的過程。

案例采用72個乳腺癌病人的數據，數據來源于文獻[3]，共包含4個變量：time，index，age和group，分別表示生存時間、刪失情況（1表示死亡，0表示刪失）、年齡和是否用藥（2表示用了抗癌藥，1表示沒用抗癌藥）。其中，年齡為定量數據，是否用藥為定性數據。現要探究患者手術時的年齡與是否用藥對乳腺癌患者術后生存狀況的影響。

經分析，這是一個典型的生存分析問題，宜采用Cox模型來研究age和group這2個協變量對生存時間time的影響。具體操作如下：

1.數據處理

由于文獻[3]的數據為sav格式，故需將其另存為逗號分隔符文件（csv格式）或文本文件（txt格式），然后再用read.csv（）或read.table（）方法進行讀取。在此，我們以另存為到桌面的csv格式為例。步驟：在SPSS的數據視圖中，通過選擇“文件”→“另存為”，然后在“保存類型”下拉列中選擇csv，以文件名“案例一”另存為到桌面即可。

2. R代碼編寫

接下來，運用R軟件進行Cox模型變量選擇。

（1）先在R中安裝survival包。可通過點擊“Packages”→“Load package”，然后在彈出的select one選擇框中選擇survival，然后點擊“OK”。此時R中會出現以下語句：

> local（{pkg <- select.list（sort（.packages（all.available = TRUE）），graphics=TRUE）

+ if（nchar（pkg）） library（pkg， character.only=TRUE）}）

Loading required package： splines

接著，用R代碼library（survival）進行加載即可。

（2）然后在R中輸入R代碼。需要指出的是，R中的所有代碼都區分大小寫，而且其標點符號都需要用半角格式書寫。Cox模型變量選擇的R代碼為：

u=read.csv（"C：/Users/Administrator/Desktop/案例一.csv"，sep="，"）

library（survival）

a=coxph（Surv（time，index）～age+group，data=u）

sa=step（a）

summary（sa）

回車后即可看到Cox模型變量選擇的輸出結果。

此外，運用R軟件還可以進行生存函數圖的繪制。在R中輸入以下R代碼：

a0=survfit（Surv（time，index）～group，data=u，type="kaplan-meier"）

plot（a0，lty=1：2，xlab="time"，ylab="survival fuction"）

title（"comparison"）

legend（"topright"，c（"group=1"，"group=2"），lty=1：2）

回車后即可看到生存函數的對比圖。

3.結果輸出及分析

（1）Cox模型變量選擇的部分輸出結果如下：

Call：

coxph（formula = Surv（time， index） ～ age + group， data = u）

n= 72， number of events= 36endprint

coef exp（coef） se（coef） z Pr（>|z|）

age 0.2076 1.2307 0.0323 6.423 1.34e-10 ***

group -1.5269 0.2172 0.4086 -3.737 0.000186 ***

---

Signif. codes： 0 ‘*** 0.001 ‘** 0.01 ‘* 0.05 ‘. 0.1 ‘ 1

Concordance= 0.883 （se = 0.055 ）

Rsquare= 0.6 （max possible= 0.978 ）

Likelihood ratio test= 65.91 on 2 df， p=4.885e-15

Wald test= 47.6 on 2 df， p=4.608e-11

Score （logrank） test = 60.62 on 2 df， p=6.872e-14

可以看到，age和group這2個協變量的系數分別為0.2076和-1.5269，其中age的系數為正，表明高齡乳腺癌患者的死亡率高于低齡的，而group的系數為負，表明用了抗癌藥的乳腺癌患者的死亡率低于沒用藥的。Cox模型表達式為：

另一方面，為0.6，表明模型的擬合效果較好。而似然比檢驗、Wald檢驗和得分檢驗的P值均小于0.001，可認為回歸系數有顯著影響。

（2）繪制的生存函數對比圖為：

圖1 group的生存函數對比圖

由圖1，我們可以比較用了抗癌藥的乳腺癌患者（group=2）和沒用抗癌藥的乳腺癌患者（group=1）的生存函數的區別。可以看到，用了抗癌藥的乳腺癌患者的生存狀況明顯高于沒用藥的，這再次表明用藥對抗癌有明顯療效。

四、結論

本文基于R軟件，對Cox模型的變量選擇進行研究分析。從案例的編程及處理上看，盡管數據結構不盡相同，但在對Cox模型進行變量選擇時，R軟件操作簡便，實用性強，準確度高，能很好地對Cox模型進行變量選擇。此外，R軟件還可生成特定變量的生存函數圖，不僅直觀形象，還可供進一步對比分析使用。

參考文獻：

[1] COX D R. Regression Models and Life Tables[J]. Journal of Royal Statistical Society， 1972， 34： 187-220.

[2] 吳喜之. 統計學：從數據到結論[M]. 北京： 中國統計出版社， 2013： 200-207.

[3] 李靜萍， 謝邦昌. 多元統計分析方法與應用[M]. 中國人民大學出版社， 2008： 209-217.

[4] 田俊. Cox回歸模型及其參數估計[J]. 福建醫學院學報. 1987， 21（2）： 115-117.

基金項目：河池學院碩士專業學位建設基金課題（2016YT004）；廣西大學生創新創業訓練計劃項目（201610605054）。

作者簡介：韋新星（1990-），女，廣西柳州人，河池學院數學與統計學院助教，碩士，研究方向：應用統計。endprint