如何使用ggplot2在Kaplan-Meier地图下放置一个“危险数字”表?

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.

我想要创建一个使用ggplot2的Kaplan-Meier图，在风险表下面显示每个组在每个时间点(即x轴标记)的风险值。风险的数字应该与相应的刻度对齐。左边的风险表的数字应该是行名称，表示风险所属的组。

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.

我写了下面的例子。我学习如何从这个问题的风险中确定数字。然而，我不知道如何在Kaplan-Meier计划下的风险表中创建一个良好的对齐的数字。一个朋友帮助我在下面的例子中创建了风险表的数量。然而，我的例子的结果是不够的。

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.

根据建议，我使用gtable的结果。我不满足于变量a的布局(来自baptiste的示例代码)，所以我尝试了其他的方法。然而，版本B有另一个缺点:标签在主情节的图层的x维范围内。

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

a)如何创建合理的layouted图形和良好对齐的风险数字。

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.

b)此外，我如何在主情节和表之间放置一个标题“风险数字”?标题“风险数字”应与标签“A组”和“B组”的左端对齐。

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?

c) tbl和相应标签“A组”和“B组”的风险编号的字体大小应该与主图中的标记标签相同。我该怎么做呢?

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?

a)如何添加标题(知道如何使用GIMP)“风险数字”的情节如下图所示?

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?

b)为什么b组在A组以上?A组的df_nums标签为1,B组为b2。在风险表中，我如何将A组置于B组之上?

> 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)

1 个解决方案

#1

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)