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(survsim) #package to simulate survival data
N=100 #number of patients
df.tf<-simple.surv.sim(#baseline time fixed
  n=N, foltime=500,
  anc.ev=c(0.68), beta0.ev=c(5.8),
  z=list(c("unif", 0.8, 1.2)),
  x=list(c("bern", 0.5),
         c("normal", 70, 13)))
nft<-sample(1:10, N,replace=TRUE)#number of follow up time points
for(n in nft){
df.td <- cbind(data.frame(id,time)[-1,],crp) #time-varying covariate
df <- tmerge(df,df.td,id=id,
df <-df[,c(1,6:11)] #dataset to be used that includes time-varying covariate 
fit.tdc <- coxph(Surv(tstart,tstop,endpt)~
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)

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

aalen(formula = Surv(tstart, tstop, endpt) ~ grp + age + crp, 
    data = df, max.time = 500, n.sim = 100)
  • 1
    $\begingroup$ On a now-deleted post, Frank Harrell commented as follows: "I should make it clear in the documentation but calibrate doesn't understand time-dependent covariates. As Therneau has stated frequently, estimation of survival probabilities in the presence of time-dependent covariates is not a simple thing to conceptualize." $\endgroup$
    – EdM
    Oct 27, 2022 at 22:00
  • 1
    $\begingroup$ You might want to look at this page and its links for an introduction to the problems with using time-varying covariates to predict survival probabilities. $\endgroup$
    – EdM
    Oct 27, 2022 at 22:05
  • $\begingroup$ Thank you for the information @EdM. So I guess there is no way to do calibrate on cox models with time-dependent covariates then? For my own research, I deal with the presence and absence of certain behaviors and so the recurrence of behaviors is possible, unlike death which is the final end point. So I think the problems with with using time-varying covariates may not apply to my study. $\endgroup$
    – cliu
    Oct 28, 2022 at 13:14
  • $\begingroup$ Look at this page for some discussion about recurrent events and problems of making predictions with time-varying covariates, particularly those that are internal to individuals rather than external. You might need to consider a joint model of covariates and events over time. $\endgroup$
    – EdM
    Oct 28, 2022 at 14:24

1 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.

  • $\begingroup$ Thank you for the explanation @EdM. Now it becomes more clear to me why time-varying covariate cox models are not supported by calibrate. Can you point me to some resources how recurrent-event modeling can be done in r? I do think the covariates in my data are external because I use the partner's behaviors to predict the occurence of behaviors of the target individual. $\endgroup$
    – cliu
    Oct 29, 2022 at 0:38
  • 1
    $\begingroup$ @cliu Section 3.2 of the survival vignette illustrates the Andersen-Gill approach. The counting-process format works for modeling and the survfit() function does allow for predictions with time-varying covariates. Amorim and Cai, Int J Epidemiol 44: 324–333 (2015), compare several approaches, with sample code in an appendix. Ozga et al, BMC Med Res Methodol 18: 2 (2018), cover composite endpoints. $\endgroup$
    – EdM
    Oct 29, 2022 at 13:51
  • 1
    $\begingroup$ @cliu the R Survival task view has a section on recurrent events. The frailtypack and reReg packages might help, beyond what you can do with coxph(). Also review the packages for joint modeling of covariates and outcomes, and for broader multi-state models if you have more than one type of event. $\endgroup$
    – EdM
    Oct 29, 2022 at 14:01
  • $\begingroup$ Are you familiar with the additive hazards model, or the additive Aalen model? It seems to be (although I am not too sure) a counting process model as it assumes intensity for the counting process N(t) conditionally on a group of covariate (Martinussen & Scheike, 2006). I was able to fit this model with the aalen function from the timereg package (see the edit of my original question) but don't quite understand the output. Where are the coefficients? How are the coefficients different from those from the cox model? $\endgroup$
    – cliu
    Oct 30, 2022 at 19:24
  • 1
    $\begingroup$ @cliu I don't use additive hazard models much, and then only with the aareg() function in the standard survival pacakge. Can't help with the timereg output. Would make sense to post that part as a new question, with a link back to this one for context. See this page for a comparison of Aalen additive and Cox models, and a limiting case in which you can find an approximate relationship between them. $\endgroup$
    – EdM
    Oct 30, 2022 at 19:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.