[SAS]存活分析(2) COX regression model

使用存活分析時，除了想要了解不同組別之間是否有顯著差異外，我們也想要了解不同的自變項(X)對存活的影響。正常情形下，我們研究樣本的病人會有不同的特徵，像是性別 、年齡、是否接受特定手術、本身的健康數據等等，我們可以使用 Cox model 知道這些變項對存活函數的影響是什麼。

本篇要介紹的就是Cox proportional hazard model ( Cox比例風險模式 )。在 Cox model中，自變項(X)可以超過兩個以上，並且不管是類別變項或是連續變項都可以。藉著Cox model，我們可以知道各變項的風險，並且調整這些潛在的解釋變數之後，比較各種治療之間的存活時間或存活機率。

example for SAS

和上一篇一樣在我們的研究範圍中追蹤20個病人 ，因為我們想使用Cox model來知道不同類別的自變項對存活函數的影響是甚麼，所以這邊加上了兩個新的變數，第一個是連續資料 weight ，第二個是類別資料 :年齡的分類 age_group ，下面是SAS code 可以直接使用

Data data_cox;
input id time status sex$ weight age_group$ @@;
datalines;
1 121.98680 1 M 58.6 40–49 2 66.33585 0 F 62.4 50–64 3 101.31992 1 M 76.5 40–49 4 120.96599 1 F 56.8 65up 5 30.89584 1 M 66.3 50–64 6 88.88542 0 M 62.7 50–64 7 114.43314 0 F 59.5 65up 8 76.36350 0 M 68.9 50–64 9 139.02362 0 F 74.5 40–49 10 42.17000 0 F 67.5 65up 11 111.80827 1 M 64.1 50–64 12 122.19009 1 F 66.8 40–49 13 80.05458 1 M 81.2 65up 14 107.82935 0 M 64.4 50–64 15 81.77838 0 F 69.4 50–64 16 85.48578 1 M 75.8 65up 17 93.72746 1 F 59.1 50–64 18 69.38994 0 M 63.5 40–49 19 53.35931 0 F 71.5 50–64 20 87.06297 1 M 65.8 40–49
;

單變項的Cox model

我們想要知道每個變項的Cox model，首先使用weight
weight是連續的資料，給定data之後再設定model 
因為是單變項的cox model ，這邊只想要知道weight對存活函數的影響
所以等號前面的第一項先放time，也就是我們的追蹤時間，乘上status，也就是我們有興趣的event ，time*status(0)這兩個組合其實就是我們的依變項(Y)，等號後面放上我們的自變項(X)，也就是weight

proc phreg data=data_cox; /*連續*/
model time*status(0)=weight/RISKLIMITS;run;

跑完之後就可以看到weight的Hazard ratio 和 95% HR CI 以及 p-value
p-value=0.5171 >0.05 ，表示沒有顯著影響的風險因子
(若p-value<0.05 ，表示有顯著影響並且每增加一公斤死亡風險增加1.031倍，注意!這邊並沒有顯著影響，只是舉例顯著時如何解釋風險函數)

class 類別資料

再來是類別資料 age_group，需要使用class 讓SAS知道他是類別，我們也需要設定一個參考組(param=ref ref=’40–49′)，這邊的參考組是40-49歲的病人

proc phreg data=data_cox;
class age_group(param=ref ref=’40–49') ;
model time*status(0)=age_group/RISKLIMITS;run;

因為我們設定了參考組，在class level information 也可以看到 design variables : 40–49 的design variables都是0 ，表示參考組

跑完之後就可以看到weight的Hazard ratio 和 95% HR CI 以及 p-value
p-value >0.05 ，表示沒有顯著影響的風險因子
(若兩個p-value都小於0.05呢 ? 表示有顯著影響
年齡50-64歲的病人的死亡風險是參考組(40-49歲的病人)的2.716倍
年齡65歲以上的病人的死亡風險是參考組(40–49歲的病人)的2.856倍
注意!這邊並沒有顯著影響，只是舉例顯著時如何解釋風險函數)

class SEX

最後一個附上性別的cox model :

proc phreg data=data_cox;
class SEX(param=ref ref=’F’) ;
model time*status(0)=SEX/RISKLIMITS;run;

同樣，p-value = 0.0662 >0.05 ，表示沒有顯著差異

多變項的Cox model

接下來就是多變項的Cox model ，我們不只想要知道單變項對於存活函數的影響，也想知道多個自變項同時對於存活函數的影響。下面是程式的範例

在依變項(Y)的地方(model 等號的右邊)依舊是 time*status(0)
現在的自變項(X)不只是一個了，分別是SEX weight age_group 
若自變項是類別的話，記得加上class 並附上參考組讓SAS判讀
plots(overlay)=survival 則是會畫上經過調整後的KM plot

proc phreg data=data_cox plots(overlay)=survival;
class SEX(param=ref ref=’F’) age_group(param=ref ref=’40–49');
model time*status(0)=SEX weight age_group/RISKLIMITS;run;

結果如下，有沒有發現數字和剛剛單變項不同呢?
可以看到SEX 的 p-value <0.05 ，表示有顯著影響，男性病人的死亡風險是女性的7.971倍。

因為Cox 可以同時校正其他干擾因子，所以畫出來的KM plot 就是調整後的KM plot

KM survival curve

接著，如果想要和上一篇一樣比較兩組的存活率，這邊也有指令可以將KM 存活曲線分成兩組。baseline survival=s /diradj group=SEX;

proc phreg data=data_cox plots(overlay)=survival;
class age_group(param=ref ref=’40–49') SEX(param=ref ref=’F’);
model time*status(0)=age_group weight SEX /RISKLIMITS;
baseline survival=s /diradj group=SEX;
run;

我們將研究的病人依照性別分組，可以看到Hazard Ratio 跟上面是一樣的，加了baseline survival=s /diradj group=SEX; 可以將KM survival plot 分成兩組

如下圖，紅色是男性的KM 存活曲線，藍色是女性的KM存活曲線，都是經過Cox 之後，調整後的KM plot。有發現跟調整前的KM plot不一樣嗎? 可以看看上一篇文章[SAS]存活分析(1) KM存活曲線

如何知道這兩組有沒有顯著差異呢? 可以看到Testing Global Null Hypothesis 的likelihood ratio p-value=0.0637>0.05 ，表示這兩組經過COX調整之後並沒有顯著差異。

存活分析最常使用到的 KM plot 和 Cox model 都在這兩篇講完啦!
接下來想介紹 : 整體存活率 vs 無病存活率