6

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术语的模型的错误（如 Terry Therneau 在 2011 年所述）。

您可以使用 提取线性预测器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系数与结结合起来以实际绘制函数。任何线索都非常感谢。

4

1 回答 1

5

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

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()

人力资源与年龄的关系图

注意这里的年龄是感兴趣的时间的年龄（即基线年龄和自研究开始以来经过的时间的总和）。它必须使用指定的范围来匹配原始模型中的参数。它也可以使用 using 的x输出来计算x = TRUE，如下所示：

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])))
于 2015-07-09T11:26:04.523 回答