Change density default bw = "nrd0" to bw = "sj"

NEWS.md

@@ -1,6 +1,9 @@
 
 # ggplot2 (development version)
 
+* (breaking) the default `bw` argument in `stat_density()` and `stat_ydensity()`
+  has changed from `"nrd0"` to `"sj"` in keeping with the recommendation given
+  in `?density` (@teunbrand, #3825). 
 * When facets coerce the faceting variables to factors, the 'ordered' class
   is dropped (@teunbrand, #5666).
 * `coord_map()` and `coord_polar()` throw informative warnings when used

R/stat-density.R

@@ -33,7 +33,7 @@
 stat_density <- function(mapping = NULL, data = NULL,
                          geom = "area", position = "stack",
                          ...,
-                         bw = "nrd0",
+                         bw = "sj",
                          adjust = 1,
                          kernel = "gaussian",
                          n = 512,
@@ -91,7 +91,7 @@ StatDensity <- ggproto("StatDensity", Stat,
 
   extra_params = c("na.rm", "orientation"),
 
-  compute_group = function(data, scales, bw = "nrd0", adjust = 1, kernel = "gaussian",
+  compute_group = function(data, scales, bw = "sj", adjust = 1, kernel = "gaussian",
                            n = 512, trim = FALSE, na.rm = FALSE, bounds = c(-Inf, Inf),
                            flipped_aes = FALSE) {
     data <- flip_data(data, flipped_aes)
@@ -110,7 +110,7 @@ StatDensity <- ggproto("StatDensity", Stat,
 
 )
 
-compute_density <- function(x, w, from, to, bw = "nrd0", adjust = 1,
+compute_density <- function(x, w, from, to, bw = "sj", adjust = 1,
                             kernel = "gaussian", n = 512,
                             bounds = c(-Inf, Inf)) {
   nx <- length(x)
@@ -220,8 +220,16 @@ reflect_density <- function(dens, bounds, from, to) {
 # Similar to stats::density.default
 # Once R4.3.0 is the lowest supported version, this function can be replaced by
 # using `density(..., warnWbw = FALSE)`.
-precompute_bw = function(x, bw = "nrd0") {
+precompute_bw = function(x, bw = "sj") {
   bw <- bw[1]
+  if (length(x) < 2) {
+    cli::cli_abort("{.arg x} must contain at least 2 elements to select a \\
+                   bandwidth automatically.")
+  }
+  if (zero_range(range(x))) {
+    # Many other bandwidth methods do not handle 0-variance input
+    return(stats::bw.nrd0(x))
+  }
   if (is.character(bw)) {
     bw <- arg_match0(bw, c("nrd0", "nrd", "ucv", "bcv", "sj", "sj-ste", "sj-dpi"))
     bw <- switch(

R/stat-ydensity.R

@@ -26,7 +26,7 @@
 stat_ydensity <- function(mapping = NULL, data = NULL,
                           geom = "violin", position = "dodge",
                           ...,
-                          bw = "nrd0",
+                          bw = "sj",
                           adjust = 1,
                           kernel = "gaussian",
                           trim = TRUE,
@@ -78,7 +78,7 @@ StatYdensity <- ggproto("StatYdensity", Stat,
 
   extra_params = c("na.rm", "orientation"),
 
-  compute_group = function(self, data, scales, width = NULL, bw = "nrd0", adjust = 1,
+  compute_group = function(self, data, scales, width = NULL, bw = "sj", adjust = 1,
                        kernel = "gaussian", trim = TRUE, na.rm = FALSE,
                        drop = TRUE, flipped_aes = FALSE, bounds = c(-Inf, Inf)) {
     if (nrow(data) < 2) {
@@ -96,7 +96,7 @@ StatYdensity <- ggproto("StatYdensity", Stat,
     }
     range <- range(data$y, na.rm = TRUE)
     modifier <- if (trim) 0 else 3
-    bw <- calc_bw(data$y, bw)
+    bw <- precompute_bw(data$y, bw)
     dens <- compute_density(
       data$y, data[["weight"]],
       from = range[1] - modifier * bw, to = range[2] + modifier * bw,
@@ -118,7 +118,7 @@ StatYdensity <- ggproto("StatYdensity", Stat,
     dens
   },
 
-  compute_panel = function(self, data, scales, width = NULL, bw = "nrd0", adjust = 1,
+  compute_panel = function(self, data, scales, width = NULL, bw = "sj", adjust = 1,
                            kernel = "gaussian", trim = TRUE, na.rm = FALSE,
                            scale = "area", flipped_aes = FALSE, drop = TRUE,
                            bounds = c(-Inf, Inf)) {
@@ -151,24 +151,3 @@ StatYdensity <- ggproto("StatYdensity", Stat,
 
   dropped_aes = "weight"
 )
-
-calc_bw <- function(x, bw) {
-  if (is.character(bw)) {
-    if (length(x) < 2) {
-      cli::cli_abort("{.arg x} must contain at least 2 elements to select a bandwidth automatically.")
-    }
-
-    bw <- switch(
-      to_lower_ascii(bw),
-      nrd0 = stats::bw.nrd0(x),
-      nrd = stats::bw.nrd(x),
-      ucv = stats::bw.ucv(x),
-      bcv = stats::bw.bcv(x),
-      sj = ,
-      `sj-ste` = stats::bw.SJ(x, method = "ste"),
-      `sj-dpi` = stats::bw.SJ(x, method = "dpi"),
-      cli::cli_abort("{.var {bw}} is not a valid bandwidth rule.")
-    )
-  }
-  bw
-}

tests/testthat/_snaps/stat-ydensity.md

@@ -4,5 +4,5 @@
 
 ---
 
-    `test` is not a valid bandwidth rule.
+    `bw` must be one of "nrd0", "nrd", "ucv", "bcv", "sj", "sj-ste", or "sj-dpi", not "test".
 

tests/testthat/test-geom-violin.R

@@ -8,8 +8,8 @@ test_that("range is expanded", {
     geom_violin(trim = FALSE) +
     facet_grid(x ~ ., scales = "free") +
     coord_cartesian(expand = FALSE)
-  expand_a <- stats::bw.nrd0(df$y[df$x == "a"]) * 3
-  expand_b <- stats::bw.nrd0(df$y[df$x == "b"]) * 3
+  expand_a <- stats::bw.SJ(df$y[df$x == "a"]) * 3
+  expand_b <- stats::bw.SJ(df$y[df$x == "b"]) * 3
   expect_equal(layer_scales(p, 1)$y$dimension(), c(0 - expand_a, 1 + expand_a))
   expect_equal(layer_scales(p, 2)$y$dimension(), c(0 - expand_b, 2 + expand_b))
 })

tests/testthat/test-stat-density.R

@@ -1,7 +1,7 @@
 test_that("stat_density actually computes density", {
   # Compare functon approximations because outputs from `ggplot()` and
   # `density()` give grids spanning different ranges
-  dens <- stats::density(mtcars$mpg)
+  dens <- stats::density(mtcars$mpg, bw = "sj")
   expected_density_fun <- stats::approxfun(data.frame(x = dens$x, y = dens$y))
 
   plot <- ggplot(mtcars, aes(mpg)) + stat_density()
@@ -19,7 +19,7 @@ test_that("stat_density can make weighted density estimation", {
   df <- mtcars
   df$weight <- mtcars$cyl / sum(mtcars$cyl)
 
-  dens <- stats::density(df$mpg, weights = df$weight, bw = bw.nrd0(df$mpg))
+  dens <- stats::density(df$mpg, weights = df$weight, bw = bw.SJ(df$mpg))
   expected_density_fun <- stats::approxfun(data.frame(x = dens$x, y = dens$y))
 
   plot <- ggplot(df, aes(mpg, weight = weight)) + stat_density()
@@ -38,7 +38,7 @@ test_that("stat_density uses `bounds`", {
   mpg_max <- max(mtcars$mpg)
 
   expect_bounds <- function(bounds) {
-    dens <- stats::density(mtcars$mpg)
+    dens <- stats::density(mtcars$mpg, bw = "sj")
     orig_density <- stats::approxfun(
       data.frame(x = dens$x, y = dens$y),
       yleft = 0,

tests/testthat/test-stat-ydensity.R

@@ -1,7 +1,7 @@
 test_that("calc_bw() requires at least two values and correct method", {
-  expect_snapshot_error(calc_bw(1, "nrd0"))
-  expect_silent(calc_bw(1:5, "nrd0"))
-  expect_snapshot_error(calc_bw(1:5, "test"))
+  expect_snapshot_error(precompute_bw(1, "nrd0"))
+  expect_silent(precompute_bw(1:5, "nrd0"))
+  expect_snapshot_error(precompute_bw(1:5, "test"))
 })
 
 test_that("`drop = FALSE` preserves groups with 1 observations", {