diff --git a/NAMESPACE b/NAMESPACE index 441effcd..5e7f1164 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -50,6 +50,7 @@ export(type_ribbon) export(type_ridge) export(type_rug) export(type_segments) +export(type_sina) export(type_spineplot) export(type_spline) export(type_summary) @@ -148,6 +149,7 @@ importFrom(stats, qt, quantile, reformulate, + runif, sd, setNames, spline, diff --git a/NEWS.md b/NEWS.md index f189be98..d6711d0b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -32,6 +32,10 @@ where the formatting is also better._ along the chosen margin by default. It also reverses the y-axis by default, so that the first row sits at the top (again matching `heatmap()`); pass an explicit `ylim` to override. (#677 @grantmcdermott) +- `type_sina()` / `"sina"` for [sina plots](https://en.wikipedia.org/wiki/Sina_plot), + a variant of the violin plot where the raw observations are displayed as points, + with each group's width bounded by its density. This makes them arguably a more + principled version of the beeswarm plot. (#734 @grantmcdermott) - While not strictly a new plot type, `type_area()` gains a new `stack` argument for drawing _stacked_ area plots, where each layer represents a discrete `by` category group. This functionality is further enhanced by two (also new) @@ -216,6 +220,21 @@ related to plot layering. See "Bug fixes" below. - `type_text()` no longer converts a categorical axis to a numeric one. (#730 @grantmcdermott) +- The `adjust` argument of `type_density()`, `type_violin()`, and + `type_ridge()` was accepted but never passed on to the underlying + `density()` call, so it silently did nothing. (#734 @grantmcdermott) +- `type_violin()` receives two further bug fixes: + - `joint.bw = "full"` computed the joint bandwidth from the `x` categories + rather than the `y` values being smoothed. (#734 @grantmcdermott) + - Dodge offsets were keyed by the `by` variable's underlying integer codes. + A numeric `by` therefore errored outright, and a factor `by` carrying an + unused level silently dropped the affected group from the plot. Offsets are + now keyed by position among the observed groups. (#734 @grantmcdermott) + - A numeric (continuous) `by` now reverts to discrete groups and a matching + discrete legend, as it already did for `"boxplot"`, `"polygon"` and the + other types that cannot render a colour gradient. Previously every violin + was drawn in the same colour while the legend showed a colourbar. + (#734 @grantmcdermott) - Fixed a bug where consecutive plots with (i) logged axes under (ii) a dynamic theme would error, due to a stale `par("xlog")`/`par("ylog")` state. We now avoid this by grabbing the log state directly from the top-level `log` diff --git a/R/by_aesthetics.R b/R/by_aesthetics.R index 0c6c3b98..8315fe2b 100755 --- a/R/by_aesthetics.R +++ b/R/by_aesthetics.R @@ -17,7 +17,7 @@ by_aesthetics = function(settings) { by_continuous = !null_by && inherits(datapoints$by, c("numeric", "integer")) # The connected line types go through segmented_lines() instead. "b" is # still excluded pending its gap handling. - if (isTRUE(by_continuous) && type %in% c("b", "ribbon", "polygon", "polypath", "boxplot", "chull")) { + if (isTRUE(by_continuous) && type %in% c("b", "ribbon", "polygon", "polypath", "boxplot", "chull", "violin")) { # Only warn if a legend would actually be drawn: the reversion to a discrete # legend still needs to happen for correct grouping, but it's not worth # flagging when the user has suppressed the legend anyway (#656). diff --git a/R/sanitize_type.R b/R/sanitize_type.R index d75f7db9..bd3bea02 100644 --- a/R/sanitize_type.R +++ b/R/sanitize_type.R @@ -46,6 +46,7 @@ sanitize_type = function(settings) { "ridge", "rug", "segments", + "sina", "spine", "spineplot", "spline", "summary", @@ -111,6 +112,7 @@ sanitize_type = function(settings) { "ridge" = type_ridge, "rug" = type_rug, "segments" = type_segments, + "sina" = type_sina, "spine" = type_spineplot, "spineplot" = type_spineplot, "spline" = type_spline, diff --git a/R/type_density.R b/R/type_density.R index e2cbf24c..e3e9a8b8 100644 --- a/R/type_density.R +++ b/R/type_density.R @@ -211,7 +211,7 @@ data_density = function(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, dens = lapply(datapoints, function(dat) { wts = if (has_weights) dat[["weights"]] / sum(dat[["weights"]]) else NULL - density(dat$x, bw = dens_bw, kernel = kernel, n = n, weights = wts) + density(dat$x, bw = dens_bw, adjust = adjust, kernel = kernel, n = n, weights = wts) }) if (length(echo.bw)) { diff --git a/R/type_ridge.R b/R/type_ridge.R index d97dee89..0330c958 100644 --- a/R/type_ridge.R +++ b/R/type_ridge.R @@ -359,7 +359,7 @@ data_ridge = function(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, } datapoints = lapply(datapoints, function(dat) { - dens = density(dat$x, bw = dens_bw, kernel = kernel, n = n) + dens = density(dat$x, bw = dens_bw, adjust = adjust, kernel = kernel, n = n) out = data.frame( by = dat$by[1], # already split facet = dat$facet[1], # already split diff --git a/R/type_sina.R b/R/type_sina.R new file mode 100644 index 00000000..abbb4d46 --- /dev/null +++ b/R/type_sina.R @@ -0,0 +1,162 @@ +#' @rdname type_violin +#' @param method character string giving how `type_sina()` spreads points +#' across the available width at each `y` value. `"quasirandom"` (the +#' default) walks a low-discrepancy sequence, which fills the width more +#' evenly than chance does and---unlike [`type_jitter`]---is deterministic, +#' so repeated calls give the same plot without setting a seed. `"random"` +#' draws the displacements uniformly at random instead. +#' @importFrom stats approx density runif weighted.mean +#' @order 2 +#' @export +type_sina = function( + bw = "nrd0", + joint.bw = c("mean", "full", "none"), + adjust = 1, + kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine"), + n = 512, + trim = FALSE, + width = 0.9, + method = c("quasirandom", "random"), + singletons = c("keep", "warn", "drop") + ) { + kernel = match.arg(kernel, c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine")) + method = match.arg(method, c("quasirandom", "random")) + singletons = match.arg(singletons, c("keep", "warn", "drop")) + if (is.logical(joint.bw)) { + joint.bw = ifelse(joint.bw, "mean", "none") + } + joint.bw = match.arg(joint.bw, c("mean", "full", "none")) + out = list( + data = data_sina(bw = bw, adjust = adjust, kernel = kernel, n = n, + joint.bw = joint.bw, trim = trim, width = width, + method = method, singletons = singletons), + draw = draw_points(), + # points, as far as the rest of the package is concerned: this gates + # `pch`, `cex` and the bubble/legend handling. See by_pch(). + name = "p" + ) + class(out) = "tinyplot_type" + return(out) +} + + +## The first `n` van der Corput values, mapped onto [-1, 1] and centred: a +## deterministic stand-in for runif(-1, 1) that spreads successive values evenly +## rather than independently, so points fill the available width instead of +## clumping and leaving gaps. +## +## The raw sequence only balances about the midpoint at n = 2^k - 1; at every +## other n it leans left (n = 2 gives c(0, -0.5), so neither point sits right of +## the tick). Centring fixes the lean, and can overshoot the envelope by ~0.1% +## in the process, so clamp it back -- points escaping the density bounds would +## defeat the whole point of the type. +centred_van_der_corput = function(n) { + u = 2 * van_der_corput(n) - 1 + u = u - mean(u) + pmax(pmin(u, 1), -1) +} + + +van_der_corput = function(n, base = 2) { + vapply( + seq_len(n), + function(i) { + out = 0 + f = 1 / base + while (i > 0) { + out = out + f * (i %% base) + i = i %/% base + f = f / base + } + out + }, + numeric(1) + ) +} + + +data_sina = function(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, + joint.bw = "none", trim = FALSE, width = 0.9, + method = "quasirandom", singletons = "keep") { + fun = function(settings, ...) { + env2env(settings, environment(), c("datapoints", "by", "facet", "ylab", "col", "bg", "null_by", "null_facet")) + + specials = dist_specials(datapoints, null_by, null_facet) + + prep = dist_prep( + datapoints, null_by = null_by, null_facet = null_facet, + specials = specials, bw = bw, joint.bw = joint.bw, width = width, + singletons = singletons, gradient = TRUE + ) + cells = prep[["cells"]] + dens_bw = prep[["dens_bw"]] + xwidth = prep[["xwidth"]] + group_offsets = prep[["group_offsets"]] + grp_levels = prep[["grp_levels"]] + offsets_axis = prep[["offsets_axis"]] + xlabs = prep[["xlabs"]] + + datapoints = lapply(cells, function(dat) { + nobs = nrow(dat) + xcat = dat[["x"]][1] + dodge = if (prep[["dodged"]]) { + group_offsets[match(dat[["by"]][1], grp_levels)] + } else { + 0 + } + + # a lone observation has no density to be displaced by, so it just + # sits on its group's tick + if (nobs < 2L) { + dat[["x"]] = xcat + dodge + if (facet_drop_levels_on(settings[["facet.args"]])) dat[[".xcat"]] = xcat + return(dat) + } + + if (trim) { + yrng = range(dat[["y"]]) + dens = density(dat[["y"]], bw = dens_bw, adjust = adjust, + kernel = kernel, n = n, from = yrng[1], to = yrng[2]) + } else { + dens = density(dat[["y"]], bw = dens_bw, adjust = adjust, + kernel = kernel, n = n) + } + + # the violin's half-width at each observation, normalized to [0, 1] + halfwidth = approx(dens[["x"]], dens[["y"]], xout = dat[["y"]], rule = 2)[["y"]] + halfwidth = halfwidth / max(dens[["y"]]) + + if (method == "random") { + u = runif(nobs, -1, 1) + } else { + # walk the sequence in y order, so that neighbouring points + # (the ones at risk of overlapping) land far apart + u = numeric(nobs) + u[order(dat[["y"]])] = centred_van_der_corput(nobs) + } + + dat[["x"]] = xcat + u * halfwidth * xwidth / 2 + dodge + # the displacement moves points off their own tick; carry the + # category position. See cat_axis_codes() + if (facet_drop_levels_on(settings[["facet.args"]])) dat[[".xcat"]] = xcat + return(dat) + }) + datapoints = do.call(rbind, datapoints) + + by = if (length(unique(datapoints[["by"]])) == 1) by else datapoints[["by"]] + facet = if (length(unique(datapoints[["facet"]])) == 1) facet else datapoints[["facet"]] + + env2env(environment(), settings, c( + "datapoints", + "by", + "facet", + "ylab", + "xlabs", + "col", + "bg", + "group_offsets", + "offsets_axis" + )) + } + return(fun) +} diff --git a/R/type_violin.R b/R/type_violin.R index 4f407628..bb2d5646 100644 --- a/R/type_violin.R +++ b/R/type_violin.R @@ -1,31 +1,68 @@ -#' Violin plot type +#' Violin and sina plot types #' #' @md -#' @description Type function for violin plots, which are an alternative to box +#' @description Type functions for violin plots, which are an alternative to box #' plots for visualizing continuous distributions (by group) in the form of -#' mirrored densities. +#' mirrored densities. `type_violin()` draws the smooth outline of each +#' density, while `type_sina()` scatters the underlying observations as points +#' within the same outline (similar to a beeswarm plot). #' @inheritParams type_density -#' @param trim logical indicating whether the violins should be trimmed to the -#' range of the data. Default is `FALSE`. +#' @param trim logical indicating whether the densities should be trimmed to +#' the range of the data. Default is `FALSE`. For `type_sina()` this only +#' affects the envelope that bounds the displacement, so it matters mainly +#' for exact alignment with a `type_violin()` layer drawn the same way. #' @param width numeric (ideally in the range `[0, 1]`, although this isn't -#' enforced) giving the normalized width of the individual violins. +#' enforced) giving the normalized width of the individual violins. For +#' `type_sina()` this is the width of the (undrawn) violin that the points +#' are scattered inside. #' @param lighten logical. Should the fills use a lighter, opaque tint of the #' series colour(s)? Default is `TRUE`, which keeps single- and multi-group #' displays consistent and lets the fill read cleanly over grid lines. Set to -#' `FALSE` to use the fully-saturated palette colour(s) instead. +#' `FALSE` to use the fully-saturated palette colour(s) instead. Only +#' applies to `type_violin()`, since `type_sina()` has no fill of its own. #' @param singletons character string indicating what to do with singleton #' groups, i.e. combinations of `x`, `by`, and `facet` that consist of only 1 -#' row. The default `"warn"` option removes any singleton cases and emits a -#' warning reporting how many there were. `"drop"` does the same thing, but -#' quietly. In either case the dropped groups may still be represented as -#' empty violins or facets in your plot. Finally, `"none"` skips all singleton +#' row. Both types accept `"warn"`, which removes any singleton cases and +#' emits a warning reporting how many there were, and `"drop"`, which does +#' the same thing quietly. In either case the dropped groups may still be +#' represented as empty violins or facets in your plot. +#' +#' The remaining option differs by type, as does the default. `type_violin()` +#' defaults to `"warn"` and also accepts `"none"`, which skips all singleton #' checks and retains the affected groups; possibly leading to an error. Note -#' that singletons require a numeric `bw`, since the data-driven bandwidth -#' rules need at least 2 observations. +#' that singletons then require a numeric `bw`, since the data-driven +#' bandwidth rules need at least 2 observations. +#' +#' `type_sina()` instead defaults to `"keep"`, which draws the lone +#' observation on its group's tick. No density can be estimated from a single +#' point, but the point itself is still worth showing, and discarding an +#' observation from what is fundamentally a scatter plot is worse than +#' discarding an unrenderable violin. There is no `"none"` for this type, +#' since `"keep"` already retains these cases without error. #' @inherit stats::density details #' @details See [`type_density`] for more details and considerations related to #' bandwidth selection and kernel types. #' +#' A sina plot (Sidiropoulos et al., 2018) is closely related to a beeswarm +#' plot, but is arguably the more principled of the two. Both spread a +#' group's observations sideways to expose its shape. A beeswarm does so by +#' packing points until they no longer collide, which makes its width an +#' artefact of the *rendering*: change the symbol size or the device and the +#' swarm changes shape. A sina instead displaces each point by the kernel +#' density at its own `y` value, so its width is a property of the *data*, +#' stable across devices and comparable between groups. The trade-off is +#' occlusion: a beeswarm guarantees that no point hides another, whereas a +#' sina accepts the occasional overlap. If you need collision-free packing, +#' use a dedicated package such as \CRANpkg{beeswarm}. +#' @references +#' Sidiropoulos, N., Sohi, S. H., Pedersen, T. L., Porse, B. T., Winther, O., +#' Rapin, N., and Bagger, F. O. (2018). \cite{SinaPlot: An Enhanced Chart for +#' Simple and Truthful Representation of Single Observations Over Multiple +#' Classes}. Journal of Computational and Graphical Statistics, 27(3), 673-676. +#' Available: https://doi.org/10.1080/10618600.2017.1366914 +#' @seealso [type_boxplot], [type_density] and [type_ridge] for the other ways +#' of displaying a distribution by group, and [type_jitter] for displacing +#' points without reference to a density. #' #' @examples #' # "violin" type convenience string @@ -52,10 +89,26 @@ #' # dodged grouped violin plot example (different dataset) #' tinyplot(len ~ dose | supp, data = ToothGrowth, type = "violin") #' -#' # note: above we relied on `...` argument passing alongside the "violin" -#' # type convenience string. But this won't work for `width`, since it will +#' # the "sina" type shows the observations themselves, rather than a smooth +#' # outline drawn around them +#' tinyplot(weight ~ feed, data = chickwts, type = "sina") +#' +#' # layering a sina on top of a violin lines up exactly, since the points are +#' # displaced by the violin's own half-width +#' tinyplot(weight ~ feed, data = chickwts, type = "violin") +#' tinyplot_add(type = "sina", pch = 16, col = "black") +#' +#' # unlike `type_violin()`, `type_sina()` supports a continuous `by` variable; +#' # it colours the points rather than splitting them into groups +#' tinyplot( +#' Sepal.Length ~ Species | Petal.Width, data = iris, +#' type = "sina", pch = 16 +#' ) +#' +#' # note: above we relied on `...` argument passing alongside the type +#' # convenience strings. But this won't work for `width`, since it will #' # clash with the top-level `tinyplot(..., width = )` arg. To ensure -#' # correct arg passing, it's safer to use the functional `type_violin()` type. +#' # correct arg passing, it's safer to use the functional type. #' tinyplot( #' len ~ dose | supp, data = ToothGrowth, #' type = type_violin(width = 0.75) @@ -63,6 +116,7 @@ #' #' @importFrom stats density weighted.mean #' @importFrom stats bw.SJ bw.bcv bw.nrd bw.nrd0 bw.ucv +#' @order 1 #' @export type_violin = function( bw = "nrd0", @@ -102,47 +156,15 @@ data_violin = function(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, env2env(settings, environment(), c("datapoints", "by", "null_palette", "facet", "ylab", "col", "bg", "log", "null_by", "null_facet")) settings[["lighten"]] = lighten - - # Handle ordering based on by and facet variables - ngrps = if (null_by) 1 else length(unique(datapoints$by)) - nfacets = if (null_facet) 1 else length(unique(datapoints$facet)) - - # catch for special cases - x_by = y_by = facet_by = FALSE - if (!null_by) { - x_by = identical(datapoints$x, datapoints$by) - y_by = identical(datapoints$y, datapoints$by) - if (!null_facet) facet_by = identical(datapoints$facet, datapoints$by) - } - + specials = dist_specials(datapoints, null_by, null_facet) + # FIXME (once we add support for gradient fill to draw_polygon) - if (y_by) { + if (specials[["y_by"]]) { warning("\n`y` == `by` is not currently supported for `type_violin`. We hope to support this in a future release, but for now `y` grouping will be turned off automatically.\n") by = NULL datapoints$by = "" - ngrps = 1 null_by = TRUE - } - - # Convert x to factor if it's not already - datapoints$x = as.factor(datapoints$x) - if (x_by) datapoints$by = datapoints$x - - # Handle factor levels and maintain order - xlvls = levels(datapoints$x) - xlabs = seq_along(xlvls) - names(xlabs) = xlvls - # xlabs = levels(datapoints$x) - datapoints$x = as.integer(datapoints$x) - - if (null_by && null_facet) { - xord = order(datapoints$x) - } else if (null_facet) { - xord = order(datapoints$by, datapoints$x) - } else if (null_by) { - xord = order(datapoints$facet, datapoints$x) - } else { - xord = order(datapoints$by, datapoints$facet, datapoints$x) + specials = dist_specials(datapoints, null_by, null_facet) } if (length(unique(datapoints[["by"]])) == 1 && null_palette) { @@ -153,85 +175,52 @@ data_violin = function(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, bg = "by" } - # Reorder x, y, ymin, and ymax based on the order determined - datapoints = datapoints[xord,] - - - datapoints = split(datapoints, list(datapoints$x, datapoints$by, datapoints$facet)) - datapoints = drop_singletons(datapoints, singletons) - - if (joint.bw == "none" || is.numeric(bw)) { - dens_bw = bw - } else { - if (joint.bw == "mean") { - # Use weighted mean of subgroup bandwidths - bws = sapply(datapoints, function(dat) bw_fun(kernel = bw, dat$y)) - ws = sapply(datapoints, nrow) - dens_bw = weighted.mean(bws, ws) - } else if (joint.bw == "full") { - dens_bw = bw_fun(kernel = bw, unlist(sapply(datapoints, `[[`, "x"))) - } - } - - # Compute group offsets for multi-group violins - if (ngrps > 1 && isFALSE(x_by) && isFALSE(facet_by)) { - xwidth_grp = width / ngrps - 0.01 - group_offsets = seq( - -((width - xwidth_grp) / 2), - ((width - xwidth_grp) / 2), - length.out = ngrps - ) - } else { - group_offsets = rep(0, max(ngrps, 1)) - } - offsets_axis = "x" + prep = dist_prep( + datapoints, null_by = null_by, null_facet = null_facet, + specials = specials, bw = bw, joint.bw = joint.bw, width = width, + singletons = singletons + ) + cells = prep[["cells"]] + dens_bw = prep[["dens_bw"]] + xwidth = prep[["xwidth"]] + group_offsets = prep[["group_offsets"]] + grp_levels = prep[["grp_levels"]] + offsets_axis = prep[["offsets_axis"]] + xlabs = prep[["xlabs"]] - datapoints = lapply(seq_along(datapoints), function(d) { - dat = datapoints[[d]] + datapoints = lapply(cells, function(dat) { if (trim) { yrng = range(dat$y) - dens = density(dat$y, bw = dens_bw, kernel = kernel, n = n, from = yrng[1], to = yrng[2]) + dens = density(dat$y, bw = dens_bw, adjust = adjust, kernel = kernel, n = n, from = yrng[1], to = yrng[2]) } else { - dens = density(dat$y, bw = dens_bw, kernel = kernel, n = n) + dens = density(dat$y, bw = dens_bw, adjust = adjust, kernel = kernel, n = n) } - + x = dens$y y = dens$x - - + if (log %in% c("y", "xy")) { if (x[1] <= 0) { warning("\nNon-positive density values have been trimmed as part of the logarthmic transformation.\n") xidx = x > 0 x = x[xidx] y = y[xidx] - } + } } - + + # mirror the density about the category position x = c(x, rev(-x)) y = c(y, rev(y)) - - xwidth = xwidth_orig = width - # dodge groups (if any) - if ((ngrps > 1) && isFALSE(x_by) && isFALSE(facet_by)) { - xwidth = xwidth_orig / ngrps - 0.01 - xcat = as.numeric(sub("^([0-9]+)\\..*", "\\1", names(datapoints)[d])) - x = rescale_num(x, to = c(0, xwidth)) - x = x + xcat - xwidth/2 - x = x + group_offsets[dat$by[1]] - } else if (nfacets > 1) { - xcat = as.numeric(sub("^([0-9]+)\\..*", "\\1", names(datapoints)[d])) - x = rescale_num(x, to = c(0, xwidth)) - x = x + xcat - xwidth/2 - } else { - xcat = d - x = rescale_num(x, to = c(0, xwidth)) - x = x + xcat - xwidth/2 + + xcat = dat$x[1] + x = rescale_num(x, to = c(0, xwidth)) + xcat - xwidth / 2 + if (prep[["dodged"]]) { + x = x + group_offsets[match(dat$by[1], grp_levels)] } - + x = c(x, NA) y = c(y, NA) - + out = data.frame( by = dat$by[1], # already split facet = dat$facet[1], # already split @@ -248,13 +237,13 @@ data_violin = function(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, by = if (length(unique(datapoints$by)) == 1) by else datapoints$by facet = if (length(unique(datapoints$facet)) == 1) facet else datapoints$facet - + # legend customizations settings$legend_args[["pch"]] = settings$legend_args[["pch"]] %||% 22 settings$legend_args[["pt.cex"]] = settings$legend_args[["pt.cex"]] %||% 3.5 settings$legend_args[["y.intersp"]] = settings$legend_args[["y.intersp"]] %||% 1.25 settings$legend_args[["seg.len"]] = settings$legend_args[["seg.len"]] %||% 1.25 - + env2env(environment(), settings, c( "datapoints", "by", @@ -269,3 +258,160 @@ data_violin = function(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, } return(fun) } + + +## Shared setup for the categorical distribution types, below. +## +## `type_violin()` and `type_sina()` both plot a categorical `x` against a +## numeric `y`, estimate one density per (x, by, facet) cell, and lay the +## result out at integer x positions with the `by` groups dodged around each +## tick. Everything up to the density is common; only what each type does +## *with* it differs -- violin traces the outline, sina scatters the +## observations inside it. That common part lives here, since `type_violin()` +## is also the home of the help page the two types share. `type_sina()` lives +## in type_sina.R and calls into these. + + +## Detect the "special" cases where `by` duplicates another aesthetic. Callers +## need these before any of the setup below, since e.g. `type_violin()` has to +## resolve a `y`-valued `by` (which it cannot render) before the split. +dist_specials = function(datapoints, null_by, null_facet) { + out = list(x_by = FALSE, y_by = FALSE, facet_by = FALSE) + if (isTRUE(null_by)) return(out) + out[["x_by"]] = identical(datapoints[["x"]], datapoints[["by"]]) + out[["y_by"]] = identical(datapoints[["y"]], datapoints[["by"]]) + if (isFALSE(null_facet)) { + out[["facet_by"]] = identical(datapoints[["facet"]], datapoints[["by"]]) + } + return(out) +} + + +## Coerce `x` to consecutive integer positions, split the data into density +## cells, pick a bandwidth, and work out the dodge offsets. +## +## `gradient` says the caller can render a continuous `by` (sina colours its +## points by it). Splitting on such a `by` would put nearly every observation +## in a cell of its own, so those cells are keyed on (x, facet) alone and the +## groups are left undodged. +## +## Returns the pieces the caller needs; type-specific concerns (fills, legend +## keys) stay with the caller. +dist_prep = function( + datapoints, + null_by, + null_facet, + specials, + bw = "nrd0", + joint.bw = "none", + width = 0.9, + singletons = "warn", + gradient = FALSE + ) { + x_by = specials[["x_by"]] + facet_by = specials[["facet_by"]] + + by_continuous = isTRUE(gradient) && isFALSE(null_by) && + inherits(datapoints[["by"]], c("numeric", "integer")) + + ngrps = if (null_by || by_continuous) 1L else length(unique(datapoints[["by"]])) + nfacets = if (null_facet) 1L else length(unique(datapoints[["facet"]])) + + ## The groups that the dodge slots below are built for, in slot order. + ## Callers must look a row's slot up by *position among these*, never by the + ## factor's integer code: a `by` carrying an unused level, or a numeric one + ## whose values are not positions at all, would otherwise index past the end + ## of the offsets and silently turn those rows' `x` into NA. + grp_levels = if (null_by || by_continuous) { + NULL + } else if (is.factor(datapoints[["by"]])) { + levels(droplevels(datapoints[["by"]])) + } else { + sort(unique(datapoints[["by"]])) + } + + ## Convert x to consecutive integer positions, keeping the labels for the axis + datapoints[["x"]] = as.factor(datapoints[["x"]]) + if (x_by) datapoints[["by"]] = datapoints[["x"]] + xlvls = levels(datapoints[["x"]]) + xlabs = seq_along(xlvls) + names(xlabs) = xlvls + datapoints[["x"]] = as.integer(datapoints[["x"]]) + + if (null_by && null_facet) { + xord = order(datapoints[["x"]]) + } else if (null_facet) { + xord = order(datapoints[["by"]], datapoints[["x"]]) + } else if (null_by) { + xord = order(datapoints[["facet"]], datapoints[["x"]]) + } else { + xord = order(datapoints[["by"]], datapoints[["facet"]], datapoints[["x"]]) + } + datapoints = datapoints[xord, ] + + ## A continuous `by` is a colour scale rather than a grouping, so it must + ## not key the split; every cell would hold a single observation. + if (by_continuous) { + cells = split(datapoints, list(datapoints[["x"]], datapoints[["facet"]])) + } else { + cells = split( + datapoints, + list(datapoints[["x"]], datapoints[["by"]], datapoints[["facet"]]) + ) + } + ## "keep" leaves the 1-row cells in place for the caller to draw as-is; it + ## still has to shed the 0-row cells that split() invents, which is what + ## "none" does. + cells = drop_singletons(cells, if (singletons == "keep") "none" else singletons) + + ## Bandwidth rules need at least 2 observations, so a singleton cell kept + ## under "keep" must not feed into them (the caller skips its density + ## anyway). Under "none" it deliberately does, since that mode's contract + ## is to run the checks the user opted out of and let them fail. + smoothable = if (singletons == "keep") { + cells[vapply(cells, nrow, integer(1)) > 1L] + } else { + cells + } + ## With every cell a retained singleton there is nothing to derive a + ## bandwidth from -- and nothing to smooth either, since the caller draws + ## those rows directly -- so leave `bw` as supplied rather than erroring. + if (joint.bw == "none" || is.numeric(bw) || length(smoothable) == 0L) { + dens_bw = bw + } else if (joint.bw == "mean") { + # Use weighted mean of subgroup bandwidths + bws = sapply(smoothable, function(dat) bw_fun(kernel = bw, dat[["y"]])) + ws = sapply(smoothable, nrow) + dens_bw = weighted.mean(bws, ws) + } else { + dens_bw = bw_fun(kernel = bw, unlist(sapply(smoothable, `[[`, "y"))) + } + + ## Dodge: groups share a tick, so each gets a narrower slot beside it + dodged = ngrps > 1 && isFALSE(x_by) && isFALSE(facet_by) + if (dodged) { + xwidth = width / ngrps - 0.01 + group_offsets = seq( + -((width - xwidth) / 2), + ((width - xwidth) / 2), + length.out = ngrps + ) + } else { + xwidth = width + group_offsets = rep(0, max(ngrps, 1)) + } + + list( + cells = cells, + xlabs = xlabs, + dens_bw = dens_bw, + group_offsets = group_offsets, + grp_levels = grp_levels, + offsets_axis = "x", + ngrps = ngrps, + nfacets = nfacets, + dodged = dodged, + xwidth = xwidth, + by_continuous = by_continuous + ) +} diff --git a/R/zzz.R b/R/zzz.R index 92244d07..53bcf30e 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -106,6 +106,7 @@ "xmax_dep", "xmin", "xmin_dep", + "xpad", "y", "y_dep", "yaxb", @@ -120,6 +121,7 @@ "ymax", "ymax_dep", "ymin", - "ymin_dep" + "ymin_dep", + "ypad" )) } diff --git a/altdoc/quarto_website.yml b/altdoc/quarto_website.yml index ccae1755..24e36155 100644 --- a/altdoc/quarto_website.yml +++ b/altdoc/quarto_website.yml @@ -134,6 +134,8 @@ website: file: man/type_rug.qmd - text: type_segments file: man/type_segments.qmd + - text: type_sina + file: man/type_violin.qmd - text: type_spineplot file: man/type_spineplot.qmd - text: type_spline diff --git a/inst/tinytest/_tinysnapshot/sina_by_continuous.svg b/inst/tinytest/_tinysnapshot/sina_by_continuous.svg new file mode 100644 index 00000000..8309bd38 --- /dev/null +++ b/inst/tinytest/_tinysnapshot/sina_by_continuous.svg @@ -0,0 +1,241 @@ + + + + + + + + + + + + + + + 0.5 + 1.0 + 1.5 + 2.0 + 2.5 + + + + + + + + + + +Petal.Width + + + + + + + +Species +Sepal.Length + + + + + + +setosa +versicolor +virginica + + + + + + + + + +4.5 +5.0 +5.5 +6.0 +6.5 +7.0 +7.5 +8.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inst/tinytest/_tinysnapshot/sina_groups.svg b/inst/tinytest/_tinysnapshot/sina_groups.svg new file mode 100644 index 00000000..fa070093 --- /dev/null +++ b/inst/tinytest/_tinysnapshot/sina_groups.svg @@ -0,0 +1,137 @@ + + + + + + + + + + + + + + + +supp +OJ +VC + + + + + + + +dose +len + + + + + + +0.5 +1 +2 + + + + + + + + +5 +10 +15 +20 +25 +30 +35 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inst/tinytest/_tinysnapshot/sina_on_violin.svg b/inst/tinytest/_tinysnapshot/sina_on_violin.svg new file mode 100644 index 00000000..cf4753c8 --- /dev/null +++ b/inst/tinytest/_tinysnapshot/sina_on_violin.svg @@ -0,0 +1,142 @@ + + + + + + + + + + + + + +feed +weight + + + + + + + +casein +horsebean +linseed +meatmeal +soybean +sunflower + + + + + + +100 +200 +300 +400 +500 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inst/tinytest/_tinysnapshot/sina_singletons_keep.svg b/inst/tinytest/_tinysnapshot/sina_singletons_keep.svg new file mode 100644 index 00000000..1102c332 --- /dev/null +++ b/inst/tinytest/_tinysnapshot/sina_singletons_keep.svg @@ -0,0 +1,136 @@ + + + + + + + + + + + + + +cyl +mpg + + + + + + + + + + + + + +4 +6 +8 + + + + + + +10 +15 +20 +25 +30 + +0 + + + + + + + + + + + + + + +4 +6 +8 + + + + + + +10 +15 +20 +25 +30 + +1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inst/tinytest/test-type_sina.R b/inst/tinytest/test-type_sina.R new file mode 100644 index 00000000..2fb8cf9f --- /dev/null +++ b/inst/tinytest/test-type_sina.R @@ -0,0 +1,42 @@ +source("helpers.R") +using("tinysnapshot") + +# core geometry, including dodged groups +f = function() { + plt(len ~ dose | supp, data = ToothGrowth, type = "sina", pch = 16) +} +expect_snapshot_plot(f, label = "sina_groups") + +# points are displaced by the violin's own half-width, so a sina layer fits +# inside a violin exactly +f = function() { + plt(weight ~ feed, data = chickwts, type = "violin") + plt_add(type = "sina", pch = 16, col = "black") +} +expect_snapshot_plot(f, label = "sina_on_violin") + +# unlike type_violin(), a continuous `by` is supported: it colours the points +# (with a gradient legend) rather than keying the density split +f = function() { + plt(Sepal.Length ~ Species | Petal.Width, data = iris, type = "sina", pch = 16) +} +expect_snapshot_plot(f, label = "sina_by_continuous") + +# cyl == 4 & vs == 0 is a single car. No density can be estimated from it, but +# unlike type_violin() the default draws the observation (on its tick) rather +# than discarding it. +f = function() { + plt(mpg ~ cyl, facet = ~vs, data = mtcars, type = "sina", pch = 16) +} +expect_snapshot_plot(f, label = "sina_singletons_keep") + +# the default "quasirandom" method never touches the RNG, so repeated calls are +# identical without a seed (c.f. `type_jitter()`); "random" does draw +set.seed(42) +seed_before = .Random.seed +plt(count ~ spray, data = InsectSprays, type = "sina") +expect_identical(.Random.seed, seed_before) +plt(count ~ spray, data = InsectSprays, type = type_sina(method = "random")) +expect_false(identical(.Random.seed, seed_before)) + +expect_error(type_sina(method = "nope")) diff --git a/inst/tinytest/test-type_violin.R b/inst/tinytest/test-type_violin.R index de32b106..80f75f8e 100644 --- a/inst/tinytest/test-type_violin.R +++ b/inst/tinytest/test-type_violin.R @@ -44,6 +44,22 @@ f = function() { } expect_snapshot_plot(f, label = "violin_groups_argpass") +# dodge offsets used to be keyed by the `by` variable's integer codes. A +# numeric `by` errored outright; it now downgrades to discrete groups like the +# other polygon-alike types, so the fills match the legend. (#734) +expect_warning( + plt(len ~ supp | dose, data = ToothGrowth, type = "violin"), + pattern = "Continuous legends not supported" +) + +# a factor `by` carrying an unused level silently dropped that group (#734) +expect_silent( + plt( + len ~ supp | factor(dose, levels = c(0.25, 0.5, 1, 2)), + data = ToothGrowth, type = "violin" + ) +) + # don't dodge if by (groups) and x are the same f = function() { plt(Sepal.Length ~ Species | Species, iris, type = "violin", legend = FALSE) diff --git a/man/type_violin.Rd b/man/type_violin.Rd index c42fca4d..6cacb924 100644 --- a/man/type_violin.Rd +++ b/man/type_violin.Rd @@ -1,8 +1,9 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/type_violin.R +% Please edit documentation in R/type_violin.R, R/type_sina.R \name{type_violin} \alias{type_violin} -\title{Violin plot type} +\alias{type_sina} +\title{Violin and sina plot types} \usage{ type_violin( bw = "nrd0", @@ -16,6 +17,19 @@ type_violin( lighten = TRUE, singletons = c("warn", "drop", "none") ) + +type_sina( + bw = "nrd0", + joint.bw = c("mean", "full", "none"), + adjust = 1, + kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", + "cosine", "optcosine"), + n = 512, + trim = FALSE, + width = 0.9, + method = c("quasirandom", "random"), + singletons = c("keep", "warn", "drop") +) } \arguments{ \item{bw}{the smoothing \code{\link[stats:bw.nrd]{bandwidth}} to be used, @@ -48,35 +62,71 @@ the version used by S.} always makes sense to specify \code{n} as a power of two. } -\item{trim}{logical indicating whether the violins should be trimmed to the -range of the data. Default is \code{FALSE}.} +\item{trim}{logical indicating whether the densities should be trimmed to +the range of the data. Default is \code{FALSE}. For \code{type_sina()} this only +affects the envelope that bounds the displacement, so it matters mainly +for exact alignment with a \code{type_violin()} layer drawn the same way.} \item{width}{numeric (ideally in the range \verb{[0, 1]}, although this isn't -enforced) giving the normalized width of the individual violins.} +enforced) giving the normalized width of the individual violins. For +\code{type_sina()} this is the width of the (undrawn) violin that the points +are scattered inside.} \item{lighten}{logical. Should the fills use a lighter, opaque tint of the series colour(s)? Default is \code{TRUE}, which keeps single- and multi-group displays consistent and lets the fill read cleanly over grid lines. Set to -\code{FALSE} to use the fully-saturated palette colour(s) instead.} +\code{FALSE} to use the fully-saturated palette colour(s) instead. Only +applies to \code{type_violin()}, since \code{type_sina()} has no fill of its own.} \item{singletons}{character string indicating what to do with singleton groups, i.e. combinations of \code{x}, \code{by}, and \code{facet} that consist of only 1 -row. The default \code{"warn"} option removes any singleton cases and emits a -warning reporting how many there were. \code{"drop"} does the same thing, but -quietly. In either case the dropped groups may still be represented as -empty violins or facets in your plot. Finally, \code{"none"} skips all singleton +row. Both types accept \code{"warn"}, which removes any singleton cases and +emits a warning reporting how many there were, and \code{"drop"}, which does +the same thing quietly. In either case the dropped groups may still be +represented as empty violins or facets in your plot. + +The remaining option differs by type, as does the default. \code{type_violin()} +defaults to \code{"warn"} and also accepts \code{"none"}, which skips all singleton checks and retains the affected groups; possibly leading to an error. Note -that singletons require a numeric \code{bw}, since the data-driven bandwidth -rules need at least 2 observations.} +that singletons then require a numeric \code{bw}, since the data-driven +bandwidth rules need at least 2 observations. + +\code{type_sina()} instead defaults to \code{"keep"}, which draws the lone +observation on its group's tick. No density can be estimated from a single +point, but the point itself is still worth showing, and discarding an +observation from what is fundamentally a scatter plot is worse than +discarding an unrenderable violin. There is no \code{"none"} for this type, +since \code{"keep"} already retains these cases without error.} + +\item{method}{character string giving how \code{type_sina()} spreads points +across the available width at each \code{y} value. \code{"quasirandom"} (the +default) walks a low-discrepancy sequence, which fills the width more +evenly than chance does and---unlike \code{\link{type_jitter}}---is deterministic, +so repeated calls give the same plot without setting a seed. \code{"random"} +draws the displacements uniformly at random instead.} } \description{ -Type function for violin plots, which are an alternative to box +Type functions for violin plots, which are an alternative to box plots for visualizing continuous distributions (by group) in the form of -mirrored densities. +mirrored densities. \code{type_violin()} draws the smooth outline of each +density, while \code{type_sina()} scatters the underlying observations as points +within the same outline (similar to a beeswarm plot). } \details{ See \code{\link{type_density}} for more details and considerations related to bandwidth selection and kernel types. + +A sina plot (Sidiropoulos et al., 2018) is closely related to a beeswarm +plot, but is arguably the more principled of the two. Both spread a +group's observations sideways to expose its shape. A beeswarm does so by +packing points until they no longer collide, which makes its width an +artefact of the \emph{rendering}: change the symbol size or the device and the +swarm changes shape. A sina instead displaces each point by the kernel +density at its own \code{y} value, so its width is a property of the \emph{data}, +stable across devices and comparable between groups. The trade-off is +occlusion: a beeswarm guarantees that no point hides another, whereas a +sina accepts the occasional overlap. If you need collision-free packing, +use a dedicated package such as \CRANpkg{beeswarm}. } \examples{ # "violin" type convenience string @@ -103,13 +153,41 @@ tinyplot(weight ~ feed | feed, data = chickwts, type = "violin", legend = FALSE) # dodged grouped violin plot example (different dataset) tinyplot(len ~ dose | supp, data = ToothGrowth, type = "violin") -# note: above we relied on `...` argument passing alongside the "violin" -# type convenience string. But this won't work for `width`, since it will +# the "sina" type shows the observations themselves, rather than a smooth +# outline drawn around them +tinyplot(weight ~ feed, data = chickwts, type = "sina") + +# layering a sina on top of a violin lines up exactly, since the points are +# displaced by the violin's own half-width +tinyplot(weight ~ feed, data = chickwts, type = "violin") +tinyplot_add(type = "sina", pch = 16, col = "black") + +# unlike `type_violin()`, `type_sina()` supports a continuous `by` variable; +# it colours the points rather than splitting them into groups +tinyplot( + Sepal.Length ~ Species | Petal.Width, data = iris, + type = "sina", pch = 16 +) + +# note: above we relied on `...` argument passing alongside the type +# convenience strings. But this won't work for `width`, since it will # clash with the top-level `tinyplot(..., width = )` arg. To ensure -# correct arg passing, it's safer to use the functional `type_violin()` type. +# correct arg passing, it's safer to use the functional type. tinyplot( len ~ dose | supp, data = ToothGrowth, type = type_violin(width = 0.75) ) } +\references{ +Sidiropoulos, N., Sohi, S. H., Pedersen, T. L., Porse, B. T., Winther, O., +Rapin, N., and Bagger, F. O. (2018). \cite{SinaPlot: An Enhanced Chart for +Simple and Truthful Representation of Single Observations Over Multiple +Classes}. Journal of Computational and Graphical Statistics, 27(3), 673-676. +Available: https://doi.org/10.1080/10618600.2017.1366914 +} +\seealso{ +\link{type_boxplot}, \link{type_density} and \link{type_ridge} for the other ways +of displaying a distribution by group, and \link{type_jitter} for displacing +points without reference to a density. +} diff --git a/vignettes/types.qmd b/vignettes/types.qmd index e417d8e3..a729445c 100644 --- a/vignettes/types.qmd +++ b/vignettes/types.qmd @@ -89,6 +89,7 @@ a convenience string (with default behaviour) or a companion `type_*()` function | `"histogram"` / `"hist"` | `type_histogram()` | Creates a histogram of a single variable. | [link](/man/type_histogram.qmd) | | `"qq"` | `type_qq()` | Creates a quantile-quantile plot. | [link](/man/type_qq.qmd) | | `"ridge"` | `type_ridge()` | Creates a ridgeline (aka joy) plot. | [link](/man/type_ridge.qmd) | +| `"sina"` | `type_sina()` | Scatters observations within a violin envelope. | [link](/man/type_violin.qmd) | | `"spineplot"` / `"spine"` | `type_spineplot()` | Creates a spine plot or spinogram. | [link](/man/type_spineplot.qmd) | | `"violin"` | `type_violin()` | Creates a violin plot. | [link](/man/type_violin.qmd) |