Cómo colocar una tabla de "número en riesgo" debajo de un diagrama de Kaplan-Meier usando ggplot2
Me gustaría crear un diagrama de Kaplan-Meier usando ggplot2 con un número en la tabla de riesgo debajo que indica el número en riesgo para cada grupo en cada punto de tiempo (es decir, la marca del eje x). El número en riesgo debe estar alineado con la marca correspondiente. A la izquierda de la tabla del número en riesgo deben aparecer los nombres de las filas que indican el grupo al que pertenecen los números en riesgo.
Escribí el siguiente ejemplo. Aprendo cómo determinar los números en riesgo de estopregunta. Sin embargo, no sé cómo crear una buena tabla de riesgo bien alineada debajo del diagrama de Kaplan-Meier. Un amigo me ayudó a crear la tabla de número de riesgos en el siguiente ejemplo. Sin embargo, la cifra resultante de mi ejemplo es insuficiente.
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)
Como sugerí, usé gtable con la figura resultante. No estaba satisfecho con el diseño de la variante a (código de ejemplo de baptiste), así que probé algo más. Sin embargo, la versión B tiene otro inconveniente: las etiquetas están dentro de las dimensiones x de la capa de trazado del trazado principal.
a) ¿Cómo puedo crear una figura distribuida razonable con números de riesgo bien alineados?
b) Además, ¿cómo puedo colocar un título "Números en riesgo" entre la trama principal y la tabla? El título "Números en riesgo" debe estar alineado con el extremo izquierdo de las etiquetas "Grupo A" y "Grupo B" detbl.
c) El tamaño de fuente de los números de riesgo en tbl y las etiquetas correspondientes "Grupo A" y "Grupo B" deben ser las mismas que las etiquetas de marca en la gráfica principal. ¿Cómo puedo hacer esto?
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)
Ahora es casi perfecto. Dos pequeñas cosas:
a) ¿Cómo puedo agregar un título (saber hecho con GIMP) "Número en riesgo" a la trama como se muestra en la figura a continuación?
b) ¿Por qué está el Grupo B en la tabla anterior al Grupo A? La etiqueta en df_nums para el Grupo A es 1 y para el Grupo B 2. ¿Cómo puedo establecer el Grupo A sobre el Grupo B en la tabla de números en riesgo?
> str(df_nums$variable)
Factor w/ 2 levels "Group.A","Group.B": 1 1 1 1 1 1 2 2 2 2 ...
Aquí el código actualizado:
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)
