Smarter axes limits around stat_density_2d polygons

[vlepori]

Axes limits are based on points rather than on the computed contour lines/shapes. This leads to some weird cropping of polygons:

faithful %>%
  mutate(gr = if_else(row_number()%%2==0,1,2)) %>% 
  ggplot(aes(x = eruptions, y = waiting)) +
  stat_density_2d(aes(alpha = after_stat(level)), geom = "polygon")+ 
  geom_density2d(color="black")+
  facet_wrap(~gr)+
  theme_classic()

A fiddly way around this is to change the limits manually with blank layers:

faithful %>%
  mutate(gr = if_else(row_number()%%2==0,1,2)) %>% 
  mutate(miny=waiting-10,maxy=waiting+10,minx = eruptions-1,maxx=eruptions+1) %>% 
  ggplot(aes(x = eruptions, y = waiting)) +
  geom_density2d(color="black")+
  stat_density_2d(aes(alpha = after_stat(level)), geom = "polygon")+ 
  facet_wrap(~gr)+
  theme_classic()+
  geom_blank(aes(y = miny)) +
  geom_blank(aes(y = maxy)) +
  geom_blank(aes(x = minx)) +
  geom_blank(aes(x = maxx))

but this becomes quickly a nightmare if one has to do it often (even without hardcoding values). It would be nice to have the axes automatically adjust to the size of the largest polygon. Alternatively, polygon cropping should be handled differently.
Thank you, and apologies if this is already well known and/or has a solution already!

[clauswilke]

The next ggplot2 release will have a geom for filled densities. I think this addresses your concerns.

library(tidyverse)

faithful %>%
  mutate(gr = if_else(row_number()%%2==0,1,2)) %>% 
  ggplot(aes(x = eruptions, y = waiting)) +
  geom_density2d_filled(aes(alpha = after_stat(level)), fill = "black") + 
  geom_density2d(color="black") +
  facet_wrap(~gr) +
  theme_classic()

^{Created on 2020-04-24 by the reprex package (v0.3.0)}

[clauswilke]

If you don't like the white bands you can turn off scale expansion.

library(tidyverse)

faithful %>%
  mutate(gr = if_else(row_number()%%2==0,1,2)) %>% 
  ggplot(aes(x = eruptions, y = waiting)) +
  geom_density2d_filled(aes(alpha = after_stat(level)), fill = "black") + 
  geom_density2d(color="black") +
  facet_wrap(~gr) +
  coord_cartesian(expand = FALSE) +
  theme_classic()

^{Created on 2020-04-24 by the reprex package (v0.3.0)}

[vlepori]

You guys are always one step ahead! This looks awesome and definitely fixes the issue, looking forward to it. Only comment: shouldn't the first level be white by default?

[thomasp85]

this can be controlled with the scale_alpha(). geoms have no control over the scale settings

[stephanbitterwolf]

Hello, I wanted to create a plot similar to the original poster where the x and y axes are expanded a bit so that I can see where the heatmap might go to 0. However, this was not possible without adding the geom_blank elements as the original poster did. I don't mind doing it manually. However, I was curious why it isn't the default for stat_density_2d?

[teunbrand]

If people want to extrapolate, they can do so, but extrapolating by default sounds like dangerous territory.

[stephanbitterwolf]

If people want to extrapolate, they can do so, but extrapolating by default sounds like dangerous territory.

Hi Teun, the default currently plots the extent of the data and does not extrapolate. It looks bad because the points are right on the border (you can run the code below to see what I mean). I wanted to zoom out a little bit but found the plot looked worse when I did that because the area remains white. It should be black. In my case, the area truly lacks points, and therefore, values should be 0 and colored black. I am OK with adding geom_blank(), but I was curious if there was a better fix for that.

`library(ggplot2)
library(viridis)
library(dplyr)

Simulated data for demonstration

set.seed(42)
pcoa_data <- data.frame(
Axis.1 = rnorm(100, 0, 0.2),
Axis.2 = rnorm(100, 0, 0.2),
genotype = sample(c("A", "B", "C", "D"), 100, replace = TRUE),
treatment = sample(c("Control", "Solvent", "Atrazine", "Estrone", "A+E"), 100, replace = TRUE)
)

Axis ranges for zooming out

x_range <- range(pcoa_data$Axis.1) + c(-0.1, 0.1)
y_range <- range(pcoa_data$Axis.2) + c(-0.1, 0.1)

Heatmap 1: Area not expanded during density calculation

Heatmap_No_Expansion <- ggplot(pcoa_data, aes(x = Axis.1, y = Axis.2)) +
stat_density_2d(
aes(fill = ..density..),
geom = "raster",
contour = FALSE,
n = 500,
h = c(0.2, 0.2)
) +
scale_fill_viridis_c(option = "magma") +
geom_point(aes(color = genotype, shape = treatment), size = 2, alpha = 0.8) +
labs(
x = "PCoA Axis 1",
y = "PCoA Axis 2",
fill = "Density",
color = "Genotype",
shape = "Treatment"
) +
coord_cartesian(expand = FALSE) +
theme_pubr() +
theme(
legend.position = "right",
plot.title = element_text(hjust = 0.5)
)

Heatmap 2: Area expanded showing enlarged white border

Heatmap_With_Expansion <- ggplot(pcoa_data, aes(x = Axis.1, y = Axis.2)) +
stat_density_2d(
aes(fill = ..density..),
geom = "raster",
contour = FALSE,
n = 500,
h = c(0.2, 0.2)
) +
scale_fill_viridis_c(option = "magma") +
geom_point(aes(color = genotype, shape = treatment), size = 2, alpha = 0.8) +
labs(
x = "PCoA Axis 1",
y = "PCoA Axis 2",
fill = "Density",
color = "Genotype",
shape = "Treatment"
) +
coord_cartesian(xlim = x_range, ylim = y_range, expand = FALSE) +
theme_pubr() +
theme(
legend.position = "right",
plot.title = element_text(hjust = 0.5)
)

Display both heatmaps

print(Heatmap_No_Expansion)
print(Heatmap_With_Expansion)`

[teunbrand]

You can set the scale limits instead of the coord limits. Scale limits can affect data, while coord limits shouldn't.