新しいstatをつくる
もっとも単純なstat
とても単純なstatをつくることから始めます。点の集合の凸包を得るstatです。まず、Statを継承するggprotoオブジェクトをつくります。
StatChull <- ggproto("StatChull", Stat,
compute_group = function(data, scales) {
data[chull(data$x, data$y), , drop = FALSE]
},
required_aes = c("x", "y")
)もっとも重要な部分は、compute_group()メソッド（計算を実行します）とrequired_aesフィールド（このstatが動くのに存在しなければいけないaestheticsをリストアップします）です。
次に、layer関数を書きます。不幸にも、早期の設計ミスにより私はこれらをstat_()やgeom_()と名付けてしまいました。layer_()関数とするべきでした。これはより実体に即した表記です。なぜなら、すべてのレイヤーがstatにもgeomにも関わっているからです。
すべてのlayer関数は同じ形式に従います。関数の引数にデフォルト引数を設定して、layer()を呼び、...をそのparam引数に渡します。 ...に含まれる引数は、geomの引数（statのラッパーの場合）かstatの引数（geomのラッパーの場合）かセットするaestheticsのいずれかです。layer()が、異なるパラメータを分離して正しい場所に格納することの面倒を見てくれます。
stat_chull <- function(mapping = NULL, data = NULL, geom = "polygon",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, ...) {
layer(
stat = StatChull, data = data, mapping = mapping, geom = geom,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(na.rm = na.rm, ...)
)
}（注：もし自分のパッケージの中でこれを書くなら、ggplot2::layer()と明示的に呼ぶか、layer()関数をパッケージの名前空間にインポートするか、いずれかを行う必要があります）
いちどlayer関数ができれば、新しいstatを試してみましょう。
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
stat_chull(fill = NA, colour = "black")
（fill=NAと毎回指定しなくてもいいようにgeomのデフォルトを変更する方法は後ほど出てきます）
一度この基本的なオブジェクトを書き上げると、ggplot2が提供するたくさんの機能を自由に使えるようになります。例えば、ggplot2はそれぞれのグループで固定のaestheticsを自動で保存してくれます。
ggplot(mpg, aes(displ, hwy, colour = drv)) +
geom_point() +
stat_chull(fill = NA)
凸包を別の方法で表示するようにデフォルトのgeomをオーバーライドすることもできます。
ggplot(mpg, aes(displ, hwy)) +
stat_chull(geom = "point", size = 4, colour = "red") +
geom_point()
Statのパラメータ
さらに複雑なstatはいくつかの計算を行います。プロットにもっともフィットする線を追加する、簡単なバージョンのgeom_smooth()を実装してみましょう。Statを継承したStatLmとlayer関数stat_lm()をつくります。
StatLm <- ggproto("StatLm", Stat,
required_aes = c("x", "y"),
compute_group = function(data, scales) {
rng <- range(data$x, na.rm = TRUE)
grid <- data.frame(x = rng)
mod <- lm(y ~ x, data = data)
grid$y <- predict(mod, newdata = grid)
grid
}
)
stat_lm <- function(mapping = NULL, data = NULL, geom = "line",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, ...) {
layer(
stat = StatLm, data = data, mapping = mapping, geom = geom,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(na.rm = na.rm, ...)
)
}
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
stat_lm()
StatLmはパラメータを持たないため柔軟性がありません。モデルの式やグリッドを生成するのに使われる点の数をユーザーが操作できるようにしたい、と思うかもしれません。そうするには、compute_group()メソッドと先ほどのラッパー関数に引数を追加します。
StatLm <- ggproto("StatLm", Stat,
required_aes = c("x", "y"),
compute_group = function(data, scales, params, n = 100, formula = y ~ x) {
rng <- range(data$x, na.rm = TRUE)
grid <- data.frame(x = seq(rng[1], rng[2], length = n))
mod <- lm(formula, data = data)
grid$y <- predict(mod, newdata = grid)
grid
}
)
stat_lm <- function(mapping = NULL, data = NULL, geom = "line",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, n = 50, formula = y ~ x,
...) {
layer(
stat = StatLm, data = data, mapping = mapping, geom = geom,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(n = n, formula = formula, na.rm = na.rm, ...)
)
}
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
stat_lm(formula = y ~ poly(x, 10)) +
stat_lm(formula = y ~ poly(x, 10), geom = "point", colour = "red", n = 20)
注意すべき点として、新しいパラメータを明示的にlayerの引数に含めることは、必ずやらなくてはいけないわけではありません。いずれにせよ...は正しい位置に渡されます。しかし、ユーザーがそれについて知ることができるようにどこかにドキュメントを書く必要はあるでしょう。 ここに短い例があります。@inheritParams ggplot2::stat_identityは、stat_identity()について定義されたすべてのパラメータについてのドキュメントを自動的に継承してくれます。
#' @inheritParams ggplot2::stat_identity
#' @param formula The modelling formula passed to \code{lm}. Should only
#' involve \code{y} and \code{x}
#' @param n Number of points used for interpolation.
stat_lm <- function(mapping = NULL, data = NULL, geom = "line",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, n = 50, formula = y ~ x,
...) {
layer(
stat = StatLm, data = data, mapping = mapping, geom = geom,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(n = n, formula = formula, na.rm = na.rm, ...)
)
}デフォルト値を使う
場合によっては、それぞれのデータに対して一度づつではなく、データセット全体に行うべき計算があります。 これは、理にかなったデフォルト値を選ぶのに役立ちます。例えば、密度推定を行いたいとき、プロット全体に対してひとつのバンド幅を選ぶのは妥当なことでしょう。 以下のStatは、それぞれのグループに最適なバンド幅の平均を全グループのバンド幅として使う（この方法に論理的な根拠はありませんが、非合理的ではなさそうに見えます）、stat_density()の一種をつくります。
そうするためには、setup_param()メソッドをオーバーライドします。これは、データとパラメータのリストを渡すと更新されたリストを返します。
StatDensityCommon <- ggproto("StatDensityCommon", Stat,
required_aes = "x",
setup_params = function(data, params) {
if (!is.null(params$bandwidth))
return(params)
xs <- split(data$x, data$group)
bws <- vapply(xs, bw.nrd0, numeric(1))
bw <- mean(bws)
message("Picking bandwidth of ", signif(bw, 3))
params$bandwidth <- bw
params
},
compute_group = function(data, scales, bandwidth = 1) {
d <- density(data$x, bw = bandwidth)
data.frame(x = d$x, y = d$y)
}
)
stat_density_common <- function(mapping = NULL, data = NULL, geom = "line",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, bandwidth = NULL,
...) {
layer(
stat = StatDensityCommon, data = data, mapping = mapping, geom = geom,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(bandwidth = bandwidth, na.rm = na.rm, ...)
)
}
ggplot(mpg, aes(displ, colour = drv)) +
stat_density_common()
#> Picking bandwidth of 0.345
ggplot(mpg, aes(displ, colour = drv)) +
stat_density_common(bandwidth = 0.5)
私はNULLをデフォルト値にすることをお薦めします。もしも重要なパラメータを自動で選ぶのであれば、message()でユーザーに伝えるのがいいでしょう。(あと、浮動小数点のパラメータを表示するときはsignif()を使って数桁だけを見せるのがいいでしょう)
変数名とデフォルトのaesthetics
以下のstatは別の重要な点を示しています。もしも他のgeomといっしょにこのstatを使いたい場合は、yではなくdensityという名前の変数を返すべきです。 そうすることで、densityをyに自動的にマップするようにdefault_aesを設定することができます。ユーザーはこのaestheticsをオーバーライドして別のgeomで使うことができます。
StatDensityCommon <- ggproto("StatDensity2", Stat,
required_aes = "x",
default_aes = aes(y = ..density..),
compute_group = function(data, scales, bandwidth = 1) {
d <- density(data$x, bw = bandwidth)
data.frame(x = d$x, density = d$y)
}
)
ggplot(mpg, aes(displ, drv, colour = ..density..)) +
stat_density_common(bandwidth = 1, geom = "point")
しかし、このstatをareaのgeomと組み合わせると正しく動作しません。面が積み重なっていません。
ggplot(mpg, aes(displ, fill = drv)) +
stat_density_common(bandwidth = 1, geom = "area", position = "stack")
これはなぜかというと、各グループの密度が別々に計算されていてxの推定値が含まれていないからです。この問題は、setup_param()で一度だけデータの範囲を計算するようにすることで解決できます。
StatDensityCommon <- ggproto("StatDensityCommon", Stat,
required_aes = "x",
default_aes = aes(y = ..density..),
setup_params = function(data, params) {
min <- min(data$x) - 3 * params$bandwidth
max <- max(data$x) + 3 * params$bandwidth
list(
bandwidth = params$bandwidth,
min = min,
max = max,
na.rm = params$na.rm
)
},
compute_group = function(data, scales, min, max, bandwidth = 1) {
d <- density(data$x, bw = bandwidth, from = min, to = max)
data.frame(x = d$x, density = d$y)
}
)
ggplot(mpg, aes(displ, fill = drv)) +
stat_density_common(bandwidth = 1, geom = "area", position = "stack")
ggplot(mpg, aes(displ, drv, fill = ..density..)) +
stat_density_common(bandwidth = 1, geom = "raster")
Exercises
stat_chullを拡張して、alphahullがやっているようにアルファシェイプを計算するようにしなさい。 新しいstatはalphaを引数に取ること。StatDensityCommonの最終バージョンに変更を加え、ユーザーがminとmaxパラメーターを指定できるようにしなさい。layer関数とcompute_group()メソッドの両方を変更する必要がある。StatLmとggplot2::StatSmoothを比較しなさい。StatSmoothをStatLmより複雑にしている重要な違いは何？
