The problem: I am hoping to visualise interaction effects between a continuous and a dichotomous variable from a cox proportional hazard regression. When graphing this with different packages (visreg and rms:ggplot(Predict()), I got bafflingly different looking graphs. What did I miss in these two approaches that generated this disagreement?
Also when I am plotting for hazard ratio (with fun=exp), I am not getting the reference point with the narrowed confidence band as suggested in this discussion ( Different prediction plot from survival coxph and rms cph) (what did I miss here?).
Example:
# Load libs library(survival) library(rms) library(visreg) # Regular survival survobj <- with(lung, Surv(time,status)) # Prepare the variables lung$sex <- factor(lung$sex, levels=1:2, labels=c("Male", "Female")) label(lung$sex) <- "Sex" label(lung$age)<- "Age" # The rms survival ddist <- datadist(lung) options(datadist="ddist") rms_surv_fit <- cph(survobj~rcs(age, 4)*sex, data=lung, x=T, y=T) p <- Predict(rms_surv_fit, age=seq(50, 70, times=20), sex=c("Male", "Female"), fun=exp) ggplot(Predict(rms_surv_fit, age , sex)) Gives this plot 
To plot the hazard ratio (where I am expecting a reference point)
ggplot(Predict(rms_surv_fit, age , sex, fun=exp)) 
# visreg coxph_fit <-coxph(survobj~rcs(age, 4)*sex, data=lung) visreg(coxph_fit, 'age', by='sex', overlay=T, partial=F, band=T, ylim=c(-3, 3) ) Gives very different graphs compared to above two:
I came across suggestions that the meaning of above conditional plot is likely questionable since baseline hazard is not estimated - a contrast plot is herhaps more appropriate as graphed below:
visreg(coxph_fit, 'age', by='sex', overlay=T, partial=F, type='contrast' ) 