机器算法验证 - 生存不为0时的Cox回归？ - 吾爱随笔录

生存不为0时的Cox回归？

机器算法验证生存 cox模型

2022-03-21 14:47:45

是否适合使用 cox 回归来研究当时间 = 无穷大时不会导致生存函数变为 0 的过程？Cox 回归对风险比感兴趣，它没有尝试对风险函数本身进行实际估计/建模，所以直觉上我觉得这应该没问题。我错过了什么吗？（如果有帮助，我可以在文献中找到很多这样做的人的例子，但这并不意味着它是正确的......）

举个例子，假设我想研究离婚率……在这种情况下，很多人永远不会离婚，所以婚姻生存率永远不会为 0。

2个回答

有相当多的关于人口异质性下的生存分析的文献，但就像你注意到的那样，我很少看到在应用研究中使用——甚至考虑过这样的模型。我会给出一些简短的直觉，希望其他人可以添加数学上更重的解释，如果那是你正在寻找的。

因此，假设我们的观察结果看起来像这样 - 四分之一的样本“治愈”并且在观察期间没有经历该事件，而其余样本是“易感的”并且慢慢死亡/失败/离婚：

整个蓝色部分将在 $t=50$ ，并且您可以像往常一样在两组中进行额外的审查（图中未显示）。

现在，Cox PH 模型假设生存时间与审查无关。但是，如果您将这样的模型拟合到图中的样本，则审查组将具有~~更高的事件发生时间和~~不同的危险函数，因为它包含整个“治愈”部分。模拟（例如2）表明 HR 估计在这些条件下可能存在很大偏差。
这个问题可以通过排除“治愈”的观察结果来解决，但这些观察结果通常无法与被审查的、但“易感”的个体区分开来。另一种分析选项是仅考虑二元结果事件与无事件，但这再次意味着排除所有被审查的个人。

一个概念上简单的解决方案是将样本视为具有权重的“固化”和“易感”分布的混合物 $p$ 和 $1-p$ , 并拟合生存的混合模型 $S(t) = 1$ 对于“固化”部分和一些减少 $S(t)$ 对于“易感”部分。这些被称为混合固化模型；还有非混合解决方案，但我不确定这些在实践中是否流行。

一些参考文献，按照从最外行到最数学的顺序：

治愈模型作为一种有用的统计工具，用于分析基于广义修正 Weibull 分布的生存混合和非混合治愈
分数模型，并应用于胃癌数据在生存分析

编辑

鉴于下面@JarleTufto 的好评，我应该澄清这个答案。我不建议仅仅在观察期结束时进行审查是不好的，或者说最终审查后的危险 $t$ 在某种程度上是相关的，但问题是由人口的潜在异质性引起的- 即具有不同风险函数的分数 $h(t)$ 审查不同。相应的 Cox 模型假设如下：

h (t | c o v a r i a t e s) = h (t | C_{i} > t, c o v a r i a t e s)

$h(t|covariates) = h(t|C_i > t, covariates)$

（例如，使用多重插补放宽 Cox 比例风险模型中的独立审查假设）

让我们从下面的答案中获取模拟代码，并添加一个“治愈”分数 $h(t)=0$ （还添加了种子并增加了样本量，以便于重现性）：

set.seed(1234)

# simulated survival times from the model
n <- 5000
x <- rnorm(n)
beta <- 0.5
# variation of inversion method
u <- runif(n)
eventtime <- 1/(1 + exp(-beta*x)*log(u)) - 1 
eventtime[eventtime < 0] <- Inf

# simulate independent right censoring points
censoringtime <- runif(n,0,20)

# compute the observed data, that is, the censoring indicator and
# the time of whichever event comes first
delta <- eventtime < censoringtime
time <- pmin(eventtime, censoringtime)

# add "cured" fraction, censored at max observation time:
n2 <- 1000
time2 <- c(time, rep(20, n2))
delta2 <- c(delta, rep(FALSE, n2))
x2 <- c(x, rnorm(n2))

所以观察到的（事件或审查）时间点看起来像：

qplot(time2)

仅使用“敏感”部分，我们得到预期的估计 $\beta=0.5$ ：

# Fit the cox proportional hazards model
library(survival)
model <- coxph(Surv(time,delta) ~ x)
model

    coef exp(coef) se(coef)    z      p
x 0.5067    1.6598   0.0198 25.6 <2e-16

Likelihood ratio test=661  on 1 df, p=0
n= 5000, number of events= 2923

但是，使用完整样本时，HR 被低估了：

model2 <- coxph(Surv(time2,delta2) ~ x2)
model2

    coef exp(coef) se(coef)    z      p
x2 0.4192    1.5207   0.0192 21.8 <2e-16

Likelihood ratio test=478  on 1 df, p=0
n= 6000, number of events= 2923

在您描述的情况下应用 Cox 模型可能完全没问题。Cox 模型不对基线危害做出任何假设 $h_0(t)$ 除了它对所有人都是非负的 $t$ . 因此，基线风险很可能会以足够快的速度变为零以使累积风险 $\int_0^t h_0(u) du$ 去一个有限的值和 $S(t)$ 到一些正值 $t\rightarrow \infty$ .

与参数生存模型不同，除了最后一次审查事件/最后一次观察到的故障之外，不会对危险函数进行任何推断。相反，就像更简单的 Kaplan-Meier 估计 $S(t)$ ，Cox 模型为您提供基线生存函数的非参数估计 $S_0(t)$ 直到最后一次审查事件/最后一次观察到的失败。因此，在这个时间点之后，危险函数的理论行为是没有任何意义的。

例如，假设真正的未知基线危害具有以下形式

h_{0} (t) = \frac{1}{(1 + t)^{2}} .

$h_0(t) = \frac1{(1+t)^2}.$ 具有协变量向量的个体的生存函数

x

$x$ 那么是

\begin{aligned} S (t; x) & = e^{- e^{β x} \int_{0}^{t} h_{0} (u) d u} \\ = e^{- e^{β x} (1 - \frac{1}{t + 1})} . \end{aligned}

$\begin{align} S(t;x) &= e^{-e^{\beta x}\int_0^t h_0(u)du} \\&= e^{-e^{\beta x}(1-\frac1{t+1})}. \end{align}$ 累积基线风险

1 - \frac{1}{t + 1}

$1-\frac1{t+1}$ 然后趋于一个极限值（一个）和生存函数

S (t; x)

$S(t;x)$ 为正概率

e^{- e^{β x}}

$e^{-e^{\beta x}}$ 作为

t \to \infty

$t\rightarrow\infty$ . 请注意，极限值

S (t; x)

$S(t;x)$ 通常必须取决于协变量

x

$x$ ，如果不是这种情况，则将违反比例风险假设。

以下是生存时间 $T$ 从这个模型（其中一些是无限的）和独立的权利审查时间 $C$ 均匀分布在 0 到 20 的区间上进行模拟，并将 Cox 比例风险模型拟合到观察到的数据（删失指标和最小 $T$ 和 $C$ ）。

# simulated survival times from the model
n <- 300
x <- rnorm(n)
beta <- .5
# variation of inversion method
u <- runif(n)
eventtime <- 1/(1 + exp(-beta*x)*log(u)) - 1 
eventtime[eventtime < 0] <- Inf

# simulate independent rigth censoring points
censoringtime <- runif(n,0,20)

# compute the observed data, that is, the censoring indicator and
# the time of whichever event comes first
delta <- eventtime < censoringtime
time <- pmin(eventtime, censoringtime)

# Fit the cox proportional hazards model
library(survival)
model <- coxph(Surv(time,delta) ~ x)

回归系数 $\beta=0.5$ ,

> model
Call:
coxph(formula = Surv(time, delta) ~ x)

    coef exp(coef) se(coef)    z       p
x 0.5463    1.7269   0.0892 6.12 9.2e-10

Likelihood ratio test=39.6  on 1 df, p=3.19e-10
n= 300, number of events= 164     Call:

和基线生存函数 $S_0(t)$ 估计就好了。

# Compare the estimated and true baseline survival functions
plot(survfit(model, newdate=data.frame(x=0)))
curve(exp(-(1-1/(t+1))),xname="t",add=TRUE,col="red")

其它你可能感兴趣的问题

上一篇引入可估计函数的概念有什么意义？下一篇这张数据前后的图表叫什么？