如何使用 ggplot2 在 Kaplan-Meier 图下方放置“风险数字”表

Question

我想使用 ggplot2 创建一个 Kaplan-Meier 图，下面有一个风险数字表，指示每个时间点每个组的风险数字（即 x 轴刻度）。处于风险中的数字应与相应的刻度对齐。风险数字表的左侧应该是指示风险数字所属组的行名。

我写了下面的例子。我从这个question 中学习了如何确定有风险的数字。但是，我不知道如何在 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)

更新 #1

按照建议，我将 gtable 与结果图一起使用。我对变体 a 的布局不满意（来自 baptiste 的示例代码），所以我尝试了其他方法。但是，版本 B 确实有另一个缺点：标签位于主图的绘图层的 x 维度内。

a) 如何创建布局合理且风险数字一致的图形。

b) 此外，如何在主图和表格之间放置标题“风险数字”？标题“有风险的数字”应与tbl 的标签“A 组”和“B 组”的左端对齐。

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)

版本 #1（风险数字与主图的 x 轴刻度对齐良好，但布局错误

版本 #2（螺旋对齐但更好的布局）

更新 #2

现在它几乎完美了。两件小事：

a) 我如何在图中添加标题（已知使用 GIMP）“有风险的数字”，如下图所示？

b) 为什么 B 组在 A 组上方的表格中？ df_nums 中 A 组的标签为 1，B 组的标签为 2。如何在风险表中将 A 组设置在 B 组之上？

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

这里是更新的代码：

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)

Answer 1

你可以的

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)