Disclaimer: I am new to R and survival analyses so apologies if I don't use the right terms or coding etiquette-figuring it out.

In the carcass package, I've used function persistence.prob to develop estimates of persistence probability for each day following placement. However, the function produces a list with estimates for only times when an event has occurred, and I need to interpolate the missing data to get daily persistence probabilities.

Note that persistence.prob argument pers.const = FALSE is the default and assumes a non-constant persistence probability, thereby using a coxph model to generate persistence estimates.

> vars <- persistence.prob(removal$siteID, removal$time, 
+                          removal$status) 
> str(vars)
List of 4
 $ time            : num [1:7] 1 2 4 6 9 14 29
 $ persistence.prob: num [1:7, 1:3] 0.782 0.614 0.542 0.47 0.397 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:7] "1" "2" "4" "6" ...
  .. ..$ : chr [1:3] "1" "2" "3"
 $ lower           : num [1:7, 1:3] 0.597 0.399 0.33 0.266 0.206 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:7] "1" "2" "4" "6" ...
  .. ..$ : chr [1:3] "1" "2" "3"
 $ upper           : num [1:7, 1:3] 1 0.943 0.888 0.829 0.766 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:7] "1" "2" "4" "6" ...
  .. ..$ : chr [1:3] "1" "2" "3"

I attempted to generate a coxph survival object with my data and extract the fitted values from it using the predict function, but the values varied significantly from those generated by function persistence.prob.

# Create a coxph survival curve

persistence.cph.model <- coxph(Surv(time, status) ~ siteID, data = removal)

# Predict missing values with predict function for all sites

newdata <- data.frame(time = 1:42, status = 1, siteID = c(1:3))

pred <- predict(persistence.cph.model, newdata = newdata, 
                type = "survival", se.fit = TRUE)
newdata$prob <- pred$fit
newdata$lower <- pred$fit - 1.96*pred$se.fit
newdata$upper <- pred$fit + 1.96*pred$se.fit
newdata

Is there any way to extract the fitted model generated by persistence.prob to interpolate the missing data? Or should I continue trying to use the predict function?