Dear R-help,

I am using R version 3.4.0 within Windows, and survival 2.41-3. I have fit a Prentice Williams and Peterson-Counting Process model to my data as shown below. This is basically an extension of the Cox model for interval censored data. My dataset, bdat5 can be found here:

https://drive.google.com/open?id=1sQSBEe1uBzh_gYbcj4P5Kuephvalc3gh
cfitcp2 <- coxph(Surv(start,stop,status)~sex+rels+factor(treat)+log(age)+log(tcrate3+0.01)+cluster(trialno)+strata(enum),data=bdat5,model=TRUE,x=TRUE,y=TRUE)

I would now like to use the model to predict the probability of zero events by two years - this is equivalent to the survival probability at 2 years I believe. This is so that I can compare the output to similar estimates obtained from negative binomial, and zero-inflated negative binomial models for the same data (albeit in a different format)

To my mind, and based on what I've read, the best way to do this is to use survfit. I want to make predictions for each individual, therefore, I have tried this code:

trialnos <- unique(bdat5$trialno)

prob0 <- function(ids,dataset,model,time){

probs <- rep(0,length(ids))

for(i in 1:length(ids)){

print(i)

sdata <- subset(dataset,trialno==ids[i])

sfit <- survfit(model,newdata=sdata)

probs[i] <-sum(summary(sfit,time)$surv)

}

return(probs)

}

prob0ests <- prob0(trialnos,bdat5,cfitcp2,730)

When I do this for the first three trial numbers I get:

0.3001021 2993.4531767 0.3445589

The unusually large "probability" arises when there is only 1 row of data for the relevant trial number. Can anyone therefore explain why there is a problem when "sdata" is only 1 row, and ideally provide a solution?

Many thanks,

Laura

Dr Laura Bonnett

NIHR Post-Doctoral Fellow

Department of Biostatistics,

Waterhouse Building, Block F,

1-5 Brownlow Street,

University of Liverpool,

Liverpool,

L69 3GL

0151 795 9686

[hidden email]
______________________________________________

[hidden email] mailing list -- To UNSUBSCRIBE and more, see

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.