R coxph 擬合比例風險回歸模型

說明

擬合 Cox 比例風險回歸模型。使用 Andersen 和 Gill 的計數過程公式合並時間相關變量、時間相關層、每個受試者的多個事件以及其他擴展。

用法

coxph(formula, data=, weights, subset, 
      na.action, init, control, 
      ties=c("efron","breslow","exact"), 
      singular.ok=TRUE, robust, 
      model=FALSE, x=FALSE, y=TRUE, tt, method=ties,
      id, cluster, istate, statedata, nocenter=c(-1, 0, 1), ...)

參數

formula	公式對象，響應位於 ~ 運算符的左側，術語位於右側。響應必須是 Surv 函數返回的生存對象。對於多狀態模型，公式可以是公式列表。
data	data.frame，用於解釋 formula 或 subset 和 weights 參數中指定的變量。
weights	案例權重向量，請參閱下麵的注釋。有關這些內容的全麵討論，請參閱 Therneau 和 Grambsch 所著的書。
subset	表達式，指示在擬合中應使用數據行的哪個子集。默認情況下包括所有觀察結果。
na.action	缺失數據過濾函數。在使用任何子集參數後，這將應用於model.frame。默認為 options()\$na.action 。
init	迭代初始值的向量。所有變量的默認初始值均為零。
control	coxph.control 類的對象，指定迭代限製和其他控製選項。默認為 coxph.control(...) 。
ties	指定平局處理方法的字符串。如果沒有固定的死亡時間，所有方法都是等效的。默認使用 Efron 近似，它在處理綁定死亡時間時更準確，並且計算效率更高。 (但請參閱下文的多狀態模型。)“exact partial likelihood” 相當於條件邏輯模型，並且適用於時間是一小組離散值的情況。
singular.ok	指示如何處理模型矩陣中的共線性的邏輯值。如果是 TRUE ，程序將自動跳過 X 矩陣中作為早期列的線性組合的列。在這種情況下，此類列的係數將為 NA，並且方差矩陣將包含零。對於輔助計算，例如線性預測器，缺失的係數被視為零。
robust	是否應該計算穩健方差。如果出現以下情況，則默認值為 TRUE：存在 cluster 參數，存在非 0 或 1 的情況權重，或者存在包含多個事件的 id 值。
id	標識主題的可選變量名稱。僅當主題在數據中可以有多行並且有多個事件類型時才需要。該變量通常可以在 data 中找到。
cluster	可選變量，用於對觀察結果進行聚類，以實現穩健的方差。如果存在，則意味著 robust 。該變量通常可以在 data 中找到。
istate	可選變量，給出每個間隔開始時的當前狀態。該變量通常可以在 data 中找到。
statedata	用於說明多狀態模型的可選數據集。
model	邏輯值：如果 TRUE ，則模型框架在組件 model 中返回。
x	邏輯值：如果 TRUE ，則 x 矩陣在組件 x 中返回。
y	邏輯值：如果 TRUE ，則響應向量在組件 y 中返回。
tt	time-transform 函數的可選列表。
method	ties 參數的備用名稱。
nocenter	X 矩陣中其值嚴格位於該集合內的列不會居中。請記住，因子變量變成一組 0/1 列。
...	其他參數將傳遞給coxph.control

細節

比例風險模型通常以每個人的單個生存時間值來表示，並可能進行審查。安徒生和吉爾將同一問題重新表述為計數過程；隨著時間的推移，我們觀察某個主題的事件，就像觀看蓋革計數器一樣。受試者的數據以多行或"observations" 的形式呈現，每一行都適用於一個觀察間隔(開始、停止]。

該例程在內部縮放和居中數據以避免指數函數的參數溢出。這些操作不會改變結果，但會帶來更高的數值穩定性。 X 矩陣中值位於 nocenter 列表內的任何列都不會居中。默認的實際結果是不將與因子相對應的虛擬變量居中。但是，偏移量的參數不會縮放，因為在某些情況下會故意使用大偏移值。然而，一般來說，用戶不應避免使用非常大的偏移數值，因為可能會損失估計的精度。

值

表示擬合的類 coxph 的對象。有關詳細信息，請參閱coxph.object 和coxphms.object。

副作用

根據調用情況， predict 、 residuals 和 survfit 例程可能需要重建 coxph 創建的 x 矩陣。此操作可能會失敗，如下麵的示例所示，其中預測函數無法找到 tform 。

  tfun <- function(tform) coxph(tform, data=lung)
  fit <- tfun(Surv(time, status) ~ age)
  predict(fit)

在這種情況下，請將 model=TRUE 選項添加到 coxph 調用中，以避免重建的需要，但代價是需要更大的 fit 對象。

箱重

案例權重被視為複製權重，即案例權重為 2 相當於該受試者的觀察有 2 個副本。當計算機較小時，將相似的主題分組在一起是節省內存的常見技巧。例如，將所有權重設置為 2 將給出相同的係數估計，但方差減半。當使用關係的 Efron 近似(默認)時，數據的複製不會給出與權重選項完全相同的係數，在這種情況下，加權擬合可以說是正確的。

當模型包含 cluster 項或 robust=TRUE 選項時，計算的方差將任何權重視為采樣權重；在這種情況下，將所有權重設置為 2 將得到與權重 1 相同的方差。

特別條款

模型方程中可以使用三個特殊術語。 strata 術語標識分層 Cox 模型；每個層都適合單獨的基線危險函數。 cluster 項用於計算模型的穩健方差。術語+ cluster(id)(其中id 的每個值都是唯一的)相當於指定robust=TRUE 參數。如果 id 變量不唯一，則假定它標識相關觀測值的聚類。穩健的估計源於許多不同的論點，因此有許多標簽。它被稱為 Huber 三明治估計量、White 估計量(線性模型/計量經濟學)、Horvitz-Thompson 估計量(調查抽樣)、工作獨立方差(廣義估計方程)、無窮小折刀和 Wei、Lin、Weissfeld (WLW)估計。

time-transform 術語允許變量隨時間動態變化。在這種情況下，tt 參數將是一個函數或一組函數(如果模型中存在多個 tt() 項)，提供適當的變換。請參閱下麵的示例。

最近出現的一個用戶錯誤是盲目地遵循一些編碼指南的建議，並將 survival:: 添加到所有內容之前，包括特殊術語，例如 survival::coxph(survival:Surv(time, status) ~ age + survival::cluster(inst), data=lung) 首先，這是不必要的：coxph 調用中的參數將是在生存命名空間內進行評估，因此其他包的 Surv 或 cluster 函數不會被注意到。 (然而，coxph 調用本身的完整資格可能具有保護作用。)其次，更重要的是，上麵的調用不會給出正確的答案。特價商品通過其名稱進行識別，survival::cluster 與 cluster 不同；上述模型將 inst 視為普通變量。在生存模型或 glm 模型中使用 stats::offset 作為術語也會出現類似的問題。

收斂

在某些數據情況下，係數的實際 MLE 估計值為無窮大，例如，其中一組沒有事件的二分變量。當這種情況發生時，相關係數會以穩定的速度增長，並且擬合例程中將存在競爭條件：或者對數似然收斂，或者信息矩陣變得有效奇異，或者 exp 的參數對於計算機硬件來說太大，或者最大值交互次數超出。 (最常見的是，數字 1 是第一個發生的。)例程嘗試檢測何時發生這種情況，但並不總是成功。對於用戶來說，主要結果是 Wald 統計量 = 係數/se(係數) 在這種情況下無效，應被忽略；然而，似然比和分數檢驗仍然有效。

領帶

處理綁定事件時間有三種可能的選擇。布雷斯洛近似是最容易編程的，因此成為幾乎所有計算機例程編碼的第一個選項。當添加其他選項以“保持向後兼容性”時，它最終成為默認選項。如果存在大量聯係，Efron 選項會更準確，並且它是此處的默認選項。在實踐中，關係的數量通常很少，在這種情況下，所有方法在統計上都無法區分。

使用“精確偏似然”方法，Cox 偏似然相當於匹配邏輯回歸的偏似然。 (clogit 函數使用 coxph 代碼進行擬合。)當時間尺度是離散的並且隻有幾個唯一值時，它在技術上是合適的，並且某些包將此稱為 "discrete" 選項。由於 Prentice，還有一個“精確邊際可能性”，但此處未實現。

準確的部分似然的計算在數值上是非常密集的。例如，假設第 7 天有 180 名受試者處於危險之中，其中 15 名受試者發生了事件；那麽代碼需要計算所有 180-choose-15 > 10^43 個大小為 15 的不同可能子集的總和。此任務有一個有效的遞歸算法，但即使這樣，計算也可能會長得令人難以忍受。對於(開始，停止)數據，情況會更糟，因為遞歸需要為每個唯一的開始時間重新開始。

多狀態模型是一個更困難的情況。首先，對 Efron 論證進行適當的擴展要困難得多，而且作者尚未完全相信所得到的算法是站得住腳的。其次，Efron 案例的當前代碼並不能一致地計算擴展邏輯(並且擴展需要對代碼進行重大更改)。由於這種複雜性，對於多狀態情況，默認值為ties='breslow'。如果選擇ties='efron'，則當前代碼實際上僅適用於相同類型的綁定轉換。

另一個問題是由於浮點不精確而導致的人為聯係。請參閱有關此主題的小插圖以獲取完整說明或 coxph.control 中的 timefix 選項。用戶可能需要添加timefix=FALSE來模擬數據集。

懲罰回歸

coxph 可以使用任意用戶定義的懲罰來最大化懲罰部分似然。提供的懲罰函數包括嶺回歸 (ridge)、平滑樣條曲線 (pspline) 和脆弱模型 (frailty)。

例子

# Create the simplest test data set 
test1 <- list(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              x=c(0,2,1,1,1,0,0), 
              sex=c(0,0,0,0,1,1,1)) 
# Fit a stratified model 
coxph(Surv(time, status) ~ x + strata(sex), test1) 
# Create a simple data set for a time-dependent model 
test2 <- list(start=c(1,2,5,2,1,7,3,4,8,8), 
              stop=c(2,3,6,7,8,9,9,9,14,17), 
              event=c(1,1,1,1,1,1,1,0,0,0), 
              x=c(1,0,0,1,0,1,1,1,0,0)) 
summary(coxph(Surv(start, stop, event) ~ x, test2)) 

#
# Create a simple data set for a time-dependent model
#
test2 <- list(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
                stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
                event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
                x    =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0) )


summary( coxph( Surv(start, stop, event) ~ x, test2))

# Fit a stratified model, clustered on patients 

bladder1 <- bladder[bladder$enum < 5, ] 
coxph(Surv(stop, event) ~ (rx + size + number) * strata(enum),
      cluster = id, bladder1)

# Fit a time transform model using current age
coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,
     tt=function(x,t,...) pspline(x + t/365.25))

參考

Andersen, P. and Gill, R. (1982). Cox's regression model for counting processes, a large sample study. Annals of Statistics 10, 1100-1120.

Therneau, T., Grambsch, P., Modeling Survival Data: Extending the Cox Model. Springer-Verlag, 2000.

也可以看看

coxph.object , coxphms.object , coxph.control , cluster , strata , Surv , survfit , pspline。