ggdensity: Highest Density Region Visualizations for ggplot2

You can install ggdensity from CRAN or the latest development version from GitHub.

From CRAN:

install.packages("ggdensity")

From GitHub (Development version):

if (!requireNamespace("remotes")) install.packages("remotes")
remotes::install_github("jamesotto852/ggdensity")

Because geom_hdr() and related functions rely on a discretization of the density estimate, the estimated HDRs are computed on a grid over the x-y aesthetic space. If the plot limits are too tight (based only on observed data points), a significant portion of the density might be excluded, biasing the HDRs.

Recommendation: If you see an HDR being truncated by the plot window:

1. Manually increase the aesthetic limits using xlim and ylim within the geom or via scale_*_continuous().
2. Use coord_cartesian() to zoom into the desired area without re-calculating the density or changing the underlying limits.

# Recommended pattern to prevent truncation
ggplot(cars, aes(speed, dist)) +
  geom_hdr(xlim = c(-20, 50), ylim = c(-40, 140)) +
  geom_point(color = "red") +
  scale_y_continuous(breaks = 25*(0:5)) +
  coord_cartesian(xlim = c(4, 25), ylim = c(-1, 120))

You can use geom_hdr_fun() to plot parametric densities estimated outside of the ggdensity framework (e.g., via Maximum Likelihood Estimation).

To do this:

1. Fit your distribution using your preferred method to obtain parameter estimates.
2. Define a function f(x, y, ...) that represents the parametric density.
3. Pass this function to geom_hdr_fun() using the fun argument and provide the estimated parameters via the args argument (as a list).

# 1. Define the parametric density function
f <- function(x, y, th) dexp(x, th[1]) * dexp(y, th[2])

# 2. Assume th_hat contains your estimated parameters
# 3. Plot using geom_hdr_fun
ggplot(df, aes(x, y)) +
  geom_hdr_fun(fun = f, args = list(th = th_hat)) +
  geom_point(shape = 21, fill = "lightgreen", alpha = .25) +
  coord_equal()

Because geom_hdr() maps probability to the alpha aesthetic, you can use fill and color aesthetics to map to variables and visualize different subpopulations (e.g., different species in a dataset).

library("palmerpenguins")
library("ggplot2")
library("ggdensity")

# Visualizing subpopulations using fill
ggplot(penguins, aes(flipper_length_mm, bill_length_mm, fill = species)) +
  geom_hdr(xlim = c(160, 240), ylim = c(30, 70)) +
  geom_point(shape = 21)

# Visualizing subpopulations using lines to reduce overplotting
ggplot(penguins, aes(flipper_length_mm, bill_length_mm, color = species)) +
  geom_hdr_lines(xlim = c(160, 240), ylim = c(30, 70)) +
  geom_point(size = 1)

If you already have a known joint Probability Density Function (PDF) as an R function, use geom_hdr_fun() or geom_hdr_lines_fun(). These functions behave similarly to ggplot2::geom_function(), accepting a fun argument that specifies the PDF to be summarized.

f <- function(x, y) dnorm(x) * dgamma(y, 5, 3)

ggplot() +
  geom_hdr_fun(fun = f, xlim = c(-4, 4), ylim = c(0, 5))

The underlying stat in geom_hdr() and geom_hdr_lines() computes a variable named probs. You can map this variable to other aesthetics (like fill or color) using ggplot2::after_stat().

Note: If you override the default alpha mapping, you must explicitly set alpha = 1 to avoid unexpected transparency behavior.

ggplot(faithful, aes(eruptions, waiting)) +
  geom_hdr(
    aes(fill = after_stat(probs)), 
    alpha = 1, 
    xlim = c(0, 8), 
    ylim = c(30, 110)
  )

The geom_hdr_points() function provides a scatterplot geom that allows individual data points to be seen simultaneously with High Density Regions (HDRs). This is particularly useful for addressing significant overplotting in datasets.

p_hdr_points <- ggplot(diamonds, aes(carat, price)) +
  geom_hdr_points()

You can access the underlying numerical summaries of estimated densities and HDRs using these functions:

• get_hdr(df, method = "kde"): Used for bivariate data. Returns a list containing:
  • $df_est: A data frame with columns x, y, fhat (density estimate), fhat_discretized, and hdr (membership indicator).
  • $breaks: Named numeric vector of probability breaks (e.g., "99%", "95%").
  • $data: A data frame containing the original data points with an hdr_membership column.
• get_hdr_1d(x, method = "kde"): Used for univariate data. Returns a similar list structure optimized for 1D data.

df <- data.frame(x = rnorm(1e3), y = rnorm(1e3))
res <- get_hdr(df, method = "kde")
# res contains $df_est, $breaks, and $data

You can control how the underlying bivariate distribution is estimated using the method argument in geom_hdr() and geom_hdr_lines().

Methods can be specified as:

1. A character vector (e.g., method = "kde").
2. A function call (e.g., method = method_kde()).

When using function calls like method_kde(), you can pass parameters such as h or adjust to control the kernel's bandwidth. For more details, see ?method_kde or vignette("method", "ggdensity").

Because geom_hdr() relies on a discretization of the density (a grid), the HDRs might be truncated if the plot window is too tight. If you see HDRs being clipped at the edges of your plot, follow these steps:

1. Manually increase the xlim and ylim within the geom_hdr() call.
2. Use coord_cartesian(xlim = ..., ylim = ...) to zoom into the region of interest without changing the underlying data scale or clipping the density computation.

ggplot(cars, aes(speed, dist)) +
  geom_hdr(xlim = c(-20, 50), ylim = c(-40, 140)) +
  geom_point() +
  coord_cartesian(xlim = c(4, 25), ylim = c(-1, 120))