使用与时间有关的系数和样条曲线绘制来自Coxph对象的估计HR [英] Plotting estimated HR from coxph object with time-dependent coefficient and splines

问题描述

在 coxph 模型具有基于时间的系数（基于样条项）的模型的情况下，我想将估计的危险比作为时间的函数进行绘制。我使用函数 tt 创建了时间相关系数，类似于本示例直接来自？coxph ：

I want to plot the estimated hazard ratio as a function of time in the case of a coxph model with a time-dependent coefficient that is based on a spline term. I created the time-dependent coefficient using function tt, analogous to this example that comes straight from ?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 术语的模型（ ）。

Calling survfit(cox) results in an error that survfit does not understand models with a tt term (as described in 2011 by Terry Therneau).

您可以使用 cox $ linear.predictors 提取线性预测变量，但我需要以某种方式提取年龄和不那么琐碎，每个时代都有。因为 tt 在事件发生时拆分数据集，所以我不能只将输入数据框的列与 coxph 输出。另外，我真的很想绘制估计函数本身，而不只是绘制观察到的数据点的预测。

You can extract the linear predictor using cox$linear.predictors, but I would need to somehow extract ages and less trivially, times to go with each. Because tt splits the dataset on event times, I can't just match up the columns of the input dataframe with the coxph output. Additionally, I really would like to plot the estimated function itself, not just the predictions for the observed data points.

有此处涉及样条线的相关问题，但不涉及 tt 。

我仍然坚持这样做。我一直在深入研究此对象：

I'm still stuck on this. I've been looking in depth at this object:

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 系数与节点结合起来以实际绘制函数。

So, I have the knots, but I am still not sure how to combine the coxph coefficients with the knots in order to actually plot the function. Any leads much appreciated.

推荐答案

我认为可以通过使用生成输入矩阵来生成所需的内容。 pspline 并将其乘以 coxph 输出的相关系数。要获得心率，则需要取指数。

I think what you need can be generated by generating an input matrix using pspline and matrix-multiplying this by the relevant coefficients from the coxph output. To get the HR, you then need to take the exponent.

ie

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 输出，通过使用 x = TRUE 计算得出，如下所示：

Note the age here is the age at the time of interest (i.e. the sum of the baseline age and the elapsed time since study entry). It has to use the range specified to match with the parameters in the original model. It could also be calculated using the x output from using x = TRUE as shown:

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对象的估计HR的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！