在 R 中创建具有时变协变量的计数过程数据集
Creating Count Process Data Set With Time-Varying Covariates in R
我正在研究 Singer 和 Willett 的 Applied Longitudinal Data-Analysis 的第 15 章，关于扩展 Cox 回归模型，但是 UCLA 网站 here 没有本章的示例 R 代码。我正在尝试重新创建关于时变协变量的部分，并且一直在研究如何从提供的人员级数据框创建计数过程数据集。我查看了 survival 包小插图，但在创建自己的数据集时遇到了问题。以下是玩具数据，该数据以 Singer 和 Willett 的示例为模型，预测 17 至 40 岁男性首次使用可卡因的时间。 id 是 ID，ageInit 是第一次使用可卡因或审查的年龄，used 是事件日志，表示第一次使用可卡因 (1) 或审查 (0).有一个时间不变的预测因子 earlyMJ，表明 17 岁之前吸食大麻。还有三个时变协变量需要告知数据集的创建：timeMJ，这是 17 岁后第一次吸食大麻的年龄，sellMJ，这是卖大麻的年龄，[=24] =]，这是第一次使用其他药物的年龄。这些预测变量中的 NAs 表示参与者在任何时候都没有执行相关操作。
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 中，他们建议每个参与者每次 anyone 都需要一行数据集经历了事件。
如能提供有关这些问题的任何指导，我们将不胜感激。
一个好的开始是 survival 包中的 Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model 小插图。
The issue is that when I run the Cox regression, using both the time-invariant, and a few of the time-varying covariates, I cannot make the model work.
乍一看你做的似乎是正确的，但可能只是你只有 6 个人和三个参数要估计的事实？
get nonsensical coefficient and the warning message
这个警告很有道理。您的两个参数估计值的标准误差大于 1000。请参见下文。
Furthermore do not even know if this is a proper count process dataframe, because in Singer and Willett they suggest that each participant needs to have a row for every time anyone in the dataset experiences the event.
这是在 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
我正在研究 Singer 和 Willett 的 Applied Longitudinal Data-Analysis 的第 15 章，关于扩展 Cox 回归模型，但是 UCLA 网站 here 没有本章的示例 R 代码。我正在尝试重新创建关于时变协变量的部分，并且一直在研究如何从提供的人员级数据框创建计数过程数据集。我查看了 survival 包小插图，但在创建自己的数据集时遇到了问题。以下是玩具数据，该数据以 Singer 和 Willett 的示例为模型，预测 17 至 40 岁男性首次使用可卡因的时间。 id 是 ID，ageInit 是第一次使用可卡因或审查的年龄，used 是事件日志，表示第一次使用可卡因 (1) 或审查 (0).有一个时间不变的预测因子 earlyMJ，表明 17 岁之前吸食大麻。还有三个时变协变量需要告知数据集的创建：timeMJ，这是 17 岁后第一次吸食大麻的年龄，sellMJ，这是卖大麻的年龄，[=24] =]，这是第一次使用其他药物的年龄。这些预测变量中的 NAs 表示参与者在任何时候都没有执行相关操作。
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 中，他们建议每个参与者每次 anyone 都需要一行数据集经历了事件。
如能提供有关这些问题的任何指导，我们将不胜感激。
一个好的开始是 survival 包中的 Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model 小插图。
The issue is that when I run the Cox regression, using both the time-invariant, and a few of the time-varying covariates, I cannot make the model work.
乍一看你做的似乎是正确的，但可能只是你只有 6 个人和三个参数要估计的事实？
get nonsensical coefficient and the warning message
这个警告很有道理。您的两个参数估计值的标准误差大于 1000。请参见下文。
Furthermore do not even know if this is a proper count process dataframe, because in Singer and Willett they suggest that each participant needs to have a row for every time anyone in the dataset experiences the event.
这是在 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