sina plot

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,

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`

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).

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,

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)) {

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

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)
+}