読者です 読者をやめる 読者になる 読者になる

論文風な Kaplan-Meier plot を書いてみよう

R ggplot2

追記

マウス操作のみで作図が可能な R commander 用プラグインが CRAN から利用できます!
RcmdrPlugin.KMggplot2_0.1-0 is on CRAN now - Triad sou.


テスト版

require(ggplot2)
require(survival)

`kmg2.theme_gray` <-
function (base_size = 12, base_family = "") {
  structure(list(
    axis.line = theme_blank(),
    axis.text.x = theme_text(family = base_family, size = base_size * 0.8,
      lineheight = 0.9, colour = "grey50", vjust = 1),
    axis.text.y = theme_text(family = base_family, size = base_size * 0.8,
      lineheight = 0.9, colour = "grey50", hjust = 1),
    axis.ticks = theme_segment(colour = "grey50"),
    axis.title.x = theme_text(family = base_family, size = base_size, vjust = -0.5),
    axis.title.y = theme_text(family = base_family, size = base_size, angle = 90,
      vjust = 0.375),
    axis.ticks.length = unit(0.15, "cm"),
    axis.ticks.margin = unit(0.2, "cm"),
    legend.background = theme_rect(fill = "grey95", colour = "white"),
    legend.key = theme_rect(fill = "grey95", colour = "white"),
    legend.key.size = unit(1.5, "lines"),
    legend.text = theme_text(family = base_family, size = base_size * 0.6),
    legend.title = theme_text(family = base_family, face = "bold", size = base_size * 0.6,
      hjust = 0),
    legend.position = c(0.95, 0.95),
    legend.direction = "vertical",
    panel.background = theme_rect(fill = "grey90", colour = NA),
    panel.border = theme_blank(),
    panel.grid.major = theme_line(colour = "white"),
    panel.grid.minor = theme_line(colour = "grey95", size = 0.25),
    panel.margin = unit(0.25, "lines"),
    strip.background = theme_rect(fill = "grey80", colour = NA),
    strip.text.x = theme_text(family = base_family, size = base_size * 0.8),
    strip.text.y = theme_text(family = base_family, size = base_size * 0.8,
      angle = -90),
    plot.background = theme_rect(colour = NA, fill = "white"),
    plot.title = theme_text(family = base_family, size = base_size * 1.2),
    plot.margin = unit(c(1, 1, 1, 1), "lines")), class = "options")
}

`kmg2.theme_natrisk` <-
function (base_size = 12, base_family = "") {
  structure(list(
    axis.line = theme_blank(),
    axis.text.x = theme_blank(),
    axis.text.y = theme_text(family = base_family, size = base_size * 0.8,
      lineheight = 0.9, colour = "transparent", hjust = 1),
    axis.ticks = theme_segment(colour = "transparent"),
    axis.title.x = theme_blank(),
    axis.title.y = theme_text(family = base_family, size = base_size, angle = 90,
      vjust = 0.375, colour = "transparent"),
    axis.ticks.length = unit(0.15, "cm"),
    axis.ticks.margin = unit(0.2, "cm"),
    legend.position = "none",
    panel.background = theme_rect(fill = "transparent", colour = "transparent"),
    panel.border = theme_blank(),
    panel.grid.major = theme_line(colour = "transparent"),
    panel.grid.minor = theme_line(colour = "transparent", size = 0.25),
    panel.margin = unit(0.25, "lines"),
    strip.background = theme_rect(fill = "transparent", colour = "transparent"),
    strip.text.x = theme_blank(),
    strip.text.y = theme_text(family = base_family, size = base_size * 0.8,
      angle = -90, colour = "transparent"),
    plot.background = theme_rect(colour = "transparent", fill = "transparent"),
    plot.title = theme_blank(),
    plot.margin = unit(c(0, 1, 0, 1), "lines")), class = "options")
}

`kmg2.theme_natrisk21` <-
function (base_size = 12, base_family = "") {
  structure(list(
    axis.line = theme_blank(),
    axis.text.x = theme_blank(),
    axis.text.y = theme_blank(),
    axis.ticks = theme_blank(),
    axis.title.x = theme_blank(),
    axis.title.y = theme_blank(),
    axis.ticks.length = unit(0, "cm"),
    axis.ticks.margin = unit(c(0, 0.2, 0, 0.2), "cm"),
    legend.position = "none",
    panel.background = theme_rect(fill = "transparent", colour = "transparent"),
    panel.border = theme_blank(),
    panel.grid.major = theme_blank(),
    panel.grid.minor = theme_blank(),
    panel.margin = unit(0, "lines"),
    strip.background = theme_rect(fill = "transparent", colour = "transparent"),
    strip.text.x = theme_blank(),
    strip.text.y = theme_blank(),
    plot.background = theme_rect(colour = "transparent", fill = "transparent"),
    plot.title = theme_blank(),
    plot.margin = unit(c(0, 0, 0, -2), "lines")), class = "options")
}

df <- data.frame( 
  fx = sort(rgamma(99, 12, 1), decreasing = TRUE), 
  ftrt = sample(c("Group-A", "Group-B", "Group-C"), 99, rep = TRUE),
  fcensor = sample(c(0, 1), 99, rep = TRUE, prob = c(0.1, 0.9))
)

kmg2.df <- data.frame(time=df$fx,surv=df$fcensor,trt=df$ftrt)
kmg2.fit <- survfit(Surv(time=time, event=surv, type='right') ~ trt, 
  kmg2.df)
kmg2.fit$strata <- sort(kmg2.df$trt)
kmg2.fit <- data.frame(time = kmg2.fit$time, surv = kmg2.fit$surv, strata = 
  kmg2.fit$strata,nrisk = kmg2.fit$n.risk, nevent = kmg2.fit$n.event, 
  ncensor= kmg2.fit$n.censor,stderr = kmg2.fit$std.err, upper = 
  kmg2.fit$upper, lower = kmg2.fit$lower)
kmg2.a <- unique(kmg2.fit$strata)
for (i in 1:(length(kmg2.a))) {kmg2.maxnrisk <- max(subset(kmg2.fit, 
  strata==kmg2.a[i])$nrisk);kmg2.fit[dim(kmg2.fit)[1]+1,] <- list(0, 1, 
  kmg2.a[i], kmg2.maxnrisk, 0, NA, NA, NA, NA)}
kmg2.cens <- subset(kmg2.fit, ncensor==1)
kmg2.n_majortick <- 4-1
kmg2.by <- signif((max(kmg2.fit$time)+kmg2.n_majortick/2)/kmg2.n_majortick, 
  -round(log10(max(kmg2.fit$time)), 0))
kmg2.natrisk <- by(kmg2.fit, kmg2.fit$strata,function(x, seq) {x <- 
  sort_df(x, x$surv); kmg2.natrisk <- NULL;for (i in (1:length(seq))) {for (j 
  in (1:length(x$surv))) {if (x$time[j] <= seq[i]) {kmg2.natrisk[i] <- 
  x$nrisk[j] - x$nevent[j];if (sum(x$nevent[(j:length(x$surv))]) == 0) 
  kmg2.natrisk[i] <- x$nrisk[j] - x$nevent[j] - 
  x$ncensor[j];}}};return(kmg2.natrisk)}, seq(0, kmg2.by * kmg2.n_majortick, 
  by = kmg2.by))
kmg2.m <- dim(kmg2.natrisk)
kmg2.label <- unlist(kmg2.natrisk, recursive = FALSE)
kmg2.x <- rep(seq(0, kmg2.by * kmg2.n_majortick, by = kmg2.by), kmg2.m)
kmg2.y <- NULL; kmg2.group <- NULL; kmg2.pmax <- 0.05+(kmg2.m-1)*0.05;
for (i in 1:kmg2.m) {kmg2.y <- append(kmg2.y, rep(kmg2.pmax-(i-1)*0.05, 
  kmg2.n_majortick+1));kmg2.group <- append(kmg2.group, 
  rep(names(kmg2.natrisk)[i], kmg2.n_majortick+1))}
kmg2.natrisk <- data.frame(label=kmg2.label, x=kmg2.x, y=kmg2.y, 
  group=kmg2.group)

kmg2.natrisk$y <-
  (-0.8 / (min(kmg2.natrisk$y) - max(kmg2.natrisk$y))) *
  kmg2.natrisk$y + (0.1 + 0.8 * min(kmg2.natrisk$y) /
  (min(kmg2.natrisk$y) - max(kmg2.natrisk$y)))
kmg2.natvaxis <- subset(kmg2.natrisk, x == 0)


plot1 <- ggplot(data=kmg2.fit, aes(y=surv, x=time, colour=strata)) +
  geom_step(size=1.5) + xlab("Months from entry") +
  ylab("Proportion of overall survival") + 
  geom_hline(yintercept=0, colour="white", size=1.2) + 
  geom_vline(xintercept=0, colour="white", size=1.2) + 
  geom_point(data=kmg2.cens, size=3) + scale_x_continuous(breaks = seq(0, 
  kmg2.by * kmg2.n_majortick, by = kmg2.by), limits=c(0, kmg2.by * 
  kmg2.n_majortick)) + scale_colour_brewer("Strata", 
  palette="Set1") + kmg2.theme_gray(30, "serif") + opts(legend.justification = c(1, 1))

plot2 <- ggplot(data=kmg2.natrisk, aes(y=y, x=x, label=label,
  colour=factor(group))) + geom_text() + scale_x_continuous(breaks = seq(0, 
  kmg2.by * kmg2.n_majortick, by = kmg2.by), limits=c(0, kmg2.by * kmg2.n_majortick)) +
  scale_y_continuous(limits=c(0, 1)) +
  scale_colour_brewer("Strata", palette="Set1") +
  kmg2.theme_natrisk(30, "serif")

plot3 <- ggplot(data=kmg2.natvaxis, aes(y=y, x=x, label=group,
  colour=factor(group))) + geom_text() + scale_x_continuous(limits=c(-5, 5)) +
  scale_y_continuous(limits=c(0, 1)) + scale_colour_brewer("Strata", palette="Set1") + 
  kmg2.theme_natrisk21(30, "serif")

png("myplot.png", height=7, width=7, units = "in", res = 300)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2,
  heights = unit(c(1, 5), c("null", "lines")),
  widths  = unit(c(7, 1), c("lines", "null")))))

print(plot1, vp=viewport(layout.pos.row=1, layout.pos.col=1:2))
print(plot2, vp=viewport(layout.pos.row=2, layout.pos.col=1:2))
print(plot3, vp=viewport(layout.pos.row=2, layout.pos.col=1  ))
grid.gedit("text", gp=gpar(fontfamily="serif"))
dev.off()

設定に依存しそうな気もしますが・・・。
grid package の勉強になったかも。

履歴

2011-05-21 なんかエラーが出たので直しておきました。