calibrate function from rms package not working for cox model with time-varying covariates?

Question

I am using the calibrate function from the rms package but it returned an error message Error in reliability[, "index.corrected"] : subscript out of bounds. The validate function, however, did work. I am running a cox model with a time-varying covariate and so I am wondering if the calibrate function cannot handle this type of model.

library(survival)
library(rms)
library(survsim) #package to simulate survival data
N=100 #number of patients
set.seed(123)
df.tf<-simple.surv.sim(#baseline time fixed
  n=N, foltime=500,
  dist.ev=c('llogistic'),
  anc.ev=c(0.68), beta0.ev=c(5.8),
  anc.cens=1.2,
  beta0.cens=7.4,
  z=list(c("unif", 0.8, 1.2)),
  beta=list(c(-0.4),c(0)),
  x=list(c("bern", 0.5),
         c("normal", 70, 13)))
names(df.tf)[c(1,6,7)]<-c("id","grp","age")
nft<-sample(1:10, N,replace=TRUE)#number of follow up time points
crp<-round(abs(rnorm(sum(nft)+N,
                     mean=100,sd=40)),1)
time<-NA
id<-NA
i=0
for(n in nft){
  i=i+1
  time.n<-sample(1:500,n)
  time.n<-c(0,sort(time.n))
  time<-c(time,time.n)
  id.n<-rep(i,n+1)
  id<-c(id,id.n)
}
df.td <- cbind(data.frame(id,time)[-1,],crp) #time-varying covariate
df<-tmerge(df.tf,df.tf,id=id,
           endpt=event(stop,status)) 
df <- tmerge(df,df.td,id=id,
             crp=tdc(time,crp))
df <-df[,c(1,6:11)] #dataset to be used that includes time-varying covariate 
fit.tdc <- coxph(Surv(tstart,tstop,endpt)~
                   grp+age+crp+cluster(id),df)
rmstvc <- cph(Surv(tstart,tstop,endpt)~
                grp+age+crp+cluster(id), x=TRUE, y=TRUE, surv = TRUE, data = df)
validate(rmstvc, method = "boot", B = 5) #validate function worked:
          index.orig training   test optimism index.corrected n
Dxy       0.1068   0.1399 0.0737   0.0662          0.0406 5
R2        0.0146   0.0277 0.0116   0.0161         -0.0015 5
Slope     1.0000   1.0000 0.6834   0.3166          0.6834 5
D         0.0060   0.0129 0.0043   0.0085         -0.0026 5
U        -0.0037  -0.0038 0.0031  -0.0068          0.0031 5
Q         0.0097   0.0166 0.0012   0.0154         -0.0057 5
g         0.2762   0.3848 0.2451   0.1397          0.1365 5
calibrate(rmstvc) #this however returned an error:
Error in reliability[, "index.corrected"] : subscript out of bounds
In addition: There were 50 or more warnings (use warnings() to see the first 50)

Does anyone have any insight of this problem? What are some workarounds?

Edit Per my discussion with @EdM about counting process models, I fit the data with the aalen function from the timereg package. Help is needed with the intepretation of the outout:

fit<-aalen(Surv(tstart, tstop, endpt) ~ grp + age + crp, df, max.time=500, n.sim = 100)
summary(fit)

Additive Aalen Model 

Test for nonparametric terms 

Test for non-significant effects 
            Supremum-test of significance p-value H_0: B(t)=0
(Intercept)                          1.61                0.60
grp                                  2.88                0.03
age                                  2.27                0.24
crp                                  1.61                0.61

Test for time invariant effects 
                  Kolmogorov-Smirnov test p-value H_0:constant effect
(Intercept)                       0.49900                        0.64
grp                               0.34300                        0.23
age                               0.00611                        0.66
crp                               0.00359                        0.20
                    Cramer von Mises test p-value H_0:constant effect
(Intercept)                      2.62e+01                        0.58
grp                              1.82e+01                        0.12
age                              3.84e-03                        0.60
crp                              9.21e-04                        0.33

  Call: 
aalen(formula = Surv(tstart, tstop, endpt) ~ grp + age + crp, 
    data = df, max.time = 500, n.sim = 100)

Answer 1

The validate.cph() function doesn't use survival predictions per se. It fits a model to each of multiple bootstrapped samples from the data and mostly uses likelihood-based evaluations of how well each model works on the corresponding bootstrap sample and the full data set. It does provide a first pass at calibration, insofar as it also evaluates slope optimism by trying to fit the outcomes of the entire data set against the linear predictors derived from each bootstrap-based model.

The more detailed model calibration provided by calibrate.cph() requires estimates of survival probabilities at a specified time. It only works with right-censored survival times, not with the counting-process data format used to handle time-varying covariates and the recurrent events that are of interest to you.

There is a good deal of reluctance in making Cox model predictions based on time-varying covariates. This page and its links provide an introduction to the issues. I suppose some of that reluctance might be alleviated by recurrent-event modeling when there is no final absorbing state, particularly if the covariates are externally imposed rather than internal to each individual (and thus potentially associated with the individual's event history).

Some software does allow predictions from Cox models with counting-process data, as noted in this answer. If your model involves individual-specific frailty terms for recurrent events, however, I don't think that the predictions will take that into account.

You might consider the type of approach used in calibrate.cph(), where model predictions are compared against "observed" survival probabilities based on highly flexible modeling; for recurrent events, you could model the number of events instead of survival probabilities. But make sure that, in your application, you are modeling something that makes sense and doesn't ultimately involve circular reasoning.