新しい stat の作成
最もシンプルなstat
まず、点の集合の凸包(c hull)を与える非常に単純な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")
)最も重要な2つのコンポーネントは、計算を行う compute_group() メソッドと、統計が機能するためにどの美学が存在しなければならないかをリストアップした required_aes フィールドです。
次にレイヤー関数を書きます。残念ながら、初期の設計ミスで、これらを stat_() または geom_() と呼んでいました。すべてのレイヤーには、stat_() と geom_() が含まれているため、より正確な表現となっています。
すべてのレイヤー関数は同じ形式で、関数の引数にデフォルトを指定してから、... を params 引数に送って layer() 関数を呼び出します。 ... の中の引数は、geom の引数（stat wrapper を作っている場合）、stat の引数（ geom wrapper を作っている場合）、または設定される美学のいずれかになります。 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() 関数をパッケージの名前空間にインポートする必要があることに注意してください)。
レイヤー関数ができたら、新しい stat を試してみましょう。
(ジオムのデフォルトを変更して、毎回 fill = NA を指定する必要がないようにする方法については、後で説明します)。
一旦この基本的なオブジェクトを書いてしまえば、ggplot2は多くのことを無料で提供してくれます。例えば、ggplot2は、各グループ内で一定の美観を自動的に保持します。
また、デフォルトの geom を上書きして、凸包を別の方法で表示することもできます。:
stat のパラメータ
より複雑な統計では、いくつかの計算を行います。簡単なバージョンの geom_smooth() を実装して、プロットにベストフィットの線を追加してみましょう。 ここでは、Stat を継承したStatLmと、レイヤ関数である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)
なお、レイヤの引数に新しいパラメータを明示的に含める 必要 はなく、 ... はいずれにせよ適切な場所に渡されます。 しかし、ユーザーがそれを知ることができるように、どこかに記録しておく必要があります。以下に簡単な例を示します。 注意点は @inheritParams ggplot2::stat_identity: これは stat_identity() で定義されたすべてのパラメーターのドキュメントを自動的に継承します。
#' @export
#' @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_lm()をエクスポートする必要があります。 また、基本的なオブジェクトを拡張したい場合には、 StatLm をエクスポートすることも考えられますが、これは慎重に行う必要があります。
デフォルトの選択
グループごとではなく、データセット全体に対して一度だけ実行すべき計算がある場合があります。 このような場合には、適切なデフォルト値を選択することが有効です。 例えば、密度の推定を行う場合、プロット全体で1つの帯域幅を選ぶのが妥当です。 以下の Stat は、各グループの “ベスト” なバンド幅の平均値を選択することにより、すべてのグループに1つのバンド幅を選択する stat_density() のバリエーションを作成します(これについての理論的な正当性はありませんが、不合理ではないと思われます)。
これを行うために、setup_params() メソッドをオーバーライドします。 このメソッドは、データとパラメータのリストを受け取り、更新されたリストを返します。
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
デフォルト値としては NULL を使うことをお勧めします。重要なパラメータを自動的に選ぶ場合は、ユーザーに message() を送ると良いでしょう（浮動小数点パラメータを印刷する場合は、signif() を使って有効数字数桁のみを表示します）。
変数名とデフォルトの美学
この stat は、もうひとつの重要なポイントを示しています。 この stat を他の geom で使用できるようにするには、y の代わりに density という変数を返す必要があります。そして、density を y に自動的にマッピングするように default_aes を設定し、ユーザーがそれをオーバーライドし、異なる geom で使用できるようにします。
StatDensityCommon <- ggproto("StatDensity2", Stat,
required_aes = "x",
default_aes = aes(y = stat(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 = stat(density))) +
stat_density_common(bandwidth = 1, geom = "point")
しかし、この統計を area geom で使用すると、うまくいきません。area がお互いに重ならないのです。:
ggplot(mpg, aes(displ, fill = drv)) +
stat_density_common(bandwidth = 1, geom = "area", position = "stack")
これは，各密度が独立して計算されているため，推定される x が揃わないためです。 この問題は，setup_params() でデータの範囲を一度計算することで解決できます．
StatDensityCommon <- ggproto("StatDensityCommon", Stat,
required_aes = "x",
default_aes = aes(y = stat(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 = stat(density))) +
stat_density_common(bandwidth = 1, geom = "raster")
演習問題
alphahull パッケージのように、alpha hull を計算するために
stat_chullを拡張します。 新しいstatは、alphaという引数を取る必要があります。最終バージョンの
StatDensityCommonを、ユーザがminとmaxを指定できるように変更する。 layer 関数とcompute_group()メソッドの両方を変更する必要があります。注意：レイヤー関数にパラメータを追加する際には注意が必要です。 次のような 名前 col, color, pch, cex, lty, lwd, srt, adj, bg, fg, min, max は、 ベースとなるグラフィカルなパラメータ名に合わせて、意図的に名前を変更しています。 グラフィカルなパラメータ名に合わせて意図的に改名されています。 例えば、レイヤーに min として渡された値は、
setup_paramsのリストには ymin と表示されます。 レイヤのパラメータにこれらの名前を使わないことをお勧めします。StatLmとggplot2::StatSmoothを比較対照してください。 どのような重要な違いが、StatSmoothをStatLmよりも複雑にしているでしょう？
