Created
September 30, 2021 20:37
-
-
Save turgeonmaxime/62153387a1699b416860931abb73cb8d to your computer and use it in GitHub Desktop.
Adding inset plots using the cowplot package
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
And here's the output:
![image](https://user-images.githubusercontent.com/6857788/135526870-701b2095-8337-4faf-8676-2a9ff77cc1c3.png)
The plot would still require a legend (orange is Kaplan-Meier; green is case-base).