Kernel density estimates for tidy and geospatial data in the eks package
Introduction
Kernel smoothers are essential tools for data analysis due to their
ability to convey complex statistical information with concise graphical
visualisations. The most widely used kernel smoother is the kernel
density estimator (KDE), though there remain some important gaps in the
implementation in R for specialised data types, most
notably for tibbles (tidy data) within the tidyverse, and for simple
features (geospatial data) within Geographical Information Systems (GIS)
analysis. The tidy_kde and st_kde functions in
the eks package fills in these gaps.
Tidy data example
The data set we focus on is the crabs data set from the
MASS package, with the variables FL frontal
lobe size (mm), CW carapace width (mm) and sp
species (B for blue, O for orange).
library(eks)
library(colorspace)
library(ggplot2)
library(dplyr)
## crabs data
data(crabs, package="MASS")
crabs2 <- select(crabs, FL, CW)
xlab <- "Frontal lobe size (mm)"
ylab <- "Carapace width (mm)"The KDE for tidy data is computed by tidy_kde. From the
output, the scatter plot of the data is generated by
geom_point_ks and the contour plot of the KDE by
geom_contour_ks. The bimodal structure of the data
distribution, corresponding to the two species, is clearly visible from
the KDE plot from tidy_kde. This is due to the optimal
choice of the matrix of smoothing parameters. This optimal smoothing
matrix is the plug-in bandwidth computed by ks::Hpi, and it
is suitable for a wide range of data sets. For further details of the
computation of the kernel density estimate and the bandwidth, see
?ks::kde and ?ks::Hpi.
## KDE contour plot + scatter plot
tkde2 <- tidy_kde(crabs2)
gkde2 <- ggplot(tkde2, aes(x=FL, y=CW)) + labs(x=xlab, y=ylab)
gkde2 + geom_point_ks(colour=8) + geom_contour_ks(colour=1)
On the other hand, the default bandwidth and the resulting KDE
computed by ggplot2::geom_density_2d leads to an
oversmoothed KDE which does not reveal the data bimodality.
## geom_density_2d KDE contour plot + scatter plot
mkde2 <- ggplot(crabs2, aes(x=FL, y=CW))
mkde2 + geom_point(colour=8) + geom_density_2d(colour=1, bins=4)
The default choice of the contour levels in the eks
package is based on probability contours. Probability contours offer an
intuitive approach to selecting the contour levels that reveal the
pertinent characteristics of the data distribution. See Chacon &
Duong (2018, Chapter 2.2). Filled contour plots, generated by
geom_contour_filled_ks, can be coloured with an appropriate
sequential colour scale. For example, a 25% contour region (dark purple
region) is the smallest region that contains 25% of the probability mass
defined by the KDE. The 50% contour region consists of the union of the
light purple region and the dark purple region, and it contains 50% of
the data points etc. Note that the 25% and 50% contour regions of the
crabs KDE are composed of separate, unconnected contour
sub-regions.
## KDE filled contour plot
gkde2 + geom_contour_filled_ks(colour=1) + scale_fill_discrete_sequential(h1=275)
As an alternative to these discretised contours, the usual
ggplot2::geom_raster generates a plot with a continuous
colour scale.
## KDE continuous colour plot
ggplot(tkde2, aes(x=FL,y=CW)) + geom_raster(aes(fill=estimate), interpolate=TRUE) +
labs(x=xlab, y=ylab) + scale_fill_continuous_sequential(h1=275)
One of the main advantages of ggplot2 is its ability to
handle multiple related plots, in this case, KDE plots for each species.
The KDE with blue contours is for the B species, and orange
contours for the O species.
crabs2g <- select(crabs, FL, CW, sp)
crabs2g <- group_by(crabs2g, sp)
tkde2g <- tidy_kde(crabs2g)
gkde2g <- ggplot(tkde2g, aes(x=FL, y=CW, group=sp)) + labs(x=xlab, y=ylab,
colour="Species") + scale_colour_manual(values=c(4, 7))
## superposed KDE contour plots + scatter plots
gkde2g + geom_point_ks(colour=8) +
geom_contour_ks(aes(colour=sp)) + guides(colour=guide_legend(title="Species"))
We can also display each species KDE on its own set of axes.
## facetted KDE contour plots + scatter plots
gkde2g + geom_point_ks(colour=8) +
guides(colour=guide_legend(reverse=FALSE)) +
geom_contour_ks(aes(colour=sp)) + facet_wrap(~sp) 
The probability contour levels computed in
geom_contour_ks and geom_contour_filled_ks are
relative to the grouping variable. So whilst the same probability 25%
level is applied to both groups KDE, the height of 25% contour region
for the blue species is 0.0414, and for the orange species it is 0.0318.
For a direct comparison, it is useful to have a set of fixed contour
heights for all KDEs . A heuristic solution is implemented in
contour_breaks. For the crabs data, this gives
0.0142, 0.0264, 0.0414. Since the KDE for the B species
exceeds the highest level 0.0414, whereas the O KDE doesn’t
reach this is level, the former KDE is more peaked.
## facetted KDE filled contour plots with fixed contour levels for all facets
bkde2g <- contour_breaks(tkde2g)
gkde2g + geom_contour_filled_ks(breaks=bkde2g, colour=1) +
scale_fill_discrete_sequential(h1=275) + facet_wrap(~sp)
Geospatial data example
GIS for geospatial data analysis in R is implemented in
the sf package, and the eks package builds on
this. To illustrate geospatial KDE, we focus on the
grevilleasf data set in the eks package. It
has 22303 rows, where each row corresponds to an observed grevillea
plant in Western Australia. In addition, we utilise wa, the
geospatial polygon for Western Australia. Both of these geospatial data
sets are in the EPSG:7850 (GDA2020/MGA zone 50) projection.
library(sf)
## Grevillea data
data(grevilleasf, package="eks")
grevilleasf <- mutate(grevilleasf, species=factor(species))
paradoxa <- filter(grevilleasf, name %in% "Grevillea paradoxa")
eryngioides <- filter(grevilleasf, name %in% "Grevillea eryngioides")
grevillea_ep <- filter(grevilleasf, name %in% c("Grevillea eryngioides",
"Grevillea paradoxa"))
grevillea_ep <- group_by(grevillea_ep, name)
xlim <- c(165479.3, 1096516.3); ylim <- c(6101931, 7255991)
## WA polygon
data(wa, package="eks")
gwa <- geom_sf(data=wa, fill=NA, colour=1)Since geospatial data can be visualised with both base R
and ggplot2 graphics engines, we provide code for both:
their outputs are similar due to the standardisation of geospatial maps
within GIS. Though these plots can’t be mixed due to fundamental
differences between the graphical rendering in base R and
ggplot2. This is the reason that, for example, adding a
scale bar requires functions from different packages.
## base R scatter plot
plot(st_geometry(wa), xlim=xlim, ylim=ylim)
plot(st_geometry(eryngioides), add=TRUE, col=3, pch=16, cex=0.5)
plot(st_geometry(paradoxa), add=TRUE, col=6, pch=17, cex=0.5)
mapsf::mf_legend(type="symb", val=c("Grevillea eryngioides", "Grevillea paradoxa"),
pal=c(3,6), pch=16:17, cex=c(1,1), title="Species", pos="topright")
mapsf::mf_scale(size=200, lwd=4)
## geom_sf scatter plot
gsc <- ggspatial::annotation_scale(data=data.frame(name="Grevillea paradoxa"),
location="br", width_hint=0.2, bar_cols=1)
theme_set(ggthemes::theme_map())
theme_update(legend.position=c(0.99,0.99), legend.justification=c(1,1))
ggplot() + gwa + gsc +
geom_sf(data=grevillea_ep, aes(colour=name, shape=name), size=0.5) +
coord_sf(xlim=xlim, ylim=ylim) + scale_colour_manual(values=c(3, 6)) +
guides(colour=guide_legend(title="Species"), shape=guide_legend(title="Species"))
The KDE for geospatial data is computed by st_kde. The
calculations of the KDE, including the bandwidth matrix fo smoothing
parameters, are the same as in tidy_kde. Though, unlike for
tidy_kde where the probability contour regions are computed
dynamically in geom_contour_filled_ks, the 1% to 99%
regions are explicitly computed as multipolygons in st_kde
since this conversion can be computationally heavy to execute for each
plot. For display, it is a matter of selecting the appropriate contour
regions. The quartile contours 25%, 50%, 75% are selected by default in
geom_contour_filled_ks for tidy data. This is also the case
for the base R plot.
skde1 <- st_kde(paradoxa)
## base R contour plot
plot(st_geometry(wa), xlim=xlim, ylim=ylim)
plot(st_geometry(paradoxa), add=TRUE, pch=16, col=8, cex=0.5)
plot(skde1, add=TRUE, col=NA, border=1, legend=FALSE)
mapsf::mf_scale(size=200, lwd=4)
On the other hand, we can’t replicate exactly the default contour
selection in geom_sf, so we first apply the auxiliary
function st_get_contour to the input of
geom_sf.
## geom_sf contour plot
gs <- ggplot(skde1) + gwa + gsc
gs + geom_sf(data=paradoxa, col=8, size=0.5) +
geom_sf(data=st_get_contour(skde1), colour=1, fill=NA, show.legend=FALSE) +
coord_sf(xlim=xlim, ylim=ylim) + ggthemes::theme_map()
To generate a filled contour plot, the only required changes are to
input an appropriate colour scale function, and for a base
R plot, to set legend=TRUE since its default
value is FALSE.
## R base filled contour plot
plot(st_geometry(wa), xlim=xlim, ylim=ylim)
plot(skde1, add=TRUE, pal=function(.) { sequential_hcl(n=., h1=275, rev=TRUE) })
mapsf::mf_scale(size=200, lwd=4)
## geom_sf filled contour
gs + geom_sf(data=st_get_contour(skde1), aes(fill=label_percent(contlabel))) +
scale_fill_discrete_sequential(h1=275) +
coord_sf(xlim=xlim, ylim=ylim) + ggthemes::theme_map()
Since the output from st_kde is compatible with
geom_sf, then it is easy to generate multiple maps of
related geospatial KDEs. For example, KDEs for each Grevillea species,
with probability contour levels or with fixed contour levels:
skde1g <- st_kde(grevillea_ep)
## facetted geom_sf filled contour
gsg <- ggplot(skde1g) + gwa + gsc
gsg + geom_sf(data=st_get_contour(skde1g), aes(fill=label_percent(contlabel))) +
scale_fill_discrete_sequential(h1=275) +
facet_wrap(~name) + coord_sf(xlim=xlim, ylim=ylim) + ggthemes::theme_map()
## facetted geom_sf filled contour with fixed contour levels for all facets
bkde1g <- contour_breaks(skde1g)
gsg + geom_sf(data=st_get_contour(skde1g, breaks=bkde1g), aes(fill=contlabel)) +
scale_fill_discrete_sequential(h1=275) + facet_wrap(~name) +
coord_sf(xlim=xlim, ylim=ylim) + ggthemes::theme_map()
References
Chacon, J. E. and Duong, T. (2018). Multivariate Kernel Smoothing and Its Applications. Chapman & Hall/CRC Press, Boca Raton.
Duong, T. (2022) Statistical visualisation for tidy and geospatial data in R via kernel smoothing methods in the eks package. https://doi.org/10.48550/arXiv.2203.01686