diff --git a/DESCRIPTION b/DESCRIPTION
index e80f0eaa..7cb6b1fd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: plot2
 Type: Package
 Title: Lightweight extension of base R plot
-Version: 0.0.3.9003
+Version: 0.0.3.9004
 Authors@R: 
   c(
     person(
diff --git a/NEWS.md b/NEWS.md
index 6f8ee341..fbd7390e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,4 @@
-# plot2 0.0.3.9003 (development version)
+# plot2 0.0.3.9004 (development version)
 
 New features:
 
@@ -7,6 +7,8 @@ convenience keyword similar to `lty` and `pch`. This is useful for plotting
 filled point characters (e.g., pch = 21), where you want a different colour for
 the fill and border. (#50 @grantmcdermott)
 
+- Support for filled density plots. (#58 @grantmcdermott)
+
 Bug fixes:
 
 - Y-label correctly prints if a function was used for the atomic plot method,
diff --git a/R/by_aesthetics.R b/R/by_aesthetics.R
index d1018ac0..ac264a7b 100755
--- a/R/by_aesthetics.R
+++ b/R/by_aesthetics.R
@@ -3,7 +3,7 @@ by_col = function(ngrps = 1L, col = NULL, palette = NULL) {
   # palette = substitute(palette, env = parent.env(environment()))
   
   # special "by" convenience keyword (will treat as NULL & handle grouping below)
-  if (!is.null(col) && length(col)==1 && col=="by") col = NULL
+  if (!anyNA(col) && !is.null(col) && length(col)==1 && col=="by") col = NULL
 
   if (is.null(col) && is.null(palette)) {
     col = seq_len(ngrps)
@@ -15,7 +15,7 @@ by_col = function(ngrps = 1L, col = NULL, palette = NULL) {
     } else if (length(col) != ngrps) {
       stop(sprintf("`col` must be of length 1 or %s.", ngrps), call. = FALSE)
     }
-    if (is.character(col)) {
+    if (anyNA(col) || is.character(col)) {
       return(col)
     }
   }
diff --git a/R/plot2.R b/R/plot2.R
index d616a7a8..86088936 100644
--- a/R/plot2.R
+++ b/R/plot2.R
@@ -30,7 +30,9 @@
 #'   and "h" for histogram-like vertical lines. "n" does not produce
 #'   any points or lines.
 #'   - Additional plot2 types: "pointrange", "errorbar", and "ribbon" for
-#'   drawing these respective plot types.
+#'   drawing these interval plot types. (Note that specifying "ribbon" for
+#'   objects of class `density` will yield a density plot with a shaded
+#'   interior.)
 #' @param xlim the x limits (x1, x2) of the plot. Note that x1 > x2 is allowed
 #'   and leads to a ‘reversed axis’. The default value, NULL, indicates that
 #'   the range of the `finite` values to be plotted should be used.
@@ -127,7 +129,12 @@
 #'   be calling `dev.off()` to reset all `par` settings to their defaults.)
 #' @param subset,na.action,drop.unused.levels arguments passed to `model.frame`
 #'   when extracting the data from `formula` and `data`.
-#' @param ymin,ymax minimum and maximum coordinates of the point range. Only used when `type="pointrange"`.
+#' @param ymin,ymax minimum and maximum coordinates of interval plot types. Only
+#'   used when the `type` argument is one of "pointrange", "errorbar", or
+#'   "ribbon".
+#' @param ribbon_alpha numeric factor modifying the opacity alpha of any ribbon
+#'   shading; typically in `[0, 1]`. Default value is 0.2. Only used when
+#'   `type = "ribbon"`.
 #' @param ... 	other `graphical` parameters (see `par` and also the "Details"
 #'   section of `plot`).
 #'   
@@ -277,6 +284,7 @@ plot2.default = function(
     par_restore = FALSE,
     ymin = NULL,
     ymax = NULL,
+    ribbon_alpha = 0.2,
     ...) {
   
   dots = list(...)
@@ -373,6 +381,13 @@ plot2.default = function(
   } else {
     bg = rep(bg, ngrps)
   }
+  if (type == "ribbon") {
+    if (!is.null(bg)) {
+      bg = adjustcolor(bg, ribbon_alpha)
+    } else if (!is.null(col)) {
+      bg = adjustcolor(col, ribbon_alpha)
+    }
+  }
   
   # Save current graphical parameters
   opar = par(no.readonly = TRUE)
@@ -436,6 +451,13 @@ plot2.default = function(
   ) {
     legend.args[["pt.cex"]] = cex
   }
+  if (type=="ribbon") {
+    if (is.null(legend.args[["pch"]])) legend.args[["pch"]] = 22
+    if (is.null(legend.args[["pt.cex"]])) legend.args[["pt.cex"]] = 3.5
+    if (is.null(legend.args[["pt.lwd"]])) legend.args[["pt.lwd"]] = 0
+    if (is.null(legend.args[["y.intersp"]])) legend.args[["y.intersp"]] = 1.25
+    if (is.null(legend.args[["seg.len"]])) legend.args[["seg.len"]] = 1.25
+  }
   
   if (legend.args[["x"]] != "none") {
     
@@ -519,6 +541,15 @@ plot2.default = function(
       
       legend.args[["horiz"]] = TRUE
       
+      # Catch for horizontal ribbon legend spacing
+      if (type=="ribbon" && isTRUE(legend.args[["horiz"]])) {
+        if (legend.args[["pt.lwd"]] == 1) {
+          legend.args[["x.intersp"]] = 1
+        } else {
+          legend.args[["x.intersp"]] = 0.5
+        }
+      }
+      
       lgnd = legend(
         0, 0,
         bty    = legend.args[["n"]],
@@ -599,7 +630,8 @@ plot2.default = function(
           graphics::polygon(
             x = c(split_data[[i]]$x, rev(split_data[[i]]$x)),
             y = c(split_data[[i]]$ymin, rev(split_data[[i]]$ymax)),
-            col = adjustcolor(col[i], 0.2),
+            # col = adjustcolor(col[i], ribbon_alpha),
+            col = bg[i],
             border = FALSE
           )
         }
@@ -829,7 +861,7 @@ plot2.formula = function(
 plot2.density = function(
     x = NULL,
     by = NULL,
-    type = "l",
+    type = c("l", "ribbon"),
     xlim = NULL,
     ylim = NULL,
     # log = "",
@@ -845,10 +877,15 @@ plot2.density = function(
     pch = NULL,
     col = NULL,
     lty = NULL,
+    bg = NULL,
     par_restore = FALSE,
     ...
     ) {
   
+  type = match.arg(type)
+  ## override if bg = "by"
+  if (!is.null(bg)) type = "ribbon"
+  
   object = x
   legend.args = list(x = NULL)
   # Grab by label to pass on legend title to plot2.default
@@ -895,6 +932,16 @@ plot2.density = function(
       xlab = paste0("N = ", n, "   Joint Bandwidth = ", bw)
     }
   }
+  if (type=="ribbon") {
+    ymin = rep(0, length(y))
+    ymax = y
+    # set extra legend params to get bordered boxes with fill
+    legend.args[["x.intersp"]] = 1.25
+    legend.args[["lty"]] = 0
+    legend.args[["pt.lwd"]] = 1
+  } else {
+    ymin = ymax = NULL
+  }
   
   ## axes range
   if (is.null(xlim)) xlim = range(x)
@@ -923,8 +970,11 @@ plot2.density = function(
     legend.args = legend.args,
     pch = pch,
     col = col,
+    bg = bg,
     lty = lty,
     par_restore = par_restore,
+    ymin = ymin,
+    ymax = ymax,
     ...
     )
 
diff --git a/README.Rmd b/README.Rmd
index 37225c2b..5b97026e 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -23,7 +23,7 @@ knitr::opts_chunk$set(
 [![R-CMD-check](https://github.com/grantmcdermott/plot2/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/grantmcdermott/plot2/actions/workflows/R-CMD-check.yaml)
 <!-- badges: end -->
 
-A lightweight extension of the base R `plot` function, with support for
+A lightweight extension of the base R `plot` system, with support for
 automatic grouping and legend handling, and several other enhancements.
 
 ## Installation
@@ -92,26 +92,27 @@ plot2(0:10, main = "plot2")
 ```
 
 Similarly, we can plot elements from a data frame using either the atomic or
-formula methods.
+formula methods. Here's a simple example using a slightly modified version of
+the `airquality` dataset that comes bundled with base R.
 
-```{r base_2}
-par(mfrow = c(2, 2))
-
-plot(airquality$Day, airquality$Temp, main = "plot")
-plot(Temp ~ Day, data = airquality, main = "plot (formula)")
-plot2(airquality$Day, airquality$Temp, main = "plot2")
-plot2(Temp ~ Day, data = airquality, main = "plot2 (formula)")
+```{r plot2_simple}
+aq = airquality
+aq$Month = factor(month.abb[aq$Month], levels = month.abb[5:9])
 
-dev.off() # reset to default (single) plot window
+# with(aq,  plot2(Temp ~ Day)) # atomic method (same as below)
+plot2(Temp ~ Day, data = aq) # formula method
 ```
 
+For the remaining examples below, we'll continue with this modified `aq`
+dataset to demonstrate some of additional features that `plot2` enables.
+
 ### Grouped data
 
-So far, so good. But where `plot2` starts to diverge from its base counterpart is
+Where `plot2` starts to diverge from its base counterpart is
 with respect to grouped data. In particular, `plot2` allows you to characterize
 groups using the `by` argument.^[At this point, experienced base plot users 
 might protest that you _can_ colour by groups using the `col` argument, e.g.
-`with(airquality, plot(Day, Temp, col = Month))`. This is true,
+`with(aq, plot(Day, Temp, col = Month))`. This is true,
 but there are several limitations. First, you don't get an automatic legend.
 Second, the base `plot.formula` method doesn't specify the grouping within the
 formula itself (not a deal-breaker, but not particularly consistent either).
@@ -122,15 +123,15 @@ plots (i.e., type="l"). Instead, you have to transpose your data and use
 old StackOverflow thread for a longer discussion.]
 
 ```{r by}
-# plot2(airquality$Day, airquality$Temp, by = airquality$Month) # same as below
-with(airquality, plot2(Day, Temp, by = Month))
+# plot2($Day, aq$Temp, by = aq$Month) # same as below
+with(aq, plot2(Day, Temp, by = Month))
 ```
 
-An even more convenient approach is to use the equivalent formula syntax. Just
-place the grouping variable after a vertical bar (i.e., `|`).
+An arguably more convenient approach is to use the equivalent formula syntax.
+Just place the grouping variable after a vertical bar (i.e., `|`).
 
 ```{r formula}
-plot2(Temp ~ Day | Month, data = airquality)
+plot2(Temp ~ Day | Month, data = aq)
 ```
 
 You can use standard base plotting arguments to adjust features of your plot.
@@ -138,8 +139,7 @@ For example, change `pch` (plot character) to get filled points.
 
 ```{r pch_16}
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   pch = 16
 )
 ```
@@ -149,43 +149,41 @@ Similarly, converting to a grouped line plot is a simple matter of adjusting the
 
 ```{r type_l}
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l"
 )
 ```
 
-Note that we can automatically adjust both `pch` and `lty` by groups using the
-`"by"` convenience keyword. This can be used in conjunction with the default
+The default behaviour of `plot2` is to represent groups through colour. However,
+note that we can automatically adjust `pch` and `lty` by groups too by passing
+the `"by"` convenience keyword. This can be used in conjunction with the default
 group colouring. Or, as a replacement for group colouring---an option that may
 be particularly useful for contexts where colour is expensive or prohibited
 (e.g., certain academic journals).
 
 ```{r by_lty}
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l",
   col = "black", # override automatic group colours
   lty = "by"     # change line type by group instead
 )
 ```
 
-### Colors
+### Colours
 
 On the subject of group colours, these are easily customized via the `palette`
-argument. The default group colours are inherited from either the "R4" or
-"Viridis" palettes, depending on the number of groups. However, all of the many
-palettes listed by `palette.pals()` and `hcl.pals()` are supported as
-convenience strings.^[See the accompanying help pages of those two functions
+argument. The default group colours are inherited from the user's default
+palette. (Most likely the "R4" set of colors; see `?palette`). However, all of
+the various palettes listed by `palette.pals()` and `hcl.pals()` are supported
+as convenience strings.^[See the accompanying help pages of those two functions
 for more details on the available palettes, or read the
 [article](https://arxiv.org/pdf/2303.04918.pdf) by Achim Zeileis and Paul
 Murrell.] For example:
 
 ```{r palette_tableau}
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l",
   palette = "Tableau 10" # or "ggplot2", "Okabe-Ito", "Set 2", "Harmonic", etc.
 )
@@ -195,64 +193,97 @@ Beyond these convenience strings, users can also supply a valid
 palette-generating function for finer control over transparency, colour order,
 and so forth. We'll see a demonstration of this further below.
 
+To underscore what we said earlier, colours are inherited from the user's
+default palette. So these can also be set globally, just as they can for the
+base `plot` function. The next code chunk will set a new default palette for the
+remainder of the plots that follow.
+
+```{r global_palette}
+# Set the default palette globally via the generic palette function
+palette("Tableau 10")
+```
+
 ### Legend
 
 In all of the preceding plots, you will have noticed that we get an automatic
-legend. The legend position and look can be customized using appropriate
-arguments. You can change (or turn off) the legend title and bounding box,
-switch the direction of the legend text, etc. Below, we particularly draw your
-attention to the trailing "!" in the `legend` position argument. This tells
-`plot2` to place the legend _outside_ the plot area.
+legend. The legend position and look can be customized with the `legend`
+argument. At a minimum, you can pass the familiar legend position keywords as a
+convenience string ("topright", "left", etc.). Moreover, a key feature of
+`plot2` is that we can easily and elegantly place the legend _outside_ the plot
+area by adding a trailing "!" to these keywords. (As you may have realised, the
+default legend position is "right!".) Let's demonstrate by moving the legend to
+the left of the plot:
 
 ```{r legend_bottom}
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l",
-  legend = legend("bottom!", title = "Month of the year", bty = "o")
+  legend = "left!"
 )
 ```
 
-Note that legend position keywords without the exclamation point (i.e., for
-inside the plot area) should still as per normal. Grouped density plot example:
+Beyond the convenience of these positional keywords, the `legend` argument also
+permits additional customization by passing an appropriate function (or, a list
+of arguments that will be passed on to the standard `legend()` function
+internally.) So you can change or turn off the legend title, remove the bounding
+box, switch the direction of the legend text to horizontal, etc. Here's a
+grouped density plot example, where we also add some shading by specifying that
+the background colour should vary by groups too.
 
 ```{r density_topright}
 with(
-  airquality,
+  aq,
   plot2(
     density(Temp),
     by = Month,
-    legend = legend("topright", bty = "o")
+    bg = "by",                           # add background fill by groups
+    grid = TRUE,                         # add background grid
+    legend = list("topright", bty = "o") # change legend features
   )
 )
 ```
 
-### Point ranges, error bars, and ribbon plots
+### Interval plots
+
+`plot2` adds supports for interval plots via the `"pointrange"`, `"errorbar"`,
+`"ribbon"` type arguments. An obvious use-case is for regression analysis and
+prediction.
 
-`plot2` adds supports for uncertainty intervals via the `"pointrange"`,
-`"errorbar"`, `"ribbon"` type arguments. An obvious use-case is for regression
-coefficient plots.
+```{r ribbon_pred}
+mod = lm(Temp ~ 0 + Month / Day, data = aq)
+aq = cbind(aq, predict(mod, interval = "confidence"))
+
+with(
+  aq,
+  plot2(
+    x = Day, y = fit, by = Month,
+    ymin = lwr, ymax = upr,
+    type = "ribbon",
+    grid = TRUE,
+    main = "Model predictions"
+  )
+)
+```
+
+Similarly, we can grab the model estimates to produce nice coefficient plots.
 
 ```{r pointrange, warning = FALSE}
-aq = airquality
-aq$mnth = factor(month.abb[aq$Month], levels = month.abb)
-mod = lm(Temp ~ 0 + mnth, aq)
 coeftab = data.frame(
-  gsub("mnth", "", names(coef(mod))),
+  gsub("Month", "", names(coef(mod))),
   coef(mod),
   confint(mod)
-  ) |>
+) |>
   setNames(c("term", "estimate", "ci_low", "ci_high"))
 
 with(
-  coeftab,
+  subset(coeftab, !grepl("Day", term)),
   plot2(
     x = term, y = estimate,
     ymin = ci_low, ymax = ci_high,
     type = "pointrange", # or: "errobar", "ribbon"
     pch = 19, col = "dodgerblue",
     grid = TRUE,
-    main = "Effect on Temperature"
+    main = "Average Monthly Effect on Temperature"
   )
 )
 ```
@@ -273,12 +304,11 @@ par(
 )
 
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "b",
   palette = palette.colors(palette = "Tableau 10", alpha = 0.5),
-  main = "Daily temperatures by month",
-  frame = FALSE, grid = TRUE
+  frame = FALSE, grid = TRUE,
+  main = "Daily temperatures by month"
 )
 ```
 
@@ -295,8 +325,7 @@ library(basetheme)
 basetheme("royal") # or "clean", "dark", "ink", "brutal", etc.
 
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "b",
   pch = "by",
   palette = "Tropic",
diff --git a/README.md b/README.md
index 4a3dfb37..186e522c 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@ badge](https://grantmcdermott.r-universe.dev/badges/plot2)](https://grantmcdermo
 [![R-CMD-check](https://github.com/grantmcdermott/plot2/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/grantmcdermott/plot2/actions/workflows/R-CMD-check.yaml)
 <!-- badges: end -->
 
-A lightweight extension of the base R `plot` function, with support for
+A lightweight extension of the base R `plot` system, with support for
 automatic grouping and legend handling, and several other enhancements.
 
 ## Installation
@@ -87,45 +87,42 @@ plot2(0:10, main = "plot2")
 <img src="man/figures/README-base_1-1.png" width="100%" />
 
 Similarly, we can plot elements from a data frame using either the
-atomic or formula methods.
+atomic or formula methods. Here’s a simple example using a slightly
+modified version of the `airquality` dataset that comes bundled with
+base R.
 
 ``` r
-par(mfrow = c(2, 2))
+aq = airquality
+aq$Month = factor(month.abb[aq$Month], levels = month.abb[5:9])
 
-plot(airquality$Day, airquality$Temp, main = "plot")
-plot(Temp ~ Day, data = airquality, main = "plot (formula)")
-plot2(airquality$Day, airquality$Temp, main = "plot2")
-plot2(Temp ~ Day, data = airquality, main = "plot2 (formula)")
+# with(aq,  plot2(Temp ~ Day)) # atomic method (same as below)
+plot2(Temp ~ Day, data = aq) # formula method
 ```
 
-<img src="man/figures/README-base_2-1.png" width="100%" />
-
-``` r
+<img src="man/figures/README-plot2_simple-1.png" width="100%" />
 
-dev.off() # reset to default (single) plot window
-#> null device 
-#>           1
-```
+For the remaining examples below, we’ll continue with this modified `aq`
+dataset to demonstrate some of additional features that `plot2` enables.
 
 ### Grouped data
 
-So far, so good. But where `plot2` starts to diverge from its base
-counterpart is with respect to grouped data. In particular, `plot2`
-allows you to characterize groups using the `by` argument.[^1]
+Where `plot2` starts to diverge from its base counterpart is with
+respect to grouped data. In particular, `plot2` allows you to
+characterize groups using the `by` argument.[^1]
 
 ``` r
-# plot2(airquality$Day, airquality$Temp, by = airquality$Month) # same as below
-with(airquality, plot2(Day, Temp, by = Month))
+# plot2($Day, aq$Temp, by = aq$Month) # same as below
+with(aq, plot2(Day, Temp, by = Month))
 ```
 
 <img src="man/figures/README-by-1.png" width="100%" />
 
-An even more convenient approach is to use the equivalent formula
+An arguably more convenient approach is to use the equivalent formula
 syntax. Just place the grouping variable after a vertical bar (i.e.,
 `|`).
 
 ``` r
-plot2(Temp ~ Day | Month, data = airquality)
+plot2(Temp ~ Day | Month, data = aq)
 ```
 
 <img src="man/figures/README-formula-1.png" width="100%" />
@@ -135,8 +132,7 @@ plot. For example, change `pch` (plot character) to get filled points.
 
 ``` r
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   pch = 16
 )
 ```
@@ -148,24 +144,24 @@ adjusting the `type` argument.
 
 ``` r
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l"
 )
 ```
 
 <img src="man/figures/README-type_l-1.png" width="100%" />
 
-Note that we can automatically adjust both `pch` and `lty` by groups
-using the `"by"` convenience keyword. This can be used in conjunction
-with the default group colouring. Or, as a replacement for group
-colouring—an option that may be particularly useful for contexts where
-colour is expensive or prohibited (e.g., certain academic journals).
+The default behaviour of `plot2` is to represent groups through colour.
+However, note that we can automatically adjust `pch` and `lty` by groups
+too by passing the `"by"` convenience keyword. This can be used in
+conjunction with the default group colouring. Or, as a replacement for
+group colouring—an option that may be particularly useful for contexts
+where colour is expensive or prohibited (e.g., certain academic
+journals).
 
 ``` r
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l",
   col = "black", # override automatic group colours
   lty = "by"     # change line type by group instead
@@ -174,18 +170,18 @@ plot2(
 
 <img src="man/figures/README-by_lty-1.png" width="100%" />
 
-### Colors
+### Colours
 
 On the subject of group colours, these are easily customized via the
-`palette` argument. The default group colours are inherited from either
-the “R4” or “Viridis” palettes, depending on the number of groups.
-However, all of the many palettes listed by `palette.pals()` and
-`hcl.pals()` are supported as convenience strings.[^2] For example:
+`palette` argument. The default group colours are inherited from the
+user’s default palette. (Most likely the “R4” set of colors; see
+`?palette`). However, all of the various palettes listed by
+`palette.pals()` and `hcl.pals()` are supported as convenience
+strings.[^2] For example:
 
 ``` r
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l",
   palette = "Tableau 10" # or "ggplot2", "Okabe-Ito", "Set 2", "Harmonic", etc.
 )
@@ -197,70 +193,106 @@ Beyond these convenience strings, users can also supply a valid
 palette-generating function for finer control over transparency, colour
 order, and so forth. We’ll see a demonstration of this further below.
 
+To underscore what we said earlier, colours are inherited from the
+user’s default palette. So these can also be set globally, just as they
+can for the base `plot` function. The next code chunk will set a new
+default palette for the remainder of the plots that follow.
+
+``` r
+# Set the default palette globally via the generic palette function
+palette("Tableau 10")
+```
+
 ### Legend
 
 In all of the preceding plots, you will have noticed that we get an
-automatic legend. The legend position and look can be customized using
-appropriate arguments. You can change (or turn off) the legend title and
-bounding box, switch the direction of the legend text, etc. Below, we
-particularly draw your attention to the trailing “!” in the `legend`
-position argument. This tells `plot2` to place the legend *outside* the
-plot area.
+automatic legend. The legend position and look can be customized with
+the `legend` argument. At a minimum, you can pass the familiar legend
+position keywords as a convenience string (“topright”, “left”, etc.).
+Moreover, a key feature of `plot2` is that we can easily and elegantly
+place the legend *outside* the plot area by adding a trailing “!” to
+these keywords. (As you may have realised, the default legend position
+is “right!”.) Let’s demonstrate by moving the legend to the left of the
+plot:
 
 ``` r
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "l",
-  legend = legend("bottom!", title = "Month of the year", bty = "o")
+  legend = "left!"
 )
 ```
 
 <img src="man/figures/README-legend_bottom-1.png" width="100%" />
 
-Note that legend position keywords without the exclamation point (i.e.,
-for inside the plot area) should still as per normal. Grouped density
-plot example:
+Beyond the convenience of these positional keywords, the `legend`
+argument also permits additional customization by passing an appropriate
+function (or, a list of arguments that will be passed on to the standard
+`legend()` function internally.) So you can change or turn off the
+legend title, remove the bounding box, switch the direction of the
+legend text to horizontal, etc. Here’s a grouped density plot example,
+where we also add some shading by specifying that the background colour
+should vary by groups too.
 
 ``` r
 with(
-  airquality,
+  aq,
   plot2(
     density(Temp),
     by = Month,
-    legend = legend("topright", bty = "o")
+    bg = "by",                           # add background fill by groups
+    grid = TRUE,                         # add background grid
+    legend = list("topright", bty = "o") # change legend features
   )
 )
 ```
 
 <img src="man/figures/README-density_topright-1.png" width="100%" />
 
-### Point ranges, error bars, and ribbon plots
+### Interval plots
 
-`plot2` adds supports for uncertainty intervals via the `"pointrange"`,
+`plot2` adds supports for interval plots via the `"pointrange"`,
 `"errorbar"`, `"ribbon"` type arguments. An obvious use-case is for
-regression coefficient plots.
+regression analysis and prediction.
+
+``` r
+mod = lm(Temp ~ 0 + Month / Day, data = aq)
+aq = cbind(aq, predict(mod, interval = "confidence"))
+
+with(
+  aq,
+  plot2(
+    x = Day, y = fit, by = Month,
+    ymin = lwr, ymax = upr,
+    type = "ribbon",
+    grid = TRUE,
+    main = "Model predictions"
+  )
+)
+```
+
+<img src="man/figures/README-ribbon_pred-1.png" width="100%" />
+
+Similarly, we can grab the model estimates to produce nice coefficient
+plots.
 
 ``` r
-aq = airquality
-aq$mnth = factor(month.abb[aq$Month], levels = month.abb)
-mod = lm(Temp ~ 0 + mnth, aq)
 coeftab = data.frame(
-  gsub("mnth", "", names(coef(mod))),
+  gsub("Month", "", names(coef(mod))),
   coef(mod),
   confint(mod)
-  ) |>
+) |>
   setNames(c("term", "estimate", "ci_low", "ci_high"))
 
 with(
-  coeftab,
+  subset(coeftab, !grepl("Day", term)),
   plot2(
     x = term, y = estimate,
     ymin = ci_low, ymax = ci_high,
     type = "pointrange", # or: "errobar", "ribbon"
     pch = 19, col = "dodgerblue",
     grid = TRUE,
-    main = "Effect on Temperature"
+    main = "Average Monthly Effect on Temperature"
   )
 )
 ```
@@ -284,12 +316,11 @@ par(
 )
 
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "b",
   palette = palette.colors(palette = "Tableau 10", alpha = 0.5),
-  main = "Daily temperatures by month",
-  frame = FALSE, grid = TRUE
+  frame = FALSE, grid = TRUE,
+  main = "Daily temperatures by month"
 )
 ```
 
@@ -309,8 +340,7 @@ library(basetheme)
 basetheme("royal") # or "clean", "dark", "ink", "brutal", etc.
 
 plot2(
-  Temp ~ Day | Month,
-  data = airquality,
+  Temp ~ Day | Month, data = aq,
   type = "b",
   pch = "by",
   palette = "Tropic",
@@ -340,14 +370,13 @@ is expensive (e.g., an R application running in a browser via
 
 [^1]: At this point, experienced base plot users might protest that you
     *can* colour by groups using the `col` argument, e.g.
-    `with(airquality, plot(Day, Temp, col = Month))`. This is true, but
-    there are several limitations. First, you don’t get an automatic
-    legend. Second, the base `plot.formula` method doesn’t specify the
-    grouping within the formula itself (not a deal-breaker, but not
-    particularly consistent either). Third, and perhaps most
-    importantly, this grouping doesn’t carry over to line plots (i.e.,
-    type=“l”). Instead, you have to transpose your data and use
-    `matplot`. See
+    `with(aq, plot(Day, Temp, col = Month))`. This is true, but there
+    are several limitations. First, you don’t get an automatic legend.
+    Second, the base `plot.formula` method doesn’t specify the grouping
+    within the formula itself (not a deal-breaker, but not particularly
+    consistent either). Third, and perhaps most importantly, this
+    grouping doesn’t carry over to line plots (i.e., type=“l”). Instead,
+    you have to transpose your data and use `matplot`. See
     [this](https://stackoverflow.com/questions/10519873/how-to-create-a-line-plot-with-groups-in-base-r-without-loops)
     old StackOverflow thread for a longer discussion.
 
diff --git a/man/figures/README-base_1-1.png b/man/figures/README-base_1-1.png
index 0fa1d112..5ecd2086 100644
Binary files a/man/figures/README-base_1-1.png and b/man/figures/README-base_1-1.png differ
diff --git a/man/figures/README-base_2-1.png b/man/figures/README-base_2-1.png
deleted file mode 100644
index f8cf1046..00000000
Binary files a/man/figures/README-base_2-1.png and /dev/null differ
diff --git a/man/figures/README-basethme_royal-1.png b/man/figures/README-basethme_royal-1.png
index b2b25bcf..f5bf1c86 100644
Binary files a/man/figures/README-basethme_royal-1.png and b/man/figures/README-basethme_royal-1.png differ
diff --git a/man/figures/README-by-1.png b/man/figures/README-by-1.png
index e09f36c5..2ac5478d 100644
Binary files a/man/figures/README-by-1.png and b/man/figures/README-by-1.png differ
diff --git a/man/figures/README-by_lty-1.png b/man/figures/README-by_lty-1.png
index 50aa8e5e..5bc49109 100644
Binary files a/man/figures/README-by_lty-1.png and b/man/figures/README-by_lty-1.png differ
diff --git a/man/figures/README-density_topright-1.png b/man/figures/README-density_topright-1.png
index fa4ab0da..73aa20db 100644
Binary files a/man/figures/README-density_topright-1.png and b/man/figures/README-density_topright-1.png differ
diff --git a/man/figures/README-desnity_topright-1.png b/man/figures/README-desnity_topright-1.png
deleted file mode 100644
index 7e33b727..00000000
Binary files a/man/figures/README-desnity_topright-1.png and /dev/null differ
diff --git a/man/figures/README-formula-1.png b/man/figures/README-formula-1.png
index e09f36c5..2ac5478d 100644
Binary files a/man/figures/README-formula-1.png and b/man/figures/README-formula-1.png differ
diff --git a/man/figures/README-hershey_plus-1.png b/man/figures/README-hershey_plus-1.png
index 953a1ab5..616625cd 100644
Binary files a/man/figures/README-hershey_plus-1.png and b/man/figures/README-hershey_plus-1.png differ
diff --git a/man/figures/README-legend_bottom-1.png b/man/figures/README-legend_bottom-1.png
index 7b2807c1..66192d70 100644
Binary files a/man/figures/README-legend_bottom-1.png and b/man/figures/README-legend_bottom-1.png differ
diff --git a/man/figures/README-palette_tableau-1.png b/man/figures/README-palette_tableau-1.png
index 5673c007..6ec52746 100644
Binary files a/man/figures/README-palette_tableau-1.png and b/man/figures/README-palette_tableau-1.png differ
diff --git a/man/figures/README-pch_16-1.png b/man/figures/README-pch_16-1.png
index c6577a60..967500a0 100644
Binary files a/man/figures/README-pch_16-1.png and b/man/figures/README-pch_16-1.png differ
diff --git a/man/figures/README-plot2_simple-1.png b/man/figures/README-plot2_simple-1.png
new file mode 100644
index 00000000..7bdd3e1c
Binary files /dev/null and b/man/figures/README-plot2_simple-1.png differ
diff --git a/man/figures/README-pointrange-1.png b/man/figures/README-pointrange-1.png
index eaad63be..add9dda2 100644
Binary files a/man/figures/README-pointrange-1.png and b/man/figures/README-pointrange-1.png differ
diff --git a/man/figures/README-ribbon_pred-1.png b/man/figures/README-ribbon_pred-1.png
new file mode 100644
index 00000000..626398f5
Binary files /dev/null and b/man/figures/README-ribbon_pred-1.png differ
diff --git a/man/figures/README-type_l-1.png b/man/figures/README-type_l-1.png
index 7e97e3af..c790748c 100644
Binary files a/man/figures/README-type_l-1.png and b/man/figures/README-type_l-1.png differ
diff --git a/man/plot2-package.Rd b/man/plot2-package.Rd
index f72eee0c..df4c0936 100644
--- a/man/plot2-package.Rd
+++ b/man/plot2-package.Rd
@@ -38,5 +38,10 @@ Authors:
   \item Achim Zeileis \email{Achim.Zeileis@R-project.org} (\href{https://orcid.org/0000-0003-0918-3766}{ORCID})
 }
 
+Other contributors:
+\itemize{
+  \item Etienne Bacher \email{etienne.bacher@protonmail.com} [contributor]
+}
+
 }
 \keyword{internal}
diff --git a/man/plot2.Rd b/man/plot2.Rd
index 75ee03fd..de1c1bf0 100644
--- a/man/plot2.Rd
+++ b/man/plot2.Rd
@@ -37,6 +37,7 @@ plot2(x, ...)
   par_restore = FALSE,
   ymin = NULL,
   ymax = NULL,
+  ribbon_alpha = 0.2,
   ...
 )
 
@@ -69,7 +70,7 @@ plot2(x, ...)
 \method{plot2}{density}(
   x = NULL,
   by = NULL,
-  type = "l",
+  type = c("l", "ribbon"),
   xlim = NULL,
   ylim = NULL,
   main = NULL,
@@ -84,6 +85,7 @@ plot2(x, ...)
   pch = NULL,
   col = NULL,
   lty = NULL,
+  bg = NULL,
   par_restore = FALSE,
   ...
 )
@@ -111,7 +113,9 @@ lines, "o" for overplotted points and lines, "s" and "S" for stair steps
 and "h" for histogram-like vertical lines. "n" does not produce
 any points or lines.
 \item Additional plot2 types: "pointrange", "errorbar", and "ribbon" for
-drawing these respective plot types.
+drawing these interval plot types. (Note that specifying "ribbon" for
+objects of class \code{density} will yield a density plot with a shaded
+interior.)
 }}
 
 \item{xlim}{the x limits (x1, x2) of the plot. Note that x1 > x2 is allowed
@@ -234,7 +238,13 @@ been altered to make space for the \code{plot2} legend. Users can opt out of
 this persistent behaviour by setting to TRUE instead. (Another option would
 be calling \code{dev.off()} to reset all \code{par} settings to their defaults.)}
 
-\item{ymin, ymax}{minimum and maximum coordinates of the point range. Only used when \code{type="pointrange"}.}
+\item{ymin, ymax}{minimum and maximum coordinates of interval plot types. Only
+used when the \code{type} argument is one of "pointrange", "errorbar", or
+"ribbon".}
+
+\item{ribbon_alpha}{numeric factor modifying the opacity alpha of any ribbon
+shading; typically in \verb{[0, 1]}. Default value is 0.2. Only used when
+\code{type = "ribbon"}.}
 
 \item{formula}{a \code{formula} that may also include a grouping variable after a
 "|", such as \code{y ~ x | z}. Note that the \code{formula} and \code{x} arguments should
