I've been following along with this thread, but have faced some unrelated problems due to my working environment (cloud computing server etc.): https://stackoverflow.com/questions/28385509/how-to-plot-a-cox-hazard-model-with-splines.
Before bothering my beloved IT professionals to fix my cloud environment, I was wondering if anyone knows another way to make the kind of graph I'm looking for, a spline of a cox regression stratified by another variable.
So for example if I have these 2 plots, with the reference group being when daily calorie intake == 1075 within each sex.
library(survival)
lung1 <- lung %>% filter(sex == 1)
fit <- coxph(Surv(time, status) ~ pspline(meal.cal, df = 5), data = lung1)
ptemp <- termplot(fit, se = TRUE, plot = FALSE)
percterm <- ptemp$meal.cal
center <- with(percterm, y[x==1075])
ytemp <- percterm$y + outer(percterm$se, c(0, -1.96,1.96), "*")
matplot(percterm$x, exp(ytemp - center), log = "y", type = "l", col = "2", lty = c(1,2,2))
lung2 <- lung %>% filter(sex == 2)
fit <- coxph(Surv(time, status) ~ pspline(meal.cal, df = 5), data = lung2)
ptemp <- termplot(fit, se = TRUE, plot = FALSE)
percterm <- ptemp$meal.cal
center <- with(percterm, y[x==1075])
ytemp <- percterm$y + outer(percterm$se, c(0, -1.96,1.96), "*")
matplot(percterm$x, exp(ytemp - center), log = "y", type = "l", col = "2", lty = c(1,2,2))
But I want 1 plot with 2 lines with one reference, for example when males have calorie intake of 1075.
Thanks for your help!