使用时间相关系数和样条从 coxph 对象绘制估计的 HR
Plotting estimated HR from coxph object with time-dependent coefficient and splines
我想在具有基于样条项的时间相关系数的 coxph
模型的情况下将估计的风险比绘制为时间的函数。我使用函数 tt
创建了时间相关系数,类似于直接来自 ?coxph
:
的这个例子
# Fit a time transform model using current age
cox = coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,
tt=function(x,t,...) pspline(x + t/365.25))
调用 survfit(cox)
导致 survfit
无法理解具有 tt
项 (as described in 2011 by Terry Therneau) 的模型。
您可以使用 cox$linear.predictors
提取线性预测变量,但我需要以某种方式提取年龄和更简单的时间来提取每个变量。因为 tt
在事件时间拆分数据集,所以我不能只将输入数据帧的列与 coxph
输出匹配。此外,我真的很想绘制估计函数本身,不只是对观察到的数据点的预测。
这里有涉及样条,但不涉及tt
.
编辑(7/7)
我仍然坚持这一点。我一直在深入研究这个对象:
spline.obj = pspline(lung$age)
str(spline.obj)
# something that looks very useful, but I am not sure what it is
# cbase appears to be the cardinal knots
attr(spline.obj, "printfun")
function (coef, var, var2, df, history, cbase = c(43.3, 47.6,
51.9, 56.2, 60.5, 64.8, 69.1, 73.4, 77.7, 82, 86.3, 90.6))
{
test1 <- coxph.wtest(var, coef)$test
xmat <- cbind(1, cbase)
xsig <- coxph.wtest(var, xmat)$solve
cmat <- coxph.wtest(t(xmat) %*% xsig, t(xsig))$solve[2, ]
linear <- sum(cmat * coef)
lvar1 <- c(cmat %*% var %*% cmat)
lvar2 <- c(cmat %*% var2 %*% cmat)
test2 <- linear^2/lvar1
cmat <- rbind(c(linear, sqrt(lvar1), sqrt(lvar2), test2,
1, 1 - pchisq(test2, 1)), c(NA, NA, NA, test1 - test2,
df - 1, 1 - pchisq(test1 - test2, max(0.5, df - 1))))
dimnames(cmat) <- list(c("linear", "nonlin"), NULL)
nn <- nrow(history$thetas)
if (length(nn))
theta <- history$thetas[nn, 1]
else theta <- history$theta
list(coef = cmat, history = paste("Theta=", format(theta)))
}
所以,我有了结,但我仍然不确定如何将 coxph
系数与结结合起来以实际绘制函数。非常感谢任何线索。
我认为您需要的可以通过使用 pspline
生成输入矩阵并将其乘以 coxph
输出的相关系数来生成。要获得 HR,您需要取指数。
即
output <- data.frame(Age = seq(min(lung$age) + min(lung$time) / 365.25,
max(lung$age + lung$time / 365.25),
0.01))
output$HR <- exp(pspline(output$Age) %*% cox$coefficients[-1] -
sum(cox$means[-1] * cox$coefficients[-1]))
library("ggplot2")
ggplot(output, aes(x = Age, y = HR)) + geom_line()
注意这里的年龄是感兴趣时间的年龄(即基线年龄和自进入研究以来经过的时间之和)。它必须使用指定的范围来匹配原始模型中的参数。也可以使用 x = TRUE
的 x
输出来计算,如下所示:
cox <- coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,
tt=function(x,t,...) pspline(x + t/365.25), x = TRUE)
index <- as.numeric(unlist(lapply(strsplit(rownames(cox$x), "\."), "[", 1)))
ages <- lung$age[index]
output2 <- data.frame(Age = ages + cox$y[, 1] / 365.25,
HR = exp(cox$x[, -1] %*% cox$coefficients[-1] -
sum(cox$means[-1] * cox$coefficients[-1])))
我想在具有基于样条项的时间相关系数的 coxph
模型的情况下将估计的风险比绘制为时间的函数。我使用函数 tt
创建了时间相关系数,类似于直接来自 ?coxph
:
# Fit a time transform model using current age
cox = coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,
tt=function(x,t,...) pspline(x + t/365.25))
调用 survfit(cox)
导致 survfit
无法理解具有 tt
项 (as described in 2011 by Terry Therneau) 的模型。
您可以使用 cox$linear.predictors
提取线性预测变量,但我需要以某种方式提取年龄和更简单的时间来提取每个变量。因为 tt
在事件时间拆分数据集,所以我不能只将输入数据帧的列与 coxph
输出匹配。此外,我真的很想绘制估计函数本身,不只是对观察到的数据点的预测。
这里有tt
.
编辑(7/7)
我仍然坚持这一点。我一直在深入研究这个对象:
spline.obj = pspline(lung$age)
str(spline.obj)
# something that looks very useful, but I am not sure what it is
# cbase appears to be the cardinal knots
attr(spline.obj, "printfun")
function (coef, var, var2, df, history, cbase = c(43.3, 47.6,
51.9, 56.2, 60.5, 64.8, 69.1, 73.4, 77.7, 82, 86.3, 90.6))
{
test1 <- coxph.wtest(var, coef)$test
xmat <- cbind(1, cbase)
xsig <- coxph.wtest(var, xmat)$solve
cmat <- coxph.wtest(t(xmat) %*% xsig, t(xsig))$solve[2, ]
linear <- sum(cmat * coef)
lvar1 <- c(cmat %*% var %*% cmat)
lvar2 <- c(cmat %*% var2 %*% cmat)
test2 <- linear^2/lvar1
cmat <- rbind(c(linear, sqrt(lvar1), sqrt(lvar2), test2,
1, 1 - pchisq(test2, 1)), c(NA, NA, NA, test1 - test2,
df - 1, 1 - pchisq(test1 - test2, max(0.5, df - 1))))
dimnames(cmat) <- list(c("linear", "nonlin"), NULL)
nn <- nrow(history$thetas)
if (length(nn))
theta <- history$thetas[nn, 1]
else theta <- history$theta
list(coef = cmat, history = paste("Theta=", format(theta)))
}
所以,我有了结,但我仍然不确定如何将 coxph
系数与结结合起来以实际绘制函数。非常感谢任何线索。
我认为您需要的可以通过使用 pspline
生成输入矩阵并将其乘以 coxph
输出的相关系数来生成。要获得 HR,您需要取指数。
即
output <- data.frame(Age = seq(min(lung$age) + min(lung$time) / 365.25,
max(lung$age + lung$time / 365.25),
0.01))
output$HR <- exp(pspline(output$Age) %*% cox$coefficients[-1] -
sum(cox$means[-1] * cox$coefficients[-1]))
library("ggplot2")
ggplot(output, aes(x = Age, y = HR)) + geom_line()
注意这里的年龄是感兴趣时间的年龄(即基线年龄和自进入研究以来经过的时间之和)。它必须使用指定的范围来匹配原始模型中的参数。也可以使用 x = TRUE
的 x
输出来计算,如下所示:
cox <- coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,
tt=function(x,t,...) pspline(x + t/365.25), x = TRUE)
index <- as.numeric(unlist(lapply(strsplit(rownames(cox$x), "\."), "[", 1)))
ages <- lung$age[index]
output2 <- data.frame(Age = ages + cox$y[, 1] / 365.25,
HR = exp(cox$x[, -1] %*% cox$coefficients[-1] -
sum(cox$means[-1] * cox$coefficients[-1])))