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
st_kde functions in the
eks package fills in these gaps.
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).
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
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.
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.
As an alternative to these discretised contours, the usual
ggplot2::geom_raster generates a plot with a continuous colour scale.
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
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, alpha=0.5) + 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.
The probability contour levels computed in
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.
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, alpha=0.1, colour=1)
Since geospatial data can be visualised with both base
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
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=alpha(3,0.5), pch=16) plot(st_geometry(paradoxa), add=TRUE, col=alpha(6,0.5), pch=17) mapsf::mf_legend(type="symb", val=c("Grevillea eryngioides", "Grevillea paradoxa"), pal=alpha(c(3,6), 0.5), pt_pch=16:17, pt_cex=c(1,1), title="Species", pos="topright") mapsf::mf_scale(size=200, lwd=4)
## geom_sf scatter plot theme_set(ggthemes::theme_map()) theme_update(legend.position=c(0.99,0.99), legend.justification=c(1,1)) gsc <- ggspatial::annotation_scale(data=data.frame(name="Grevillea paradoxa"), location="br", width_hint=0.2, bar_cols=1) ggplot() + gwa + gsc + geom_sf(data=grevillea_ep, aes(colour=name, shape=name), alpha=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
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
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
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:
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