I'm performing a cross validation on a competing risks proportional hazards model. With help from the mstate pacakge, I've prepared my data and am fitting it with survival::coxph. I get a fitted Cox model object for my training data, but I want to evaluate the partial likelihood of my trained coefficients with my test data.
If I need to, I'll write the partial likelihood function myself, but I'd rather not (though it would probably be good for me). The survival package calculates in this C code, but the likelihood calculation is embedded in the fitting function. Maybe there's a way to fix parameters, or some other tools to easily get at the partial likelihood?
Minimum Working Exmaple
# 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")
Data looks like this:
> 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
Now, the ccr5 variable could be modeled as transition-specific, or as a having equal proportional effect for all transitions. The models are:
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)
Now I would like to use the test data to evaluate the variable selection
on whether or not ccr5 should be transition-specific or not.
I have a large data set and many variables--mostly but not all categorical--that could go either way. The evaluation is where I'm stuck.
# 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")
I'm hoping to calculate the partial likelihood for each of the models on the test data with the training coefficient estimates. Maybe I should move the question to Cross Validated and ask if the sum of the linear predictors (or the sum of the linear predictors excluding censored cases) is close enough to an equivalent measure.
