Plot a quadratic relation for a predictor of a cox

I need to plot the relative risk for a quadratic effect in a cox regression. My model looks like this:

cox_mod <- coxph(Surv(time, status) ~
               ph.karno + pat.karno + meal.cal + meal.cal_q,
             data = lung)

Where meal.cal_q is defined as:

lung$meal.cal_q <- lung$meal.cal^2

The plot should consider the coefficients of meal.cal and meal.cal_q and show the relative risk on the y-axis and the meal.cal values on the x-axis. The relative risk should be defined as the risk at a given meal.cal value compared to all of the predictors being at their mean. Additionaly, the plot should include the 95% confidence intervals. The plot should look something like this: Expected plot

If possible, the plot should be a ggplot object so that I can customize it.

I have been reading for hours on the web, but can not figure out how make the described plot and hope someone can help me. I tried it for example with the predict() function:

meal.cal_new <- seq(min(lung$meal.cal, na.rm= TRUE), max(lung$meal.cal, na.rm= TRUE), by= 1)

meal.cal_q_new <- meal.cal_new^2

n <- length(meal.cal_new)

lung_new <- data.frame(ph.karno= rep(mean(lung$ph.karno, na.rm= TRUE), n), pat.karno= rep(mean(lung$pat.karno, na.rm= TRUE), n), meal.cal= meal.cal_new, meal.cal_q = meal.cal_q_new)

predicted_rel_risk <- predict(cox_mod, lung_new,  interval = "confidence")

print(predicted_rel_risk)

Firstly, the predicted values do not include the 95% confidence itnervals. And in addition there are negative values in predicted_rel_risk which in my opinien should not be the case since the minimal relative risk should be zero. Therefore I can not get the desired plot. So all I can do is this:

lung_new$predicted_rel_risk <- predicted_rel_risk

ggplot(lung_new, aes(meal.cal, predicted_rel_risk)) +
geom_smooth(se= TRUE)

The resulting plot does not include the confidence intervals and shows neagtive relative risk. Here is what I get:

Thank you a lot in advance!

The prediction includes negative values since you did not specify that you want to obtain the relative risk (as you stated). Try the following code

predicted_rel_risk <- predict(cox_mod, lung_new,  interval = "confidence",
type= "risk")

This gives you the following plot:

Plot without negative values

In order to get the confidence intervalls as well, you can use bootstrapping. To put it short, this means that a random sample will be drawn from your data and the relative risk will be calculated. This procedure will be repeated 10,000 times, for example. This gives you 10,000 different relative risk values for each value of your predictor. You get the main line for your plot by calculating the mean relative risk for each value of your predictor. To get the condidence intervall, you need to order the relative risks from the smallest to the greatest for each value of your predictor. The 250th (9,750th) relative risk value gives you your lower (upper) ci. Again, it is the 250th (9,750th) value of each predictor value.

Hope this helps.