R中的时间相关系数 - 怎么做?

机器算法验证 r 回归 生存 cox模型
2022-01-21 12:11:39

更新:很抱歉再次更新,但我找到了一些可能的分数多项式解决方案和竞争风险包,我需要一些帮助。


问题

我找不到在 R 中进行时间相关系数分析的简单方法。我希望能够将变量系数转换为时间相关系数(非变量),然后绘制随时间变化的变化:

βmy_variable=β0+β1t+β2t2...

可能的解决方案

1)拆分数据集

我看过这个例子(Se part 2 of the lab session)但是创建一个单独的数据集似乎很复杂,计算成本很高而且不是很直观......

2) 降阶模型 - coxvc 包

coxvc提供了一种处理问题的优雅方法——这里有一个手册问题是作者不再开发包(最新版本是自 2007 年 5 月 23 日以来),经过一些电子邮件对话后,我已经让包工作了,但在我的数据集上运行一次需要 5 小时(140 000条目)并在期末给出极端估计。你可以在这里找到一个稍微更新的——我大部分只是更新了绘图功能。

这可能只是一个调整的问题,但由于该软件不容易提供置信区间并且该过程非常耗时,我现在正在寻找其他解决方案。

3)timereg包

令人印象深刻的timereg 包也解决了这个问题,但我不确定如何使用它,它并没有给我一个流畅的情节。

4)分数多项式时间(FPT)模型

我发现 Anika Buchholz 的优秀论文“评估时间变化的治疗和预后因素的长期影响”在涵盖不同模型方面做得非常出色。她得出结论,Sauerbrei 等人提出的 FPT似乎最适合时间相关系数:

FPT 非常擅长检测时变效应,而 Reduced Rank 方法导致模型过于复杂,因为它不包括时变效应的选择。

这项研究似乎非常完整,但对我来说有点遥不可及。我也有点想知道,因为她恰好和 Sauerbrei 一起工作。不过这似乎很合理,我想可以使用mfp 包进行分析,但我不确定如何。

5) cmprsk 包

我一直在考虑进行竞争风险分析,但计算非常耗时,所以我改用常规的 cox 回归。crr具有时间相关协变量的选项

....
cov2        matrix of covariates that will be multiplied 
            by functions of time; if used, often these 
            covariates would also appear in cov1 to give 
            a prop hazards effect plus a time interaction
....

有二次示例,但我不太了解时间实际出现的位置,我不确定如何显示它。我还查看了 test.R 文件,但那里的示例基本相同......

我的示例代码

这是我用来测试不同可能性的示例

library("survival")
library("timereg")
data(sTRACE)

# Basic cox regression    
surv <- with(sTRACE, Surv(time/365,status==9))
fit1 <- coxph(surv~age+sex+diabetes+chf+vf, data=sTRACE)
check <- cox.zph(fit1)
print(check)
plot(check, resid=F)
# vf seems to be the most time varying

######################################
# Do the analysis with the code from #
# the example that I've found        #
######################################

# Split the dataset according to the splitSurv() from prof. Wesley O. Johnson
# http://anson.ucdavis.edu/~johnson/st222/lab8/splitSurv.ssc
new_split_dataset = splitSuv(sTRACE$time/365, sTRACE$status==9, sTRACE[, grep("(age|sex|diabetes|chf|vf)", names(sTRACE))])

surv2 <- with(new_split_dataset, Surv(start, stop, event))
fit2 <- coxph(surv2~age+sex+diabetes+chf+I(pspline(stop)*vf), data=new_split_dataset)
print(fit2)

######################################
# Do the analysis by just straifying #
######################################
fit3 <- coxph(surv~age+sex+diabetes+chf+strata(vf), data=sTRACE)
print(fit3)

# High computational cost!
# The price for 259 events
sum((sTRACE$status==9)*1)
# ~240 times larger dataset!
NROW(new_split_dataset)/NROW(sTRACE)

########################################
# Do the analysis with the coxvc and   #
# the timecox from the timereg library #
########################################
Ft_1 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=3))
fit_coxvc1 <- coxvc(surv~vf+sex, Ft_1, rank=2, data=sTRACE)

fit_coxvc2 <- coxvc(surv~vf+sex, Ft_1, rank=1, data=sTRACE)

Ft_3 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=5))
fit_coxvc3 <- coxvc(surv~vf+sex, Ft_3, rank=2, data=sTRACE)

layout(matrix(1:3, ncol=1))
my_plotcoxvc <- function(fit, fun="effects"){
    plotcoxvc(fit,fun=fun,xlab='time in years', ylim=c(-1,1), legend_x=.010)
    abline(0,0, lty=2, col=rgb(.5,.5,.5,.5))
    title(paste("B-spline =", NCOL(fit$Ftime)-1, "df and rank =", fit$rank))
}
my_plotcoxvc(fit_coxvc1)
my_plotcoxvc(fit_coxvc2)
my_plotcoxvc(fit_coxvc3)

# Next group
my_plotcoxvc(fit_coxvc1)

fit_timecox1<-timecox(surv~sex + vf, data=sTRACE)
plot(fit_timecox1, xlab="time in years", specific.comps=c(2,3))

代码生成以下图表: coxvc 的不同设置以及coxvc 和 timecox图的比较。我想结果还可以,但我认为我无法解释 timecox 图 - 它似乎很复杂......

我的(当前)问题

  • 如何在 R 中进行 FPT 分析?
  • 如何在 cmprsk 中使用时间协变量?
  • 如何绘制结果(最好使用置信区间)?
3个回答

@mpiktas 提供了一个可行的模型,但是在 time=t 中需要用于二次方的术语是I(t^2))之所以如此,是因为在 R 中,“^”的公式解释创建了交互作用并且不执行幂运算,因此“t”与“t”的交互作用只是“t”。(这不应该用 [r] 标签迁移到 SO 吗?)

对于这个过程的替代方案,在我看来这对于推理目的有点可疑,并且可能符合您对使用 Harrell 的 rms/Hmisc 包中的支持功能的兴趣,请参阅 Harrell 的“回归建模策略”。他提到(但只是顺便提及,尽管他确实引用了他自己的一些论文)构造样条拟合到时间尺度来模拟浴缸形危险。他关于参数生存模型的章节描述了各种绘图技术,用于检查比例风险假设和检查估计的对数风险效应在时间尺度上的线性度。

编辑:另一个选项是使用coxph'tt' 参数,描述为“时间转换函数的可选列表”。

我已经改变了这个答案,因为@DWin 或@Zach 的答案都没有完全回答如何对时变系数建模。我最近写了一篇关于这个的帖子。这是它的要点。

Cox 回归模型中的核心概念风险函数h (它被定义为:h(t)

h(t)=f(t)S(t)

其中是在任何给定时间发生事件的风险,而是在该票价中幸存的概率。因此,该数字是理论范围从的分数。f(t)S(t)0

危险函数的一个有趣特征是,我们可以包括在之外的其他时间点的观察,例如,如果“彼得”在英国接受髋关节置换术,1 年后到达瑞典,他已经活了 1 年我们选择包括他。请注意,在此之前已经死亡的任何患者都不会引起我们的注意,因此在查看 1 年前的危害时一年后,我们可能会包括彼得。time0S(t)

当允许受试者在其他时间点进入时,我们必须将Survfrom更改Surv(time, status)Surv(start_time, end_time, status)虽然end_time与结果密切相关,但start_time现在可以作为交互项使用(正如原始建议中所暗示的那样)。在常规设置中,start_time除了一些稍后出现的主题之外,它是 0,但是我们将每个观察分成几个时期,我们突然有很多非零的开始时间。唯一的区别是每个观察都发生多次,除了最后一个观察之外,所有观察都可以选择非审查结果。

实践中的时间分割

我刚刚在 CRAN 上发布了Greg包,它使这个时间分割变得容易。首先,我们从一些理论观察开始:

library(Greg)
test_data <- data.frame(
  id = 1:4,
  time = c(4, 3.5, 1, 5),
  event = c("censored", "dead", "alive", "dead"),
  age = c(62.2, 55.3, 73.7, 46.3),
  date = as.Date(
    c("2003-01-01", 
      "2010-04-01", 
      "2013-09-20",
      "2002-02-23"))
)

我们可以用 * 表示事件以图形方式显示:

在此处输入图像描述

如果我们应用timeSplitter如下:

library(dplyr)
split_data <- 
  test_data %>% 
  select(id, event, time, age, date) %>% 
  timeSplitter(by = 2, # The time that we want to split by
               event_var = "event",
               time_var = "time",
               event_start_status = "alive",
               time_related_vars = c("age", "date"))

我们得到以下信息:

在此处输入图像描述

如您所见,每个对象已被拆分为多个事件,其中最后一个时间跨度包含实际事件状态。这使我们现在可以构建具有简单:交互项的模型(不要使用 ,*因为这会扩展time + var + time:var,我们对时间本身不感兴趣)。不需要使用该I()函数,尽管如果您想检查非线性与时间,我经常创建一个单独的时间交互变量,我添加一个样条曲线,然后使用rms::contrast. 无论如何,您的回归调用现在应该如下所示:

coxp(Surv(start_time, end_time, event) ~ var1 + var2 + var2:time, 
     data = time_split_data)

使用生存包的tt功能

还有一种方法可以使用该tt函数直接在生存包中对时间相关系数进行建模。Therneau 教授在他的小插图中对该主题进行了全面介绍。不幸的是,在大型数据集中,由于内存限制,这很难做到。似乎该tt函数将时间分成非常精细的片段,在此过程中生成了一个巨大的矩阵。

您可以使用PerformanceAnalytics中的apply.rolling函数通过滚动窗口运行线性回归,这将允许您的系数随时间变化。

例如:

library(PerformanceAnalytics)
library(quantmod)
getSymbols(c('AAPL','SPY'), from='01-01-1900')
chart.RollingRegression(Cl(AAPL),Cl(SPY), width=252, attribute='Beta')
#Note: Alpha=y-intercept, Beta=regression coeffient

这也适用于其他功能。