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

ggplot2-2.0.0 の拡張

R ggplot2

ggproto と export

ggplot2 の v2.0.0 では OO の機構が ggproto というパッケージ内に含まれた独自のものに変更され、他のパッケージからの拡張が容易になったようです。
詳しくは、Extending ggplot2 に公式の解説があるので、これを熟読すると良いでしょう。
GeomXXX, geom_xxx, StatXXX, stat_xxx などが export されたので、かなり簡単に拡張が書けるようです。

独自の geom

私も自分のパッケージで独自の geom を作ってみました。
Kaplan-Meier plot という階段状の図を作図するコンポーネントで、信頼区間も階段状に帯の様に表示したかったので、geom_ribbon を継承し、前処理を少し加えたものを作成しました。

#' Step ribbon plots.
#'
#' \code{geom_stepribbon} is an extension of the \code{geom_ribbon}, and
#' is optimized for Kaplan-Meier plots with pointwise confidence intervals
#' or a confidence band.
#'
#' @section Aesthetics:
#' \Sexpr[results=rd,stage=build]{ggplot2:::rd_aesthetics("geom", "ribbon")}
#'
#' @seealso
#'   \code{\link[ggplot2:geom_ribbon]{geom_ribbon}} \code{geom_stepribbon}
#'   inherits from \code{geom_ribbon}.
#' @inheritParams ggplot2:::geom_ribbon
#' @param kmplot If \code{TRUE}, missing values are replaced by the previous
#' values. This option is needed to make Kaplan-Meier plots if the last
#' observation has event, in which case the upper and lower values of the
#' last observation are missing. This processing is optimized for results
#' from the survfit function.
#' @examples
#' huron <- data.frame(year = 1875:1972, level = as.vector(LakeHuron))
#' h <- ggplot(huron, aes(year))
#' h + geom_stepribbon(aes(ymin = level - 1, ymax = level + 1), fill = "grey70") +
#'     geom_step(aes(y = level))
#' h + geom_ribbon(aes(ymin = level - 1, ymax = level + 1), fill = "grey70") +
#'     geom_line(aes(y = level))
#' @rdname geom_stepribbon
#' @importFrom ggplot2 layer GeomRibbon
#' @export
geom_stepribbon <- function(
  mapping = NULL, data = NULL, stat = "identity", position = "identity",
  na.rm = FALSE, show.legend = NA, inherit.aes = TRUE, kmplot = FALSE, ...) {
  layer(
    data = data,
    mapping = mapping,
    stat = stat,
    geom = GeomStepribbon,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      na.rm = na.rm,
      kmplot = kmplot,
      ...
    )
  )
}

#' @rdname geom_stepribbon
#' @format NULL
#' @usage NULL
#' @export
GeomStepribbon <- ggproto(
  "GeomStepribbon", GeomRibbon, 
  
  extra_params = c("na.rm", "kmplot"),
  
  draw_group = function(data, panel_scales, coord, na.rm = FALSE) {
    if (na.rm) data <- data[complete.cases(data[c("x", "ymin", "ymax")]), ]
    data <- rbind(data, data)
    data <- data[order(data$x), ]
    data$x <- c(data$x[2:nrow(data)], NA)
    data <- data[complete.cases(data["x"]), ]
    GeomRibbon$draw_group(data, panel_scales, coord, na.rm = FALSE)
  },
  
  setup_data = function(data, params) {
    if (params$kmplot) {
      data <- data[order(data$PANEL, data$group, data$x), ]
      tmpmin <- tmpmax <- NA
      for (i in 1:nrow(data)) {
        if (is.na(data$ymin[i])) {
          data$ymin[i] <- tmpmin
        }
        if (is.na(data$ymax[i])) {
          data$ymax[i] <- tmpmax
        }
        tmpmin <- data$ymin[i]
        tmpmax <- data$ymax[i]
      }
    }
    data
  }
)  
geom_stepribbon

継承しているため、ほぼ geom_ribbon のコピーですが、kmplot という引数を追加しています。
これは、survival パッケージの survfit という関数を内部で利用していて、この関数の出力を使う場合の例外処理の有無を指定できるようにしたかったので入れました。

GeomStepribbon

実際に処理している本体の部分です。

  • ggproto(`_class` = NULL, `_inherit` = NULL, ...) となっていて、ひとつ目の引数がクラス名、ふたつ目が継承するオブジェクトなので、GeomRibbon を継承しています、geom を一から作成する場合は Geom を継承すれば良いと思います。
  • extra_params は追加引数名の指定です。geom_xxx の方で追加引数を入れた場合はここで指定していないとエラーになります。
  • draw_group は描画処理の部分ですが、今回はここで前処理をして GeomRibbon$draw_group を呼んでいます。
  • setup_data はデータのセットアップ処理の部分です、今回はここに例外処理を入れています。

実際には、Geom 内部で順番に継承している他の処理と上述の処理が呼ばれていくのですが、今回は一部で事足りてしまいました。
今後は今回触れていない部分を勉強してみようと思っています。

出力例

出力は、geom_stepribbon または以下のようになりました。

require("ggplot2")
data(dataKm, package = "RcmdrPlugin.KMggplot2")

.df <- na.omit(data.frame(x = dataKm$time, y = dataKm$event, z = dataKm$trt))
.df <- .df[do.call(order, .df[, c("z", "x"), drop = FALSE]), , drop = FALSE]
.fit <- survival::survfit(
  survival::Surv(time = x, event = y, type = "right") ~ z, .df)
.fit <- data.frame(x = .fit$time, y = .fit$surv, nrisk = .fit$n.risk, 
  nevent = .fit$n.event, ncensor= .fit$n.censor, upper = .fit$upper,
  lower = .fit$lower)
.df <- .df[!duplicated(.df[,c("x", "z")]), ]
.df <- .fit <- data.frame(.fit, .df[, c("z"), drop = FALSE])
.df <- .fit <- rbind(unique(data.frame(x = 0, y = 1, nrisk = NA, nevent = NA,
  ncensor = NA, upper = 1, lower = 1, .df[, c("z"), drop = FALSE])), .fit)
.cens <- subset(.fit, ncensor == 1)

ggplot(data = .fit, aes(x = x, y = y, colour = z)) + 
  RcmdrPlugin.KMggplot2::geom_stepribbon(data = .fit,
    aes(x = x, ymin = lower, ymax = upper, fill = z), alpha = 0.25,
    colour = "transparent", show.legend = FALSE, kmplot = TRUE) +
  geom_step(size = 1.5) +
  geom_linerange(data = .cens, aes(x = x, ymin = y, ymax = y + 0.02),
    size = 1.5) + 
  scale_x_continuous(breaks = seq(0, 21, by = 7), limits = c(0, 21)) + 
  scale_y_continuous(limits = c(0, 1), expand = c(0.01, 0)) + 
  scale_colour_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  xlab("Time from entry") +
  ylab("Proportion of survival") +
  labs(colour = "trt") +
  theme_bw(base_size = 14, base_family = "sans") + 
  theme(legend.position = "right")

因みに、geom_ribbon では帯は折れ線状につながるので geom_line と同じ感じです。
今回の、geom_stepribbon はある意味で geom_step と合わせて使うために作っています。
あと、geom_stepribbon の引数で colour = "transparent" とすると帯の端が透明にならないため、次のバージョンアップで colour = NA にでも変えようかと思っています。

geom-.r ファイル

また、geom-.r というファイルを roxygen2 のコメントを含めて読むと、仕組みがよくわかりました。

RcmdrPlugin.KMggplot2_0.2-3 now on CRAN

R-bloggers

New version (v0.2-3) of 'RcmdrPlugin.KMggplot2' (an Rcmdr plug-in; a GUI for 'ggplot2') released.

NEWS

Changes in version 0.2-3 (2015-12-30)
  • New geom_stepribbon().
  • Pointwise confidence intervals of Kaplan-Meier plots with band (Thanks to Dr. M. Felix Freshwater).
  • Notched box plots (Thanks to Dr. M. Felix Freshwater).
  • Stratified (colored) box plots (Thanks to Dr. M. Felix Freshwater).
Changes in version 0.2-2 (2015-12-22)
  • New ggplot2's theme (theme_linedraw, theme_light).
  • New ggthemes's themes (theme_base, theme_par).
  • Fixed a bug was caused by new ggplot2's theme.
  • Fixed a bug related to ggthemes 3.0.0.
  • Fixed a bug related to windows fonts.
  • Fixed a bug related to file saving.
  • Added a vignette for dataset requirements to make Kaplan-Meier plot.
  • Added a vignette for extrafont.

geom_stepribbon()

The geom_stepribbon is an extension of the geom_ribbon, and is optimized for Kaplan-Meier plots with pointwise confidence intervals or a confidence band.

Usage
geom_stepribbon(mapping = NULL, data = NULL, stat = "identity",
  position = "identity", na.rm = FALSE, show.legend = NA,
  inherit.aes = TRUE, kmplot = FALSE, ...)

The additional argument

  • kmplot
    • If TRUE, missing values are replaced by the previous values. This option is needed to make Kaplan-Meier plots if the last observation has an event, in which case the upper and lower values of the last observation are missing. This processing is optimized for results from the survfit function.

Other arguments are the same as the geom_ribbon.

Example
require("ggplot2")
data(dataKm, package = "RcmdrPlugin.KMggplot2")

.df <- na.omit(data.frame(x = dataKm$time, y = dataKm$event, z = dataKm$trt))
.df <- .df[do.call(order, .df[, c("z", "x"), drop = FALSE]), , drop = FALSE]
.fit <- survival::survfit(
  survival::Surv(time = x, event = y, type = "right") ~ z, .df)
.fit <- data.frame(x = .fit$time, y = .fit$surv, nrisk = .fit$n.risk, 
  nevent = .fit$n.event, ncensor= .fit$n.censor, upper = .fit$upper,
  lower = .fit$lower)
.df <- .df[!duplicated(.df[,c("x", "z")]), ]
.df <- .fit <- data.frame(.fit, .df[, c("z"), drop = FALSE])
.df <- .fit <- rbind(unique(data.frame(x = 0, y = 1, nrisk = NA, nevent = NA,
  ncensor = NA, upper = 1, lower = 1, .df[, c("z"), drop = FALSE])), .fit)
.cens <- subset(.fit, ncensor == 1)

ggplot(data = .fit, aes(x = x, y = y, colour = z)) + 
  RcmdrPlugin.KMggplot2::geom_stepribbon(data = .fit,
    aes(x = x, ymin = lower, ymax = upper, fill = z), alpha = 0.25,
    colour = "transparent", show.legend = FALSE, kmplot = TRUE) +
  geom_step(size = 1.5) +
  geom_linerange(data = .cens, aes(x = x, ymin = y, ymax = y + 0.02),
    size = 1.5) + 
  scale_x_continuous(breaks = seq(0, 21, by = 7), limits = c(0, 21)) + 
  scale_y_continuous(limits = c(0, 1), expand = c(0.01, 0)) + 
  scale_colour_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  xlab("Time from entry") +
  ylab("Proportion of survival") +
  labs(colour = "trt") +
  theme_bw(base_size = 14, base_family = "sans") + 
  theme(legend.position = "right")

Pointwise confidence intervals of Kaplan-Meier plots with band


Notched box plots and stratified (colored) box plots


RcmdrPlugin.KMggplot2_0.2-1 is on CRAN now

R-bloggers

I posted a new version of the "RcmdrPlugin.KMggplot2" package, which is an Rcmdr plug-in for a "ggplot2" GUI front-end.

NEWS

Changes in version 0.2-1 (2015-12-14)
  • Kaplan-Meier plot: fixed a bug was caused by a tie data handling. (Thanks to Dr. M. Felix Freshwater).
  • Slovenian translation by Matjaž Jeran (Thanks to Matjaž Jeran and their friends).
  • New ggplot2's theme (theme_dark).
  • New ggthemes's themes (theme_calc, theme_fivethirtyeight, theme_gdocs, theme_hc, theme_pander).
  • Preview button (Thanks to Prof. Erich Neuwirth).
  • GUI refinement (Thanks to Prof. Erich Neuwirth).
  • Updates for release of the next version ggplot2.
Preview button

The new version has a preview mode.
You can make plots step-by-step.

Slovenian translation

themes
  • theme_dark (ggplot2)


  • theme_calc (ggthemes)


  • theme_fivethirtyeight (ggthemes)


  • theme_gdocs (ggthemes)


  • theme_hc (ggthemes)


  • theme_pander (ggthemes)

ggplot2の最新版(1.1.0?)ではまった点のメモ

R ggplot2

そろそろ ggplot2 の最新版が CRAN にリリースされるようなので、自分のパッケージのメンテナンスをしました。
メンテナンス中に困った点をメモしておこうと思います。

stat_summary

fun.args でパラメータを渡さないとエラー。

# <= 1.0.1
stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", conf.int = .95, width = 0.1)

# > 1.0.1
stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", width = 0.1, fun.args = list(conf.int = 0.95))

stat_boxplot

stat_paramsgeom_params を使うとエラーになる。
最新版ではパラメータが直接渡せるようになり、少し便利に。

# <= 1.0.1
stat_boxplot(geom = "errorbar", stat_params = list(width = 0.5), geom_params = list())

# > 1.0.1
stat_boxplot(geom = "errorbar", width = 0.5)

show_guide

show_guide を使うと警告 (deprecated)。
最新版では show.legend を使う様になった。

# <= 1.0.1
geom_xxx(..., show_guide = FALSE, ...)

# > 1.0.1
geom_xxx(..., show.legend = FALSE, ...)

grid::unit

grid::unitggplot2 上にインポートされるようになったため、呼び出しが楽になった。
grid::arrowscales::alpha も同様にインポートされるようになった。

# <= 1.0.1
grid::unit(...)

# > 1.0.1
unit(...)

SAS/STAT 14.1 が出たようだ

SAS

Enhancements in SAS/STAT 14.1 Software


以下の変更点が特に気になった。

  • LIFETEST Procedureが Gray's test + CIF のノンパラメトリック推定量 (Gray 1988) に対応 (かなりうれしい!)
  • MIXED Procedure でも ddfm = KENWARDROGER2 が追加された、Kenward & Roger (2009) の改良版自由度推定量が利用できるようになった
  • POWER Procedure が Cox 回帰のサンプルサイズ設計に対応 (生存時間の非劣性にはまだ対応しないらしい、残念)
  • MI Procedure の多重補完法のデフォルト補完回数が 5 から 25 へ変更された(もっと多くてもよさそう)。やっぱ 5 は少なすぎる・・・
  • GEE Procedure の正式版リリース
  • GLIMMIX Procedure の積分計算方法の改善 (FASTQUAD option; Pinheiro and Chao, 2006)
  • HPSPLIT Procedure が追加され、SAS でも CART が使える様になった、個人的にグラフが綺麗だと思う
  • MCMC Procedure で ODE への対応 (PKモデルなど)、NLMIXED 同様の GENERAL 関数に対する積分への対応 (周辺モデルなど)、すごいな。

文献

  • Gray RJ. A class of K-sample tests for comparing the cumulative incidence of a competing risk. Annals of Statistics 1988; 16: 1141–1154.
  • Kenward MG, Roger JH. An improved approximation to the precision of fixed effects from restricted maximum likelihood. Computational Statistics and Data Analysis 2009; 53: 2583–2595.
  • Pinheiro JC, Chao EC. Efficient Laplacian and adaptive Gaussian quadrature algorithms for multilevel generalized linear mixed models. Journal of Computational and Graphical Statistics 2006; 15: 58–81.

SAS/STAT 13.2

SAS

またバージョンが上がったいたらしい!
かなり使い勝手がよくなりそうな機能拡張が含まれているので、結構期待しています(特にNLMIXED)。

GEE Procedure (評価版)

ついにGEEのProcedureが出たらしい。
Observation-specific and subject-specific weighted estimating equationsを実装とのこと。
Liang and Zeger (1986)のGEEとあるので、GEE1だけっぽい?

追記

CAUSALTRT procedure (評価版?)

GEEに加え (Weighted Methods for Analyzing Missing Data with the GEE and CAUSALTRT Procedures) doubly robust estimation approachについて実装したCAUSALTRT procedureも追加されるらしい。
昨今注目が集まっているためか、missing data関係の解析ソフトウェアが整備されていくのはうれしい。

ICPHREG Procedure

13.1のICLIFETESTに続き、ICPHREGも出た。
区間打ち切りの生存時間解析もやりやすくなりましたね。

機能拡張

  • NLMIXED procedureが複数のRANDOM statementをサポート。これで、階層構造を持った変量効果のモデリングが簡単にできるようになるはず。
  • FREQ procedureがオッズ比と相対リスクに対するスコア信頼区間をサポート。
  • MCMC procedureがカテゴリカルな確率分布をサポート(MODEL, RANDOM, PRIOR statment)。


その他にも色々更新されていました (SAS/STAT(R) 13.2 User's Guide)。
無料版のSAS University Editionではまだ使えないようでした (PROC PRODUCT_STATUS; run;でバージョン確認できます)。