during package development I have come across peculiarities with the 'predict.coxph' function from the survival package in R. It appears 'predict.coxph' is not exported from the survival package, and that to call it I must use 'survival:::predict.coxph' (note the triple ':::'). I have been using 'predict' from the stats package instead of 'predict.coxph' for a long time, and this has made me think about the implications of both.
The 'predict' function from stats must work out that the object is an coxph object and therefore generate the appropriate linear predictor, but there is no mention of using 'predict.coxph' on the documentation page for 'predict'.
Question 1), how does it do this? Does it somehow find 'predict.coxph' from the survival package, or is it completely seperate?
The 'predict.coxph' is not exported from the survival package, which seems odd, given it would be used widely.
Question 2) why is predict.coxph not exported? Is it a bad idea to use 'survival:::predict.coxph'?
This is of interest because when developing a package reliant on 'predict.coxph', I get a warning message about the use of ':::', so plan to use 'predict' instead, but would like to fully understand how it works.
Reproducible code below. Note that coxph.pred1 and coxph.pred4 give the same answers, coxph.pred2 and coxph.pred3 cannot find the function.
library(survival)
### Load data
myeloma <- myeloma
### Add predictor
myeloma$x1 <- rnorm(nrow(myeloma), 0, 1)
### Fit cox proportional hazards
coxph.obj <- coxph(Surv(futime, death) ~ x1, data = myeloma)
### Check model output
summary(coxph.obj)
### Create predictions for individuals in myeloma dataset
coxph.pred1 <- predict(coxph.obj, type = "lp")
coxph.pred2 <- predict.coxph(coxph.obj, type = "lp")
Error in predict.coxph(coxph.obj, type = "lp") :could not find function "predict.coxph"
coxph.pred3 <- survival::predict.coxph(coxph.obj, type = "lp")
Error: 'predict.coxph' is not an exported object from 'namespace:survival'
coxph.pred4 <- survival:::predict.coxph(coxph.obj, type = "lp")
### Print how many elements are not the same in coxph.pred1 and coxph.pred4
sum(coxph.pred1 != coxph.pred4)
The predict function is generic.
predict
# function (object, ...)
# UseMethod("predict")
# <bytecode: 0x13c66f4c8>
# <environment: namespace:stats>
you can tell that because it calls UseMethod which triggers S3 dispatching in R. This means that predict will call different functions depending on the class() of the object you pass to it. These functions are all called predict.<classname>. You can see all the different types of predict function with methods(predict)
By convention, packages that provide methods for predict to not export the "raw" predict.<classname> functions. Instead they register them as S3 methods. (See here in the survival NAMESPACE file.) So you are meant to call just the predict() function with an object of class "coxph" and R will call predict.coxph for you. You should not use ::: form of functions. Just use
coxph.pred1 <- predict(coxph.obj, type = "lp")