1

我正在阅读Singer 和 Willett的Applied Longitudinal Data-Analysis的第 15 章，关于扩展 Cox 回归模型，但这里的 UCLA 网站没有本章的示例 R 代码。我正在尝试重新创建关于时变协变量的部分，并坚持如何从提供的人员级别数据框创建计数过程数据集。我查看了survival包小插图，但在创建自己的数据集时遇到了麻烦。这是以 Singer 和 Willett 示例为模型的玩具数据，该示例预测了 17 至 40 岁男性首次使用可卡因的时间。id是 ID，ageInit是第一次使用可卡因或审查的年龄，used是事件日志，表示第一次使用可卡因 ( 1) 或审查 (0）。有一个时间不变预测因子earlyMJ，表明 17 岁之前吸食大麻。还有三个时变协变量需要告知数据集的创建：timeMJ，这是 17 岁之后第一次使用大麻sellMJ的年龄，这是销售大麻的年龄，以及odFirst，这是第一次使用其他药物的年龄。NA这些预测变量中的 s 表示参与者在任何时候都没有执行相关操作。

set.seed(1356)
df <- data.frame(id = 1:6,
                 ageInit = c(25,34,40,29,27,40),
                 used = c(0,1,1,0,1,1), 
                 earlyMJ = c(0,0,1,0,1,1),
                 timeMJ = c(18,27,22,21,22,19),
                 sellMJ = c(NA,NA,25,NA,35,NA),
                 odFirst = c(19,22,35,NA,22,34))

遵循 Therneau 等人在survival小插图中的过程，我们创建了第二个数据集，在这种情况下，将结果变量 和随时间变化的预测变量重新表达ageInit为从研究开始的年数（即 17 岁）。 

tdata <- with(df, data.frame (id = id,
                              usedTime = ageInit - 17,
                              timeToFirstMJ = timeMJ - 17,
                              timeToSellMJ = sellMJ - 17,
                              timeToFirstOD = odFirst - 17,
                              usedCocaine = used))

我们将这些新数据与原始数据集合并，通过event()调用创建两个新列，表示每个协变量的时间间隔。我们还通过tdc()调用为每个参与者每次经历一个协变量事件时创建一个新行。

sdata <- tmerge(df, tdata, 
                id=id, 
                firstUse = event(futime, usedCocaine), 
                t1MJ   =  tdc(timeToFirstMJ),
                t1SMJ  =  tdc(timeToSellMJ),
                t1OD  =  tdc(timeToFirstOD),
                options= list(idname="subject"))
attr(sdata, "tcount")

问题是当我运行 Cox 回归时，同时使用时不变和一些时变协变量，我无法使模型工作。 

coxph(Surv(tstart, tstop, firstUse) ~ earlyMJ + t1SMJ + t1OD, data= sdata, ties="breslow")

并获得无意义的系数和警告信息

In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Loglik converged before variable  1,3 ; beta may be infinite.

此外，我什至不知道这是否是一个正确的计数过程数据框，因为在 Singer 和 Willett 中，他们建议每个参与者在每次数据集中的任何人经历事件时都需要有一行。

在这些问题上的任何指导将不胜感激。

4

1 回答 1

4

一个好的开始是在包中的 Cox 模型小插图中使用时间相关协变量和时间相关系数survival。

问题是当我运行 Cox 回归时，同时使用时不变和一些时变协变量，我无法使模型工作。

乍一看，您所做的似乎是正确的，但可能只是您只有 6 个人和三个参数要估计的事实？

得到无意义的系数和警告信息

这个警告很有道理。您的两个参数估计值的标准误差大于 1000。见下文。

此外，甚至不知道这是否是一个正确的计数过程数据框，因为在 Singer 和 Willett 中，他们建议每个参与者在每次数据集中的任何人经历事件时都需要有一行。

这是在内部处理的coxph。

这是我的代码来重现您的结果

#####
# setup data
df <- data.frame(id = 1:6,
                 ageInit = c(25,34,40,29,27,40),
                 used =    c( 0, 1, 1, 0, 1, 1), 
                 earlyMJ = c( 0, 0, 1, 0, 1, 1),
                 timeMJ  = c(18,27,22,21,22,19),
                 sellMJ =  c(NA,NA,25,NA,35,NA),
                 odFirst = c(19,22,35,NA,22,34))
shift_cols <- c("ageInit", "timeMJ", "sellMJ", "odFirst")
df[shift_cols] <- lapply(df[shift_cols], "-", 17)

library(survival)
est_df <- df[, c("id", "ageInit", "earlyMJ", "used")]
est_df <- tmerge(
  est_df, est_df, id = id, start_using = event(ageInit, used))
est_df <- tmerge(
  est_df, df, id = id, 
  t1MJ =  tdc(timeMJ), t1SMJ = tdc(sellMJ), t1OD = tdc(odFirst))

#####
# fit model
fit <- coxph(
  Surv(tstart, tstop, start_using) ~ earlyMJ + t1SMJ + t1OD, data = est_df)
#R> Warning message:
#R> In fitter(X, Y, strats, offset, init, control, weights = weights,  :
#R>   Loglik converged before variable  1,3 ; beta may be infinite. 

summary(fit)
#R> Call:
#R> coxph(formula = Surv(tstart, tstop, start_using) ~ earlyMJ + 
#R>     t1SMJ + t1OD, data = est_df)
#R> 
#R>   n= 17, number of events= 4 
#R> 
#R>              coef exp(coef)  se(coef)     z Pr(>|z|)
#R> earlyMJ 2.047e+01 7.792e+08 2.791e+04 0.001    0.999
#R> t1SMJ   4.744e-16 1.000e+00 1.414e+00 0.000    1.000
#R> t1OD    4.157e+01 1.136e+18 3.883e+04 0.001    0.999
#R> 
#R>         exp(coef) exp(-coef) lower .95 upper .95
#R> earlyMJ 7.792e+08  1.283e-09   0.00000       Inf
#R> t1SMJ   1.000e+00  1.000e+00   0.06255     15.99
#R> t1OD    1.136e+18  8.806e-19   0.00000       Inf
#R> 
#R> Concordance= 1  (se = 0.258 )
#R> Rsquare= 0.273   (max possible= 0.33 )
#R> Likelihood ratio test= 5.42  on 3 df,   p=0.1437
#R> Wald test            = 0  on 3 df,   p=1
#R> Score (logrank) test = 4.14  on 3 df,   p=0.2463
于 2018-03-11T21:44:30.817 回答