How to place a "number at risk" table beneath a Kaplan-Meier plot using ggplot2

I would like to create a Kaplan-Meier plot using ggplot2 with a number at risk table beneath indicating the number at risk for each group at each time point (i.e. x-axis tick). The number at risk should be aligned to the corresponding tick. Left to the number at risk table should be row names indicating the group to which the numbers at risk belong.

I wrote the following example. I learn how to determine the numbers at risk from this question. However, I do not know how to create a nice, well aligned number at risk table beneath the Kaplan-Meier plot. A friend helped me to create the number of risk table in the following example. However, the resulting figure of my example is insufficient.

library(survival)
library(reshape2)
data(colon)
library(Hmisc)

d <- colon[, Cs(time, status, rx)]
rm(colon)
names(d) <- c("days", "event", "group")
d$group <- ifelse(d$group == "Obs", 1, 2)

fit <- survfit(Surv(days,event)~group, data=d)
diff <- survdiff(Surv(days,event)~group, data=d)

risksets <- with(na.omit(d[, Cs(days, event, group)]), table(group, cut(days, seq(0, max(days), by=365) ) ))
number.at.risk <- sapply(1:nrow(risksets), function(i) Reduce("-",  risksets[i,], init=rowSums(risksets)[i], accumulate=TRUE))
number.at.risk <- data.frame(number.at.risk)
names(number.at.risk) <- c("Group.A", "Group.B")
number.at.risk

###
p.value <- round(1 - pchisq(diff$chisq, 1), digits=4)
p.value <- ifelse(p.value < 0.001, "<0.001", paste("= ", p.value))

d.mortality <- data.frame(time=fit$time, surv=fit$surv, strata=summary(fit, censored=T)$strata)
zeros <- data.frame(time=0, surv=1, strata=unique(d.mortality$strata))
d.mortality <- rbind(d.mortality, zeros)
levels(d.mortality$strata) <- c("Group A", "Group B")
d.mortality$surv <- (1-d.mortality$surv)*100 # event free to events and in %
###
g <- ggplot(d.mortality, aes(time, surv, group=strata)) + 
     geom_step(aes(colour=strata), size=1) +
     theme_bw() + # white background
     theme(
          plot.background = element_blank(), 
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          legend.position="none",
          axis.line = element_line(color = 'black'),
          axis.text.x = element_text(size=15),
          axis.text.y = element_text(size=15),
          axis.title.x = element_text(size=17, hjust=.5, vjust=.25, face="bold"),
          axis.title.y = element_text(size=17, hjust=.5, vjust=1.5, face="bold"),
          plot.title = element_text(size=20, hjust=-.1, vjust=1, face="bold")
     ) +
     scale_y_continuous("Cumulative event rate [%]", limits=c(0, 60)) + 
     scale_x_continuous("Time [years]", limits=c(0, 1825), breaks=seq(0, 1825, 365), labels=c(0, 1, 2, 3, 4, 5)) +
     annotate("text", x = 1000, y = 45, label = "Group A") +
     annotate("text", x = 1000, y = 30, label = "Group B") +
     annotate("text", x = 1000, y = 55, label = paste("P ", p.value, "by log-rank test", collapse=""))

number.at.risk = number.at.risk[1:6,]
df_nums = melt(number.at.risk)
df_nums$year = 1:6
tbl = ggplot(df_nums, aes(x = year, y = factor(variable), colour = variable,label=value)) +
     geom_text(size = 3.5) + theme(panel.grid.major = element_blank(), legend.position = "none") +      theme_bw() + 
     theme(
          plot.background = element_blank(), 
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          legend.position="none",
          axis.line = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks=element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          plot.title = element_blank()
     ) + scale_y_discrete(breaks=c("Group.B","Group.A"), labels=c("Number at Risk\nGroup B", "Group A"))

Layout <- grid.layout(nrow = 2, ncol = 1, heights = unit(c(2, 0.55), c("null", "null")))
 grid.show.layout(Layout)
 vplayout <- function(...) {
    grid.newpage()
    pushViewport(viewport(layout = Layout))
}

subplot <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
mmplot <- function(a, b) {
     vplayout()
     print(a, vp = subplot(1, 1))
     print(b, vp = subplot(2, 1))
 }

 dev.new()
mmplot(g, tbl)

UPDATE #1

As suggested I used gtable with the resulting figure. I was not satisfied with the layout of variant a (example code from baptiste), so I tried something else. However, version B does have another drawback: the labels are within the x-dimensions of the plot layer of the main plot.

a) How can I create reasonable layouted figure with well aligned risk numbers.

b) Moreover, how can I place a title "Numbers at risk" between the main plot and the table? The title "Numbers at risk" should be aligned with the left end of the labels "Group A" and "Group B" of tbl.

c) The font size of the risk numbers in tbl and the corresponding labels "Group A" and "Group B" should be the same as the tick labels in the main plot. How can I do this?

library(survival)
library(reshape2)
data(colon)
library(Hmisc)

d <- colon[, Cs(time, status, rx)]
rm(colon)
names(d) <- c("days", "event", "group")
d$group <- ifelse(d$group == "Obs", 1, 2)

fit <- survfit(Surv(days,event)~group, data=d)
diff <- survdiff(Surv(days,event)~group, data=d)

risksets <- with(na.omit(d[, Cs(days, event, group)]), table(group, cut(days, seq(0, max(days), by=365) ) ))
number.at.risk <- sapply(1:nrow(risksets), function(i) Reduce("-",  risksets[i,], init=rowSums(risksets)[i], accumulate=TRUE))
number.at.risk <- data.frame(number.at.risk)
names(number.at.risk) <- c("Group.A", "Group.B")
number.at.risk

###
p.value <- round(1 - pchisq(diff$chisq, 1), digits=4)
p.value <- ifelse(p.value < 0.001, "<0.001", paste("= ", p.value))

d.mortality <- data.frame(time=fit$time, surv=fit$surv, strata=summary(fit, censored=T)$strata)
zeros <- data.frame(time=0, surv=1, strata=unique(d.mortality$strata))
d.mortality <- rbind(d.mortality, zeros)
levels(d.mortality$strata) <- c("Group A", "Group B")
d.mortality$surv <- (1-d.mortality$surv)*100 # event free to events and in %
###
g <- ggplot(d.mortality, aes(time, surv, group=strata)) + 
     geom_step(aes(colour=strata), size=1) +
#           theme_bw() + # white background
     theme(
          plot.background = element_blank(), 
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          legend.position="none",
          axis.line = element_line(color = 'black'),
          axis.text.x = element_text(size=15),
          axis.text.y = element_text(size=15),
          axis.title.x = element_text(size=17, hjust=.5, vjust=.25, face="bold"),
          axis.title.y = element_text(size=17, hjust=.5, vjust=4, face="bold"),
          plot.title = element_text(size=20, hjust=-.1, vjust=1, face="bold")
     ) +
     scale_y_continuous("Cumulative event rate [%]", limits=c(0, 60)) + 
     scale_x_continuous("Time [years]", limits=c(0, 1825), breaks=seq(0, 1825, 365), labels=c(0, 1, 2, 3, 4, 5)) +
     annotate("text", x = 1000, y = 45, label = "Group A") +
     annotate("text", x = 1000, y = 30, label = "Group B") +
     annotate("text", x = 1000, y = 55, label = paste("P ", p.value, "by log-rank test", collapse=""))

number.at.risk = number.at.risk[1:6,]
df_nums = melt(number.at.risk)
df_nums$year = 1:6
str(df_nums)

tbl <- ggplot(df_nums, aes(x = year, y = factor(variable), colour = variable, label=value)) +
     geom_text() +
#           theme_bw() + 
     theme(
          panel.grid.major = element_blank(), 
          legend.position = "none",
          plot.background = element_blank(), 
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          legend.position="none",
          axis.line = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks=element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          plot.title = element_blank()
     ) + 
     scale_y_discrete(breaks=c("Group.B","Group.A"), labels=c("Group B", "Group A"))

library(gtable)

# Version A
both = rbind(ggplotGrob(g), ggplotGrob(tbl), size="last")
grid.newpage()
grid.draw(both)

# Version B
a <- gtable(unit(15, c("cm")), unit(c(10,3), "cm"))
a <- gtable_add_grob(a, ggplotGrob(g), 1, 1)
a <- gtable_add_grob(a, ggplotGrob(tbl), 2, 1)
grid.newpage()
grid.draw(a)

Version #1 (risk numbers well-aligned to x-axis ticks of main plot but bad layout

Version #2 (screwed alignement but better layout)

UPDATE #2

Now it's almost perfect. Two small things:

a) How can I add a the title (know done with GIMP) "Number at risk" to the plot as shown in the figure below?

b) Why is Group B in the table above Group A? The label in df_nums for Group A is 1 and for Group B 2. How can I set Group A above Group B in the number at risk table?

> str(df_nums$variable)
 Factor w/ 2 levels "Group.A","Group.B": 1 1 1 1 1 1 2 2 2 2 ...

Here the updated code:

library(survival)
library(reshape2)
data(colon)
library(Hmisc)

d <- colon[, Cs(time, status, rx)]
rm(colon)
names(d) <- c("days", "event", "group")
d$group <- ifelse(d$group == "Obs", 1, 2)

fit <- survfit(Surv(days,event)~group, data=d)
diff <- survdiff(Surv(days,event)~group, data=d)

risksets <- with(na.omit(d[, Cs(days, event, group)]), table(group, cut(days, seq(0, max(days), by=365) ) ))
number.at.risk <- sapply(1:nrow(risksets), function(i) Reduce("-",  risksets[i,], init=rowSums(risksets)[i], accumulate=TRUE))
number.at.risk <- data.frame(number.at.risk)
names(number.at.risk) <- c("Group.A", "Group.B")
number.at.risk

###
p.value <- round(1 - pchisq(diff$chisq, 1), digits=4)
p.value <- ifelse(p.value < 0.001, "<0.001", paste("= ", p.value))

d.mortality <- data.frame(time=fit$time, surv=fit$surv, strata=summary(fit, censored=T)$strata)
zeros <- data.frame(time=0, surv=1, strata=unique(d.mortality$strata))
d.mortality <- rbind(d.mortality, zeros)
levels(d.mortality$strata) <- c("Group A", "Group B")
d.mortality$surv <- (1-d.mortality$surv)*100 # event free to events and in %
###
g <- ggplot(d.mortality, aes(time, surv, group=strata)) + 
     geom_step(aes(colour=strata), size=1) +
#           theme_bw() + # white background
     theme(
          plot.background = element_blank(), 
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          legend.position="none",
          axis.line = element_line(color = 'black'),
          axis.text.x = element_text(size=15),
          axis.text.y = element_text(size=15),
          axis.title.x = element_text(size=17, hjust=.5, vjust=.25, face="bold"),
          axis.title.y = element_text(size=17, hjust=.5, vjust=4, face="bold"),
          plot.title = element_text(size=20, hjust=-.1, vjust=1, face="bold")
     ) +
     scale_y_continuous("Cumulative event rate [%]", limits=c(0, 60)) + 
     scale_x_continuous("Time [years]", limits=c(0, 1825), breaks=seq(0, 1825, 365), labels=c(0, 1, 2, 3, 4, 5)) +
     annotate("text", x = 1000, y = 45, label = "Group A") +
     annotate("text", x = 1000, y = 30, label = "Group B") +
     annotate("text", x = 1000, y = 55, label = paste("P ", p.value, "by log-rank test", collapse=""))

number.at.risk = number.at.risk[1:6,]
df_nums = melt(number.at.risk)
str(df_nums$variable)
df_nums
df_nums$year = 1:6
str(df_nums)

tbl <- ggplot(df_nums, aes(x = year, y = factor(variable), colour = variable, label=value)) +
     geom_text() +
#           theme_bw() + 
     theme(
          panel.grid.major = element_blank(), 
          legend.position = "none",
          plot.background = element_blank(), 
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          legend.position="none",
          axis.line = element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_text(size=15, face="bold", color = 'black'),
          axis.ticks=element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          plot.title = element_blank()
     ) + 
     scale_y_discrete(breaks=c("Group.A", "Group.B"), labels=c("Group A", "Group B"))

library(gtable)

# Version C
both = rbind(ggplotGrob(g), ggplotGrob(tbl), size="last")
panels <- both$layout$t[grep("panel", both$layout$name)]
both$heights[panels] <- list(unit(1,"null"), unit(2, "lines"))
grid.newpage()
grid.draw(both)

Solution

you could do

both = rbind(ggplotGrob(g), ggplotGrob(tbl), size="last")
panels <- both$layout$t[grep("panel", both$layout$name)]
both$heights[panels] <- list(unit(1,"null"), unit(2, "lines"))
both <- gtable_add_rows(both, heights = unit(1,"line"), 8)
both <- gtable_add_grob(both, 
                        textGrob("Number at risk", hjust=0, x=0), 
                        t=9, l=2, r=4)
grid.newpage()
grid.draw(both)