I am having problems with the "survfit" method to calculate survival probabilities following fitting an Anderson-Gill (AG) model for recurrent events using the "cph" method in the "rms" package. The problem is that I request the survival probabilities to be calculated at the observed event times, but not all the times and their probabilities are returned, and conversely some times and their probabilities are returned for times that do not exist in the input dataset. It appears a similar question has been asked in this forum before here https://stats.stackexchange.com/questions/172562/problems-with-the-surv-vector-from-survfit-in-r-mising-aleatory-values1 but not resolved.
My data are in counting process format (below are fictitious data)
where the above data show one subject with an outcome event at $T=18$ and $T=55$ and who was censored at $T=365$. The records with FALSE are either the censoring record for that subject, or "created records" corresponding to the observed censoring times or event times of the whole sample (here $T=10,15,32$ but in real applications this list becomes very long as the number of recurrent events and/or subjects increases).
The purpose of creating these records is so that I can calculate survival probabilities for each subject at all the outcome and censoring event times for the whole sample for use in an Inverse Probability Censoring Weighting (IPCW) analysis. Thus I simply wish to calculate survival probabilities for each "tstop" time for each subject where the data are in the above format.
However I might get a survival probability estimated at $T=7$ which doesn't exist, or the survival probabilities at all $T> 32$ where clearly observed times exist beyond this time. I was wondering if perhaps this might be due to the method employed to deal with tied event times. In this respect the AG model is estimated by default using the "efron" method, and I believe the call to "Survfit" then uses this same method which it detects through the cox object passed as an argument.
My R code to fit the model is as follows where cov1-cov3 are my covariates and subjid is my subject index;
AG_mod = cph(Surv(tstart,tstop,censored)~cov1+cov2+cov3+cluster(subjid),x=TRUE,y=TRUE,data=mydata,surv=TRUE,singular.ok=TRUE)
Then my pseudo-code to estimate the survival probabilities is
for each subject #get the data for ith subject datai = data for subject i #get the time intervals as a survival object intervals = Surv(dati$tstart,dati$tstop,dati$censored) #get the covariate data covsi =dati[,c("subjid","cov1","cov2","cov3"),with=FALSE] #combine intervals with covariate data covsi = data.frame(covsi,intervals) #estimate the survival curve for this subject surv_fit = survfit(AG_mod,newdata=covsi,conf.type="none",se.fit=FALSE,id=covsi$subjid,censor=TRUE) end
Finally I do have time-varying covariates in my AG model which might somehow be causing the problem here since this appears to be non-standard. In general I wish to know how to correctly get R to calculate survival probabilities in this situation (be that with the "rms" package or otherwise). In this respect I currently have an open question with a bounty here Using R to calculate survival probabilities with time-varying covariates using an Andersen-Gill model