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。