Adding inset plots using the cowplot package

casebase_inset.R

library(casebase)
library(survival)
library(splines)
library(tidyverse)
library(cowplot)

# 1. Fit casebase with splines----
data("ERSPC")
ERSPC <- mutate(ERSPC, ScrArm = factor(ScrArm, 
                                       levels = c(0,1), 
                                       labels = c("Control group", "Screening group"))
)

fit <- fitSmoothHazard(DeadOfPrCa ~ bs(Follow.Up.Time) + ScrArm, 
                       data = ERSPC, ratio = 100)

# Absolute risk
new_data <- data.frame(ScrArm = c("Control group", 
                                  "Screening group"))
new_time <- seq(0, 14, by = 0.1)

risk_mp <- absoluteRisk(fit, time = new_time, 
                        newdata = new_data)

# Confidence bands
conf_ints <- confint(risk_mp, fit)

# 2. Fit Kaplan-Meier----
# Note: Done separately to get separate objects

# Control group
km_fit <- survfit(Surv(Follow.Up.Time, DeadOfPrCa) ~ 1, 
                  data = ERSPC,
                  subset = ScrArm == "Control group",
                  conf.type = "log-log")

km_ctrl <- tibble(
    time = km_fit$time,
    estimate = 1 - km_fit$surv,
    conf.low = 1 - km_fit$lower,
    conf.upper = 1 - km_fit$upper
) |> 
    pivot_longer(cols = c(estimate, conf.low, conf.upper),
                 names_to = "type",
                 values_to = "estimate") |> 
    filter(time <= 14)

# Screening group
km_fit <- survfit(Surv(Follow.Up.Time, DeadOfPrCa) ~ 1, 
                  data = ERSPC,
                  subset = ScrArm == "Screening group",
                  conf.type = "log-log")

km_scr <- tibble(
    time = km_fit$time,
    estimate = 1 - km_fit$surv,
    conf.low = 1 - km_fit$lower,
    conf.upper = 1 - km_fit$upper
) |> 
    pivot_longer(cols = c(estimate, conf.low, conf.upper),
                 names_to = "type",
                 values_to = "estimate") |> 
    filter(time <= 14)

# 3. Create plots with insets----
range_y <- c(0, 0.015)
myCols <- c("#009E73", "#D55E00")

# Control group
full_ctrl <- conf_ints |> 
    filter(cov_prof == "Control group") |> 
    ggplot() +
    geom_ribbon(aes(x = time,
                    ymin = conf.low,
                    ymax = conf.high),
                fill = "lightgray") +
    geom_step(data = km_ctrl,
              aes(x = time, y = estimate, group = type,
                  linetype = type),
              colour = myCols[2]) +
    geom_line(aes(x = time, y = estimate),
              colour = myCols[1]) +
    theme_half_open() +
    scale_y_continuous(labels = scales::percent) +
    scale_x_continuous(breaks = seq(0, 14, by = 2)) +
    scale_linetype_manual(values = c("dotted", "dotted", "solid"),
                          guide = "none") +
    labs(title = "Control group",
         x = "Years since randomization",
         y = "Cumulative incidence (%)") +
    coord_cartesian(ylim = range_y)

inset_ctrl <- conf_ints |> 
    filter(cov_prof == "Control group",
           time >= 8, time <= 14) |> 
    ggplot() +
    geom_ribbon(aes(x = time,
                    ymin = conf.low,
                    ymax = conf.high),
                fill = "lightgray") +
    geom_step(data = filter(km_ctrl,
                            time >= 8, 
                            time <= 14),
              aes(x = time, y = estimate, group = type,
                  linetype = type),
              colour = myCols[2]) +
    geom_line(aes(x = time, y = estimate),
              colour = myCols[1]) +
    theme_minimal_hgrid(10) +
    scale_y_continuous(labels = scales::percent,
                       expand = c(0, 0)) +
    scale_x_continuous(breaks = seq(8, 14, by = 2)) +
    scale_linetype_manual(values = c("dotted", "dotted", "solid"),
                          guide = "none") +
    labs(x = "", y = "") + expand_limits(y = 0) +
    coord_cartesian(ylim = range_y)

panel_ctrl <- ggdraw(full_ctrl +
                         draw_plot(inset_ctrl, 0, 0.005, 
                                   8, 0.01)
                     )

# Screening group
full_scr <- conf_ints |> 
    filter(cov_prof == "Screening group") |> 
    ggplot() +
    geom_ribbon(aes(x = time,
                    ymin = conf.low,
                    ymax = conf.high),
                fill = "lightgray") +
    geom_step(data = km_scr,
              aes(x = time, y = estimate, group = type,
                  linetype = type),
              colour = myCols[2]) +
    geom_line(aes(x = time, y = estimate),
              colour = myCols[1]) +
    theme_half_open() +
    scale_y_continuous(labels = scales::percent) +
    scale_x_continuous(breaks = seq(0, 14, by = 2)) +
    scale_linetype_manual(values = c("dotted", "dotted", "solid"),
                          guide = "none") +
    labs(title = "Screening group",
         x = "Years since randomization",
         y = "") +
    coord_cartesian(ylim = range_y)

inset_scr <- conf_ints |> 
    filter(cov_prof == "Screening group",
           time >= 8, time <= 14) |> 
    ggplot() +
    geom_ribbon(aes(x = time,
                    ymin = conf.low,
                    ymax = conf.high),
                fill = "lightgray") +
    geom_step(data = filter(km_scr,
                            time >= 8, 
                            time <= 14),
              aes(x = time, y = estimate, group = type,
                  linetype = type),
              colour = myCols[2]) +
    geom_line(aes(x = time, y = estimate),
              colour = myCols[1]) +
    theme_minimal_hgrid(10) +
    scale_y_continuous(labels = scales::percent,
                       expand = c(0, 0)) +
    scale_x_continuous(breaks = seq(8, 14, by = 2)) +
    scale_linetype_manual(values = c("dotted", "dotted", "solid"),
                          guide = "none") +
    labs(x = "", y = "") + expand_limits(y = 0) +
    coord_cartesian(ylim = range_y)

panel_scr <- ggdraw(full_scr +
                         draw_plot(inset_scr, 0, 0.005, 
                                   8, 0.01)
)

# Combine plots
plot_grid(panel_ctrl, panel_scr)

And here's the output:

The plot would still require a legend (orange is Kaplan-Meier; green is case-base).