title: Hazard ratio plot with confidence "waist" in ggplot2
body: <p>When fitting a cox model that includes spline terms for a continuous covariate, I would like to be able to produce a plot of the hazard ratio across range of that covariate (relative to a fixed reference value) using ggplot2.</p>
  
  <p>I have adapted an example from Terry Therneau's splines vignette <a href="https://cran.r-project.org/web/packages/survival/vignettes/splines.pdf" rel="nofollow noreferrer">here</a> (see page 3). The only issue with this approach is the lack of a "waist" in the confidence interval at the reference value, as in this plot:</p>
  
  <p><a href="https://i.stack.imgur.com/AzAzC.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/AzAzC.png" alt="enter image description here"></a> </p>
  
  <p>The example below produces the following plot, without the narrowing of the CI at the reference value.</p>
  
  <pre><code>library(survival)
  library(splines)
  library(ggplot2)
  
  # colon cancer death dataset
  ccd &lt;- na.omit(subset(colon, etype == 2))
  
  # fit model with ns() term for age
  cox &lt;- coxph(Surv(time, status) ~ rx + sex + ns(age, knots = c(20, 50, 70)), data = ccd)
  
  # get data for plot
  tp &lt;- termplot(cox, se = TRUE, plot = FALSE)
  
  # hazard ratio plot for natural spline of age, with reference @ 50 yrs
  ref &lt;- tp$age$y[tp$age$x == 50]
  
  ggplot() +
  
    geom_line(data = tp$age, aes(x = x, y = exp(y - ref))) + 
    geom_line(data = tp$age, aes(x = x, y = exp(y - 1.96 * se - ref)), linetype = 2) + 
    geom_line(data = tp$age, aes(x = x, y = exp(y + 1.96 * se - ref)), linetype = 2) +
  
    geom_hline(aes(yintercept = 1), linetype = 3) +
  
    geom_rug(data = ccd, aes(x = age), sides = "b") +
  
    labs(x = "Age at baseline, years",
         y = "Hazard Ratio (95% CI) vs. 50 years",
         title = "Mortality hazard ratio as a function of age",
         subtitle = "Natural spline: knots at 20, 50, and 70 years")
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/k7hvp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/k7hvp.png" alt="enter image description here"></a></p>
  
  <p>I am aware that there are features in the <code>rms</code> package and the <code>smoothHR</code>package that produce these types of plots, but I am looking for a solution that is amenable to <code>ggplot2</code> graphics and the <code>coxph()</code> function in the <code>survival</code> package. My question therefore boils down to:</p>
  
  <ul>
  <li>Is there a way to adapt the output of <code>termplot()</code> to produce a plot with a "waist" at the reference value?</li>
  <li>If <code>termplot()</code> cannot be used, how can I obtain the relevant plotting data by other means?</li>
  </ul>
  
  <hr>
  
  <p><strong>Edit 1:</strong> As the first comment suggested, this can be accomplished using <code>rms</code> and <code>ggplot2</code> together. For example:</p>
  
  <pre><code>    library(rms)
  dd &lt;- datadist(ccd)
  dd$limits$age[2] &lt;- 50
  options(datadist = "dd")
  
  cph &lt;- cph(Surv(time, status) ~ rx + sex + rcs(age, c(20, 50, 70)), data = ccd, x = TRUE, y = TRUE)
  
  pdata &lt;- Predict(cph, age, ref.zero = TRUE, fun = exp)
  
  ggplot(data = pdata) +
  
    geom_hline(aes(yintercept = 1), linetype = 3) +
  
    labs(x = "Age at baseline, years",
         y = "Hazard Ratio (95% CI) vs. 50 years",
         title = "Mortality hazard ratio as a function of age",
         subtitle = "Natural spline: knots at 20, 50, and 70 years")
  </code></pre>
  
  <p>Which produces a plot very close to what I am after:</p>
  
  <p><a href="https://i.stack.imgur.com/9x86A.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/9x86A.png" alt="enter image description here"></a></p>
  
  <p>However, <strong>I would still like to know if there is a way to do this using <code>coxph()</code> and <code>ns()</code></strong>. Not that I have anything against the <code>rms</code> package, I just have a bunch of old code based on <code>survival</code>functionality.</p>