在之前的文章中，筆者介紹了Cox proportional hazard model（Cox 比例危險模式，接下來簡稱 Cox model）關於比例危險假設（proportional hazard assumption; 以下簡稱 PH assumption）的兩種檢定方式。第一種為檢視個別解釋變項的 Schoenfeld residual 與遞增排序後的事件存活時間的相關性。第二種為在 Cox model 中，增加解釋變項與存活時間的交互作用項，通常存活時間會先經過自然對數的轉換，亦即取 ln (survival time)。解釋變項與存活時間的交互作用項的作法很簡單，只要在資料中先創造好解釋變項與存活時間（自然對數轉換過後）的交互作用變項。或甚至是不需要在資料中另外創造變項，而是直接在統計軟體的語法中增加交互作用變項，通常是用「*」這個符號來將兩個變項作交乘的動作。因此，筆者預計在此篇文章示範如何以 R 統計軟體進行 Schoenfeld residual 檢定方式的操作與解讀。首先，會需要使用到兩個套件，分別為「survival」跟「survminer」。
安裝套件後，先呼叫，然後將資料匯入。示範資料為第四期頭頸癌病人的數據，一共 159 名病人。變項說明如下表：
首先使用「coxph」指令（survival 套件的指令）進行 Cox PH model，解釋變項包含 Age、Stage_4gr 以及 ECS。其中由於 Stage_4gr 是多個水準的類別變項，因此在方程式中要宣告「factor (Stage_4gr)」，並將該模式的結果指定為「fit1」的物件。
以下是 fit1 的報表內容
一、scaled Schoenfeld residual 的檢驗
接著使用「cox.zph」指令（survival 套件的指令）進行 scaled Schoenfeld residual 的檢驗，並儲存為「test.ph」的物件。目前裡頭的設定都是套件的預設，其中「transform」表示對於存活時間的轉換方式，有三種選項，分別為 “km”（預設）、”rank”（排序）以及 ”identity”（原本的存活時間）。「terms」預設為 TRUE，表示如果 k 個水準（k 大於或等於 3）的類別型解釋變項時，列出聯合檢定的結果（自由度為 k-1），而不是列出個別迴歸係數（例如虛擬變數）的檢定結果（自由度為 1）。此外，當指令下 TRUE 的時候，之後顯示的圖形也會以線性預測值（linear predictor）呈現，即使該類別型解釋變項有 k 個水準，也只會列出一條線。「singledf」是指是否將 k 個水準（k 大於或等於 3）的自由度由原本的 k-1 變成 1，預設是 FALSE；「global」是指整個模式所包含的全部解釋變項的聯合檢定，預設是 TRUE。接著將物件列印出來：print(test.ph)
「test.ph」報表內容可看到有三個解釋變項的各自的檢定結果，此卡方檢定是在檢驗該解釋變項的 scaled Schoenfeld residual 是否與存活時間為獨立的分佈，亦即為 Grambsch and Therneau 的檢定法。由以下報表可知 Age 與 Stage_4gr 此兩個變項的卡方檢定未達顯著，注意 Stage_4gr 的自由度是 3，這是因為我們稍早之前在語法設定為「terms=TRUE」，會以聯合檢定的方式呈現。另外 ECS 的卡方檢定則是達顯著，表示 ECS 的 PH assumption 可能是不成立的。
接下來執行「plot(test.ph)」以繪製圖形，如以下三張圖所示。Y 軸是每一個解釋變項的 scaled Schoenfeld residual，X 軸是存活時間。由圖可知 Age 與 Stage_4gr 的 scaled Schoenfeld residual 大致與存活時間呈現無相關的分佈型態，然而 ECS 的分佈則明顯是一個負相關的型態。
接著我們將「terms=FALSE」，從聯合檢定改成單獨效果（各個虛擬變項）的檢定。
可注意到「Stage_4gr」從本來自由度為 3 的聯合檢定，目前變成三個虛擬變項各自自由度為 1 的檢定，可明顯看到是「factor(Stage_4gr)3」這一項的檢定特別顯著，此為第三組比上第一組的結果（T4N23 vs. T4N0）。
圖形也會變成每一個虛擬變項皆有一張圖，不過從三張虛擬變項各自的圖來看，似乎不是很明顯有明顯的正相關或負相關的型態。因此，應需以聯合檢定的結果為主。
也可以使用「survminer」套件來達成 scaled Schoenfeld residual 的檢驗，使用「cox.zph」的指令即可，會得到與上述完全一樣的檢定結果。接著用「ggcoxzph」繪製圖形，先只繪製「ECS」的圖形，否則一次匯出三張會在同一個視窗堆疊，結果如下圖。
上述內容即為在 R 實現 scaled Schoenfeld residual 的檢驗，以 ECS 這個變項而言，沒有通過 PH assumption 的檢驗，通常有幾種處理方式，最簡單的是分層分析（stratum by ECS），另外一種則是將 ECS 視為時間相依（time-dependent）的變數，或是使用其他延伸的 Cox model，例如 restricted mean survival time（RMST）等替代模型。
未來有機會筆者會再針對違反 PH assumption 時的因應處理方式，有興趣的讀者可以先參考此篇文章（https://ekja.org/upload/pdf/kja-19183.pdf）。
留言列表
{{ article.title }}