Creating a new stat
The simplest stat
We’ll start by creating a very simple stat: one that gives the convex
hull (the c hull) of a set of points. First we create a new
ggproto object that inherits from Stat:
StatChull <- ggproto("StatChull", Stat,
compute_group = function(data, scales) {
data[chull(data$x, data$y), , drop = FALSE]
},
required_aes = c("x", "y")
)The two most important components are the
compute_group() method (which does the computation), and
the required_aes field, which lists which aesthetics must
be present in order for the stat to work.
Next we write a layer function. Unfortunately, due to an early design
mistake I called these either stat_() or
geom_(). A better decision would have been to call them
layer_() functions: that’s a more accurate description
because every layer involves a stat and a geom.
All layer functions follow the same form - you specify defaults in
the function arguments and then call the layer() function,
sending ... into the params argument. The
arguments in ... will either be arguments for the geom (if
you’re making a stat wrapper), arguments for the stat (if you’re making
a geom wrapper), or aesthetics to be set. layer() takes
care of teasing the different parameters apart and making sure they’re
stored in the right place:
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, ...)
)
}(Note that if you’re writing this in your own package, you’ll either
need to call ggplot2::layer() explicitly, or import the
layer() function into your package namespace.)
Once we have a layer function we can try our new stat:
(We’ll see later how to change the defaults of the geom so that you
don’t need to specify fill = NA every time.)
Once we’ve written this basic object, ggplot2 gives a lot for free. For example, ggplot2 automatically preserves aesthetics that are constant within each group:
We can also override the default geom to display the convex hull in a different way:
Stat parameters
A more complex stat will do some computation. Let’s implement a
simple version of geom_smooth() that adds a line of best
fit to a plot. We create a StatLm that inherits from
Stat and a layer function, 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 is inflexible because it has no parameters. We
might want to allow the user to control the model formula and the number
of points used to generate the grid. To do so, we add arguments to the
compute_group() method and our wrapper function:
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)
Note that we don’t have to explicitly include the new
parameters in the arguments for the layer, ... will get
passed to the right place anyway. But you’ll need to document them
somewhere so the user knows about them. Here’s a brief example. Note
@inheritParams ggplot2::stat_identity: that will
automatically inherit documentation for all the parameters also defined
for 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() must be exported if you want other people to
use it. You could also consider exporting StatLm if you
want people to extend the underlying object; this should be done with
care.
Picking defaults
Sometimes you have calculations that should be performed once for the
complete dataset, not once for each group. This is useful for picking
sensible default values. For example, if we want to do a density
estimate, it’s reasonable to pick one bandwidth for the whole plot. The
following Stat creates a variation of the stat_density()
that picks one bandwidth for all groups by choosing the mean of the
“best” bandwidth for each group (I have no theoretical justification for
this, but it doesn’t seem unreasonable).
To do this we override the setup_params() method. It’s
passed the data and a list of params, and returns an updated list.
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
I recommend using NULL as a default value. If you pick
important parameters automatically, it’s a good idea to
message() to the user (and when printing a floating point
parameter, using signif() to show only a few significant
digits).
Variable names and default aesthetics
This stat illustrates another important point. If we want to make
this stat usable with other geoms, we should return a variable called
density instead of y. Then we can set up the
default_aes to automatically map density to
y, which allows the user to override it to use with
different geoms:
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")
#> Warning: `stat(density)` was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(density)` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
However, using this stat with the area geom doesn’t work quite right. The areas don’t stack on top of each other:
ggplot(mpg, aes(displ, fill = drv)) +
stat_density_common(bandwidth = 1, geom = "area", position = "stack")
This is because each density is computed independently, and the
estimated xs don’t line up. We can resolve that issue by
computing the range of the data once in 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")
Exercises
Extend
stat_chullto compute the alpha hull, as from the alphahull package. Your new stat should take analphaargument.Modify the final version of
StatDensityCommonto allow the user to specify theminandmaxparameters. You’ll need to modify both the layer function and thecompute_group()method.Note: be careful when adding parameters to a layer function. The following names col, color, pch, cex, lty, lwd, srt, adj, bg, fg, min, and max are intentionally renamed to accomodate base graphical parameter names. For example, a value passed as min to a layer appears as ymin in the
setup_paramslist of params. It is recommended you avoid using these names for layer parameters.Compare and contrast
StatLmtoggplot2::StatSmooth. What key differences makeStatSmoothmore complex thanStatLm?
