追記
最新バージョン(R commander 用プラグイン) が CRAN から利用できます！
RcmdrPlugin.KMggplot2_0.1-0 is on CRAN now - Triad sou.
旧バージョン
- RcmdrPlugin.KMggplot2_0.0-2.zip: Windows 用バイナリ (Dropbox)
- RcmdrPlugin.KMggplot2_0.0-2.tar.gz: Linux & Mac 用ソース (Dropbox)
いったんまとまったので、アップロードしてみました。
今月は大変忙しいので、更新できなさそう。
コメントをいただければ何か考えてみます。
まだ Kaplan-Meier plot しか描けませんが、いずれは描けるグラフを増やして CRAN に登録か・・・？
追記
R Commander で実行すると、以下のようなコードが自動生成されます。
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) plot1 <- qplot(time, surv, data=kmg2.fit, geom="step", colour = strata) plot2 <- plot1 + 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)) + geom_text(data=kmg2.natrisk, aes(y=y, x=x, label=label, colour=factor(group)), legend=FALSE) + scale_colour_brewer("Strata", palette="Set1") + kmg2.theme_gray(30, "serif") plot2
そして、今見ると何故か qplot を使っているのに気づきました。
次のバージョンでは直すようにしてみます。
最近知ったのですが Deducer という R の対話式データ解析用パッケージでは、グラフ作成に ggplot2 を使ってくれるみたいです。JGR ベースらしいです。
関数で使えると便利というコメントをいただいたので、とりあえず変数だけ指定できる関数を作ってみました。
`kmg2.theme_gray` <- function (base_size = 12, family = "") { structure(list( axis.line = theme_blank(), axis.text.x = theme_text(family = family, size = base_size * 0.8, lineheight = 0.9, colour = "grey50", vjust = 1), axis.text.y = theme_text(family = 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 = family, size = base_size, vjust = -0.5), axis.title.y = theme_text(family = family, size = base_size, angle = 90, vjust = 0.375, hjust = 0.25), 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 = family, size = base_size * 0.6), legend.title = theme_text(family = family, face = "bold", size = base_size * 0.6, hjust = 0), legend.position = c(0.95, 0.95), legend.justification = c(1, 1), 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 = family, size = base_size * 0.8), strip.text.y = theme_text(family = family, size = base_size * 0.8, angle = -90), plot.background = theme_rect(colour = NA, fill = "white"), plot.title = theme_text(family = family, size = base_size * 1.2), plot.margin = unit(c(1, 1, 1, 1), "lines")), class = "options") } `kmg.plot` <- function(time, surv, trt, natrisk=TRUE) { kmg2.df <- data.frame(time=time,surv=surv,trt=trt) 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)) if (natrisk) { 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) } plot1 <- qplot(time, surv, data=kmg2.fit, geom="step", colour = strata) plot2 <- plot1 + 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") if (natrisk) { plot2 <- plot2 + geom_text(data=kmg2.natrisk, aes(y=y, x=x, label=label, colour=factor(group)), legend=FALSE) } return(plot2) }
kmg.plot(time=df$fx, surv=df$fcensor, trt=df$ftrt)
