【发布时间】：2014-11-03 19:27:00
【问题描述】：
我正在对竞争风险比例风险模型进行交叉验证。在mstate pacakge 的帮助下，我已经准备好我的数据并正在使用survival::coxph 进行拟合。我为我的训练数据获得了一个合适的 Cox 模型对象，但我想用我的测试数据评估我训练的系数的部分可能性。
如果需要，我会自己编写部分似然函数，但我宁愿不写（尽管它可能对我有好处）。生存包计算in this C code，但似然计算嵌入在拟合函数中。也许有一种方法可以修复参数，或者其他一些工具可以轻松获得部分可能性？
最低工作示例
# Adapted from examples in the mstate vignette
# http://cran.r-project.org/web/packages/mstate/vignettes/Tutorial.pdf
# beginning at the bottom of page 28
library(mstate)
library(survival)
# Get data. I add a second explanatory variable (badx) for illustration
# Also divide the data by subject into training and test sets.
data(aidssi)
si <- aidssi # Just a shorter name
si$badx <- sample(c("A", "B"), size = nrow(si), replace = TRUE)
si$fold <- sample(c("train", "test"), size = nrow(si), replace = TRUE, prob = c(0.7, 0.3))
tmat <- trans.comprisk(2, names = c("event-free", "AIDS", "SI"))
si$stat1 <- as.numeric(si$status == 1)
si$stat2 <- as.numeric(si$status == 2)
# Convert the data to a long competing risks format
silong <- msprep(time = c(NA, "time", "time"),
status = c(NA,"stat1", "stat2"),
data = si, keep = c("ccr5", "badx", "fold"), trans = tmat)
silong <- na.omit(silong)
silong <- expand.covs(silong, c("ccr5", "badx"))
train.dat <- subset(silong, fold == "train")
test.dat <- subset(silong, fold == "test")
数据如下所示：
> head(silong)
An object of class 'msdata'
Data:
id from to trans Tstart Tstop time status ccr5 badx fold ccr5WM.1 ccr5WM.2 badxB.1 badxB.2
1 1 1 2 1 0 9.106 9.106 1 WW A train 0 0 0 0
2 1 1 3 2 0 9.106 9.106 0 WW A train 0 0 0 0
3 2 1 2 1 0 11.039 11.039 0 WM B train 1 0 1 0
4 2 1 3 2 0 11.039 11.039 0 WM B train 0 1 0 1
5 3 1 2 1 0 2.234 2.234 1 WW B train 0 0 1 0
6 3 1 3 2 0 2.234 2.234 0 WW B train 0 0 0 1
现在，ccr5 变量可以建模为特定于转换，或对所有转换具有相等的比例效应。模型是：
train.mod.equal <- coxph(Surv(time, status) ~ ccr5 + badx + strata(trans),
data = train.dat)
train.mod.specific <- coxph(Surv(time, status) ~ ccr5WM.1 + ccr5WM.2 + badx + strata(trans),
data = train.dat)
现在我想使用测试数据来评估变量选择
关于ccr5 是否应该是特定于转换的。
我有一个庞大的数据集和许多变量——大部分但不是所有的分类变量——这两种方式都可以。评估是我卡住的地方。
# We can fit the same models to the test data,
# this yields new parameter estimates of course,
# but the model matrices might be useful
test.mod.equal <- coxph(Surv(time, status) ~ ccr5 + badx + strata(trans),
data = test.dat)
test.mod.specific <- coxph(Surv(time, status) ~ ccr5WM.1 + ccr5WM.2 + badx + strata(trans),
data = test.dat)
test.eq.mm <- model.matrix(test.mod.equal)
test.sp.mm <- model.matrix(test.mod.specific)
# We can use these to get the first part of the sum of the partial likelihood:
xbeta.eq <- test.eq.mm[test.dat$status == 1, ] %*% coef(train.mod.equal)
xbeta.sp <- test.sp.mm[test.dat$status == 1, ] %*% coef(train.mod.specific)
# We can also get linear predictors
lp.eq <- predict(train.mod.equal, newdata = test.dat, type = "lp")
lp.sp <- predict(train.mod.specific, newdata = test.dat, type = "lp")
我希望使用训练系数估计值计算测试数据上每个模型的部分似然性。也许我应该将问题移至 Cross Validated 并询问线性预测变量的总和（或不包括审查案例的线性预测变量的总和）是否足够接近等效度量。
【问题讨论】：
-
问题是随着风险集的缩小，部分可能性随着时间变量的变化而变化。如果您解释您实际尝试的内容，最好用一个小示例集，可能会提供具体细节。
-
@BondedDust 我正在研究一个例子。我意识到每次观察对可能性的贡献取决于哪些其他变量仍在风险集中，但总和部分对数似然仍然只是给定数据的模型参数的函数。我有模型参数（来自训练拟合）和数据（来自测试集），我只是缺少计算部分对数似然的函数。
-
因此，您想计算 predict(mdl, type="lp") 与具有审查和结果变量的一组数据的距离（在对数似然量表上）。您能否使用包含使用 beta 估计的偏移量的公式计算“新模型”，然后使用 summary(mdl) 为您完成繁重的工作？您甚至可以使用
predict.coxph计算偏移量。 -
嗯，我得考虑一下。我喜欢新模式的想法。我要参加一个会议，但今晚我会带着 MWE 回来（如果我真的很幸运，我会发布作为答案）。
-
MWE 补充说，“新模型”方法不走运，至少我自己不行。
标签： r cross-validation survival-analysis cox-regression