I’m running into an error when trying to obtain survival curves from a Cox PH model fitted using timeline-type data with a Surv2 formula.

Error in aeqSurv(Y) : argument is not a Surv object

Here is a barebones example, using the mgus data and the exact code from the vignette.

# Preparing the dataframe in timeline format, per the vignette
ctime <- with(mgus2, ifelse(pstat==1, ptime, futime))
cstat <- with(mgus2, ifelse(pstat==1, 1, 2*death))
cstat <- factor(cstat, 0:2, c("censor", "PCM", "death"))
tdata <- data.frame(id=mgus2$id, days=ctime, cstat=cstat)
mdata2 <- data.frame(mgus2[,1:7], days=0)
mdata2 <- merge(mdata2, tdata, all=TRUE)

# Fitting the cox model
cph <-  coxph(Surv2(days, cstat) ~ age + sex, id=id, mdata2)

# Attempting to obtain survival curves for two new records
ndata <- expand.grid(sex=c("F", "M"), age=c(60, 80))
survfit(cph, newdata=ndata) 
> Error in aeqSurv(Y) : argument is not a Surv object

I have traced the error inside survfit.coxph to the variable mf not taking the value it’s supposed to take (survfit.coxph.R:83):

mf <- stats::model.frame(object)

When called with a coxph object that uses a Surv-type formula, mf is as expected. When the formula is a Surv2, however, it appears to return an empty frame, which then dominoes into the error I’m experiencing.

Am I making a mistake on my side when calling survfit.coxph? Or does it currently not support Surv2, making it an internal package problem?