Como colocar uma tabela "número em risco" abaixo de um gráfico de Kaplan-Meier usando ggplot2
Gostaria de criar um gráfico de Kaplan-Meier usando ggplot2 com um número em tabela de risco abaixo, indicando o número em risco de cada grupo em cada momento (isto é, marca do eixo x). O número em risco deve ser alinhado ao tick correspondente. À esquerda da tabela de número em risco deve haver nomes de linhas indicando o grupo ao qual os números em risco pertencem.
Eu escrevi o seguinte exemplo. Eu aprendo como determinar os números em risco com issoPergunta, questão. No entanto, não sei como criar um número agradável e bem alinhado na tabela de risco abaixo do gráfico de Kaplan-Meier. Um amigo me ajudou a criar o número da tabela de riscos no exemplo a seguir. No entanto, a figura resultante do meu exemplo é 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 sugerido, usei gtable com a figura resultante. Eu não estava satisfeito com o layout da variante a (código de exemplo do baptiste), então tentei outra coisa. No entanto, a versão B tem outra desvantagem: os rótulos estão dentro das dimensões x da camada de plotagem do plot principal.
a) Como posso criar uma figura de layout razoável com números de risco bem alinhados.
b) Além disso, como posso colocar um título "Números em risco" entre o gráfico principal e a tabela? O título "Números em risco" deve estar alinhado com a extremidade esquerda dos rótulos "Grupo A" e "Grupo B" detbl.
c) O tamanho da fonte dos números de risco em tbl e os rótulos correspondentes "Grupo A" e "Grupo B" devem ser iguais aos rótulos de marca na plotagem principal. Como posso fazer isso?
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)
Agora está quase perfeito. Duas pequenas coisas:
a) Como posso adicionar um título (conhecido feito com o GIMP) "Número em risco" ao gráfico, como mostra a figura abaixo?
b) Por que o Grupo B está na tabela acima do Grupo A? O rótulo em df_nums para o Grupo A é 1 e para o Grupo B 2. Como posso definir o Grupo A acima do Grupo B no número na tabela de risco?
> str(df_nums$variable)
Factor w/ 2 levels "Group.A","Group.B": 1 1 1 1 1 1 2 2 2 2 ...
Aqui o código atualizado:
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)
