r - Fail when combining ggplots with plots produced from another package: why and how to succeed?

Question

I want to combine two ggplots with two calibration plots made in riskRegression.

I tried many approaches, e.g. par(mfrow) (see below), but failed. I reckon there must be a way to combine ggplots with other plots in a par(mfrow=c(2,2))-manner.

Question: why is this approach invalid, and how can it be fixed?

Data sample o can be found below.

I use

library(riskRegression)
library(survival)
library(tidyverse) 

First produce the calibration

i <- Score(list("WHO w/o TERT" = coxph(Surv(event.tid, event==1) ~ WHO, data=o, x=T, y=T), 
                "WHO +TERT   " = coxph(Surv(event.tid, event==1) ~ WHO + tert.mut, data=o, x=T, y=T)), 
           Surv(event.tid, event==1)~1,
           data=o,
           times=c(60,90),
           plots = "calibration",
           summary="risks")

Then the two calibration plot

plotCalibration(i, 
                cens.method = "local", 
                legend.x=.215,
                legend.y=.225,
                cex=1,
                times=60,
                col=c("#6DBCC3","#8B3A62"),
                ylab="Observed frequency",
                xlab="Predicted risk of recurrence",
                auc.in.legend = TRUE)

And

plotCalibration(i, 
                cens.method = "local", 
                legend.x=.215,
                legend.y=.225,
                cex=1,
                times=90,
                col=c("#6DBCC3","#8B3A62"),
                ylab="Observed frequency",
                xlab="Predicted risk of recurrence",
                auc.in.legend = TRUE)

And then the two ggplot

ggplot(arrange(o, event.tid) %>% 
         transmute(ID = 1:nrow(o),
                   WHO = WHO,
                   tert.mut = tert.mut) %>%
         left_join(i$Calibration$plotframe, 
                   ., 
                   by="ID") %>%
         select(ID, times, model, risk, WHO, tert.mut) %>%
         pivot_wider(names_from = model, values_from = risk) %>%
         rename(., "woTERT" = "WHO w/o TERT", "wTERT" = "WHO +TERT   ") %>%
         filter(times==60),
       aes(x=woTERT, y=wTERT, shape=WHO, col=tert.mut, fill=tert.mut)) +
  geom_point(size = 7)

And

ggplot(arrange(o, event.tid) %>% 
         transmute(ID = 1:nrow(o),
                   WHO = WHO,
                   tert.mut = tert.mut) %>%
         left_join(i$Calibration$plotframe, 
                   ., 
                   by="ID") %>%
         select(ID, times, model, risk, WHO, tert.mut) %>%
         pivot_wider(names_from = model, values_from = risk) %>%
         rename(., "woTERT" = "WHO w/o TERT", "wTERT" = "WHO +TERT   ") %>%
         filter(times==90),
       aes(x=woTERT, y=wTERT, shape=WHO, col=tert.mut, fill=tert.mut)) +
  geom_point(size = 7) 

Attempt #1

par(mfrow=c(2,2), mai=c(.65,.75,.3,.3))
# the first two goes well
plotCalibration(i, 
                cens.method = "local", 
                legend.x=.215,
                legend.y=.225,
                cex=1,
                times=60,
                col=c("#6DBCC3","#8B3A62"),
                ylab="Observed frequency",
                xlab="Predicted risk of recurrence",
                auc.in.legend = TRUE)

plotCalibration(i, 
                cens.method = "local", 
                legend.x=.215,
                legend.y=.225,
                cex=1,
                times=90,
                col=c("#6DBCC3","#8B3A62"),
                ylab="Observed frequency",
                xlab="Predicted risk of recurrence",
                auc.in.legend = TRUE)

# here it goes wrong

ggplot(arrange(o, event.tid) %>% 
         transmute(ID = 1:nrow(o),
                   WHO = WHO,
                   tert.mut = tert.mut) %>%
         left_join(i$Calibration$plotframe, 
                   ., 
                   by="ID") %>%
         select(ID, times, model, risk, WHO, tert.mut) %>%
         pivot_wider(names_from = model, values_from = risk) %>%
         rename(., "woTERT" = "WHO w/o TERT", "wTERT" = "WHO +TERT   ") %>%
         filter(times==60),
       aes(x=woTERT, y=wTERT, shape=WHO, col=tert.mut, fill=tert.mut)) +
  geom_point(size = 7) 

ggplot(arrange(o, event.tid) %>% 
         transmute(ID = 1:nrow(o),
                   WHO = WHO,
                   tert.mut = tert.mut) %>%
         left_join(i$Calibration$plotframe, 
                   ., 
                   by="ID") %>%
         select(ID, times, model, risk, WHO, tert.mut) %>%
         pivot_wider(names_from = model, values_from = risk) %>%
         rename(., "woTERT" = "WHO w/o TERT", "wTERT" = "WHO +TERT   ") %>%
         filter(times==90),
       aes(x=woTERT, y=wTERT, shape=WHO, col=tert.mut, fill=tert.mut)) +
  geom_point(size = 7) 

Attemp #2 was producing a ggarrange to combine with par(mfrow), which also failed.

Expected output

Data sample

o <- structure(list(event.tid = c(56, 127, 29.6, 18, 13.2, 41, 96, 
146.83, 45.37, 1, 105.29, 55, 13.84, 12, 26, 10.8, 1, 7.2, 22.68, 
59, 16, 8.4, 75.3, 25.7, 13, 21, 19, 22, 48, 93, 26, 72, 73.5, 
125, 31, 63, 2.02, 28.2, 2, 8, 60, 46.7, 53, 4, 85, 12, 75.9, 
1, 5, 2.5), event = c(0, 0, 1, 1, 1, 0, 0, 2, 1, 2, 0, 0, 1, 
1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 2, 0, 
0, 0, 1, 0, 0, 1, 2, 0, 0, 1, 0, 1, 0, 2, 0, 1), WHO = structure(c(3L, 
3L, 3L, 2L, 3L, 2L, 3L, 2L, 3L, 3L, 3L, 1L, 2L, 1L, 2L, 3L, 2L, 
3L, 2L, 3L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 
3L, 2L, 1L, 3L, 2L, 2L, 2L, 1L, 2L, 3L, 2L, 2L, 3L, 1L, 2L, 3L, 
2L), .Label = c("WHO-I", "WHO-II", "WHO-III"), class = "factor"), 
    tert.mut = structure(c(1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L), .Label = c("0", 
    "1"), class = "factor")), row.names = c(NA, -50L), class = "data.frame")

question from:https://stackoverflow.com/questions/65861547/fail-when-combining-ggplots-with-plots-produced-from-another-package-why-and-ho

Answer 1

Here is one approach using cowplot.

You can use as_grob to convert your plotCalibration plots to grobs that can be combined with your 2 ggplots.

You can further customize your calibration plots by adding base R graphics statements to your functions.

With plot_grid from cowplot you can put your 2x2 plots in arranged form, in this case with 2 rows + 2 columns.

library(riskRegression)
library(survival)
library(tidyverse) 
library(cowplot)

i <- Score(list("WHO w/o TERT" = coxph(Surv(event.tid, event==1) ~ WHO, data=o, x=T, y=T), 
                "WHO +TERT   " = coxph(Surv(event.tid, event==1) ~ WHO + tert.mut, data=o, x=T, y=T)), 
           Surv(event.tid, event==1)~1,
           data=o,
           times=c(60,90),
           plots = "calibration",
           summary="risks")

make_pc1 <- function() {
  par(mgp=c(4,1,0), mar=c(5,6,4,2)+0.1, cex.lab=1.4, cex.axis=1.5)
  plotCalibration(i, 
                  cens.method = "local", 
                  legend.x=.215,
                  legend.y=.225,
                  cex=1,
                  times=60,
                  col=c("#6DBCC3","#8B3A62"),
                  ylab="Observed frequency",
                  xlab="Predicted risk of recurrence",
                  auc.in.legend = TRUE)
  mtext("10-year prediction", adj=0.05, line = -1.4, cex=1, font=2)
}

pc1 <- as_grob(make_pc1)

make_pc2 <- function() {
  par(mgp=c(4,1,0), mar=c(5,6,4,2)+0.1, cex.lab=1.4, cex.axis=1.5)
  plotCalibration(i, 
                  cens.method = "local", 
                  legend.x=.215,
                  legend.y=.225,
                  cex=1,
                  times=90,
                  col=c("#6DBCC3","#8B3A62"),
                  ylab="Observed frequency",
                  xlab="Predicted risk of recurrence",
                  auc.in.legend = TRUE)
  mtext("10-year prediction", adj=0.05, line = -1.4, cex=1, font=2)
}

pc2 <- as_grob(make_pc2)

gg1 <- ggplot(arrange(o, event.tid) %>% 
         transmute(ID = 1:nrow(o),
                   WHO = WHO,
                   tert.mut = tert.mut) %>%
         left_join(i$Calibration$plotframe, 
                   ., 
                   by="ID") %>%
         select(ID, times, model, risk, WHO, tert.mut) %>%
         pivot_wider(names_from = model, values_from = risk) %>%
         rename(., "woTERT" = "WHO w/o TERT", "wTERT" = "WHO +TERT   ") %>%
         filter(times==60),
       aes(x=woTERT, y=wTERT, shape=WHO, col=tert.mut, fill=tert.mut)) +

  geom_point(size = 7)

gg2 <- ggplot(arrange(o, event.tid) %>% 
         transmute(ID = 1:nrow(o),
                   WHO = WHO,
                   tert.mut = tert.mut) %>%
         left_join(i$Calibration$plotframe, 
                   ., 
                   by="ID") %>%
         select(ID, times, model, risk, WHO, tert.mut) %>%
         pivot_wider(names_from = model, values_from = risk) %>%
         rename(., "woTERT" = "WHO w/o TERT", "wTERT" = "WHO +TERT   ") %>%
         filter(times==90),
       aes(x=woTERT, y=wTERT, shape=WHO, col=tert.mut, fill=tert.mut)) +
  geom_point(size = 7) 

# Draw all 4 plots
plot_grid(pc1, pc2, gg1, gg2, ncol = 2)

Plot

Tweaking the spacing of the calibration plots a bit (e.g., using title for the axis labels, and adjusting mgp for margin line of axis title), you can do:

make_pc1 <- function() {
  par(mar=c(5,6,4,2)+0.1, cex.lab=1.4, cex.axis=1.5)
  plotCalibration(i, 
                  cens.method = "local", 
                  legend.x=.215,
                  legend.y=.225,
                  cex=1,
                  times=60,
                  col=c("#6DBCC3","#8B3A62"),
                  ylab="",
                  xlab="",
                  auc.in.legend = TRUE)
  mtext("10-year prediction", adj=0.05, line = -1.4, cex=1, font=2)
  title(ylab = "Observed frequency", mgp = c(4.5,1,0))
  title(xlab = "Predicted risk of recurrence", mgp = c(2.5,1,0))
}

pc1 <- as_grob(make_pc1)

make_pc2 <- function() {
  par(mar=c(5,6,4,2)+0.1, cex.lab=1.4, cex.axis=1.5)
  plotCalibration(i, 
                  cens.method = "local", 
                  legend.x=.215,
                  legend.y=.225,
                  cex=1,
                  times=90,
                  col=c("#6DBCC3","#8B3A62"),
                  ylab="",
                  xlab="",
                  auc.in.legend = TRUE)
  mtext("10-year prediction", adj=0.05, line = -1.4, cex=1, font=2)
  title(ylab = "Observed frequency", mgp = c(4.5,1,0))
  title(xlab = "Predicted risk of recurrence", mgp = c(2.5,1,0))
}

pc2 <- as_grob(make_pc2)

# Draw all 4 plots
plot_grid(pc1, pc2, gg1, gg2, ncol = 2)

Plot