boxplots.tex

% zip code.zip boxplots-*.r images/* lvplot/*
\documentclass[oneside]{article}
\usepackage{fullpage}
\usepackage{mathtools}
\usepackage[pdftex]{graphicx}
\DeclareGraphicsExtensions{.png,.pdf}
\graphicspath{{images/}}
\usepackage{amssymb}
\usepackage{hyperref}
\usepackage{verbatim}
\usepackage{booktabs}
\usepackage[round,sort]{natbib}
\bibliographystyle{abbrvnat}
\renewcommand\rmdefault{bch}
\usepackage[small]{caption}
\usepackage[small]{titlesec}
\linespread{2.07} 

\title{40 years of boxplots}
\author{Hadley Wickham and Lisa Stryjewski}

\raggedbottom

\pdfminorversion=4
\begin{document}
\maketitle 

\begin{abstract} 
  
  The boxplot plot has been around for over 40 years. This paper summarises the improvements, extensions and variations since Tukey first introduced his ``schematic plot'' in 1970. We focus particularly on richer displays of density and extensions to 2d.

\end{abstract}

\section{Introduction}

% Add reference to Schmid book 178.  Thank Glen Barnett,
% http://pubs.amstat.org/doi/pdf/10.1198/tas.2011.11012

John Tukey introduced the box and whiskers plot as part of his toolkit for exploratory data analysis \citep{tukey:1970}, but it did not become widely known until formal publication \citep{tukey:1977}. The boxplot is a compact distributional summary, displaying less detail than a histogram or kernel density, but also taking up less space. Boxplots use robust summary statistics that are always located at actual data points, are quickly computable (originally by hand), and have no tuning parameters. They are particularly useful for comparing distributions across groups.

Today, over 40 years later, the boxplot has become one of the most frequently used statistical graphics, and is one of the few plot types invented in the 20th century that has found widespread adoption. Due to their elegance and practicality, boxplots have spawned a wealth of variations and enhancement. This paper pulls these together in one place, showing how the boxplot has evolved.

We begin with a review of Tukey's definition and an overview of minor variations to both the underlying summary statistics and their visual representation. Section~\ref{sec:density} describes the richer displays of density facilitated by widespread desktop computing, and Section~\ref{sec:2d} explores how the boxplot has been extended to deal with 2d data. We conclude with some comments on the state of boxplot research and describe where future contributions are most needed.

The online supplementary materials include all R code \citep{R} used to create plots in this paper, and features original code for four boxplots (vase plot, quelplot, rotational boxplot, and bivariate clockwise boxplot) that previously lacked publicly available implementation.

\section{Tukey's boxplot}
\label{sec:tukey}

The basic graphic form of the boxplot, the range-bar, was established in the early 1950's \citet[pg. 164]{spear:1952}. Tukey's contribution was to think deeply about appropriate summary statistics that worked for a wide range of data and to connect those to the visual components of the range bar. Today, what we call a boxplot is more closely related to what Tukey called a schematic plot, a box and whiskers plot with some special restrictions on the summary statistics used. % Tukey suggested using dashed lines for the more precise schematic plot, and solid lines for the box-and-whisker plot.

The boxplot is made up of five components, carefully chosen to give a robust summary of the distribution of a dataset:

\begin{itemize}

\item the \emph{median},

\item two \emph{hinges}, the upper and lower fourths (quartiles),

\item the data values adjacent to the upper and lower \emph{fences}, which lie 1.5 times the inter-fourth range from the median,

\item two \emph{whiskers} that connect the hinges to the fences, and

\item (potential) \emph{out-liers}, individual points further away from the median than the extremes.

\end{itemize}

\noindent These elements are summarised in Figure~\ref{fig:construction}. Our notation follows Tukey's, except where we can be more precise or where common usage has changed over the last 40 years. 

\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.5]{components}
  \caption{Construction of a boxplot.  Labels on the left give names for graphic elements, labels on the right give the corresponding summary statistics.}
  \label{fig:construction}
\end{figure}

There are a number of variations of these basic definitions. As well as variations in the definition of a quantile \citep{hyndman:1996}, some boxplots replace the extremes with fixed quantiles (e.g. min and max, 2\% and 98\%) or use multipliers other than 1.5 for the whiskers \citep{frigge:1989}. Others use the semi-interquartile ranges (e.g.\ $Q_1 - Q_2$) for asymmetric whiskers \citep{rousseuw:1999}, explicit adjustments to the extremes to account for skewness \citep{hubert:2008}, alternative definitions of fences \citep{dumbgen:2007} or alternative definitions of outliers \citep{schwertman:2004,carter:2009}. Others have used additional graphical elements to display distributional features like kurtosis \citep{aslam:1991}, skewness and multimodality \citep{choonpradub:2005}, and mean and standard error \citep{marmolejo-ramos:2010}.

One of the appealing attributes of the boxplot is that if you have a rank function for the type of data you are dealing with, you can generate a boxplot. This makes it easy extend to the boxplot to work with weighted data, as described by \citet{survey,korn:1998} for survey weights, by \citet{willmott:2007} for spatial area weights, and by \citet{dykes:2007} for distance weights.

In an effort to improve the data-ink ratio of the boxplot, \citep{tufte:2001} proposed the midgap plot. As shown in Figure~\ref{fig:tufte}, the box is removed and the median line replaced with a dot. No information is lost, and the boxplot becomes substantially more compact. However, perceptual studies \citep{stock:1991} have found Tufte's variation to be substantially less accurate than the original. \citet{carr:1994a} proposed a colourful variation, also shown in Figure~\ref{fig:tufte}. This variation is designed to be tightly perceptually linked, so that each boxplot appears a single object, not a collection of lines. No perceptual testing has been performed on this variant.

\begin{figure}[htbp]
  \centering
  \includegraphics[scale = 0.5, angle = 270]{midgap}
  \caption{Tukey's original boxplot (top) compared to Tufte's box-less (middle) and Carr's colourful (bottom) variations. When colour is available, Carr suggests using red for components above the median and blue for colours below.}
  \label{fig:tufte}
\end{figure}

Another variation aims to overcome an important problem with the boxplot: there is visual display of group size, and hence no way of assessing if the differences are significant. The variable-width and notched boxplots \citep{mcgill:1978} add inferential detail. As the name suggests, the box widths of the variable-width boxplot vary according to the number of points in the group. The notched boxplot goes one step further by displaying confidence intervals around the medians, supporting visual assessment of statistical significance. The length of the confidence interval is determined heuristically so that non-overlapping intervals imply (approximately) a difference at the 5\% level, regardless of the underlying distribution.

\begin{figure}[htbp]
  \centering
  \includegraphics[width = 0.25\linewidth]{width-boxplot}% 
  \includegraphics[width = 0.25\linewidth]{width-variable}% 
  \includegraphics[width = 0.25\linewidth]{width-notched}
  \caption{Boxplot variations showing 100, 1000, 10000, and 100000 numbers drawn from a standard normal distribution.  (Left) In a regular boxplot the only hint that the groups are different sizes is the number of outliers. (Middle) A variable-width boxplot shows the differences in group size.  (Right) The notched boxplots displays an inferentially meaningful quantity: the error associated with the estimate of the median.}
  \label{fig:width}
\end{figure}

Other more unusual variations are an adaption for circular variables \citep{abuzaid:2011}, and an adaption to make boxplots more suitable for display as glyphs \citet{carr:1998}, particularly when overlaid on maps to display how data distribution varies in space.

There have been some perceptual studies on boxplots. \citet{behrens:1990} found evidence of significant bias when reading the length of the whiskers: whisker length was overestimated when whiskers were shorter than boxes and underestimated when whiskers were longer than boxes. There is a similar bias for reading the length of boxes: box length is overestimated when boxes are shorter than whiskers and vice-versa. Notched plots appear to suffer from similar problems \citep{wells:1996}.

\section{Richer displays of density}
\label{sec:density}

One of the original constraints on the boxplot was that it was designed to be computed and drawn by hand. As every statistician now has a computer on their desk, this constraint can be relaxed, allowing variations of the boxplot that are substantially more complex. These variations attempt to display more information about the distribution, maintaing the compact size of the boxplot, but bringing in the richer distributional summary of the histogram or density plot. These plots can overcome problems in the original such as the failure to display multi-modality, or the excessive number of ``outliers'' when $n$ is large.

The first variation to display a density estimate was the \emph{vase plot} \citep{benjamini:1988}, where the box is replaced with a symmetrical display of estimated density. \emph{Violin plots} \citep{hintze:1998} are very similar, but display the density for all data points, not just the middle half. The \emph{bean plot} \citep{kampstra:2008} is a further enhancement that adds a rug that showing every value and a line that shows the mean. The name is inspired by the appearance of the plot: the shape of the density looks like the outside of a bean pod, and the rug plot looks like the seeds within. \citet{kampstra:2008} also suggests a way of comparing two groups more easily: use the left and right sides of the bean to display different distributions. A related idea is the raindrop plot \citep{barrowman:2003}, but its focus is on the display of error distributions from complex models.

Figure~\ref{fig:density} demonstrates these density boxplots applied to 100 numbers drawn from each of four distributions with mean 0 and standard deviation 1: a standard normal, a skew-right distribution (Johnson distribution with skewness 2.2 and kurtosis 13), a leptikurtic distribution (Johnson distribution with skewness 0 and kurtosis 20) and a bimodal distribution (two normals with mean -0.95 and 0.95 and standard deviation 0.31). Richer displays of density make it much easier to see important variations in the distribution: multi-modality is particularly important, and yet completely invisible with the boxplot.

\begin{figure}[htbp]
  \centering
  \includegraphics[width = 0.25\linewidth]{four-box}%
  \includegraphics[width = 0.25\linewidth]{four-vase}%
  \includegraphics[width = 0.25\linewidth]{four-violin}%
  \includegraphics[width = 0.25\linewidth]{four-bean}
  \caption{From left to right: box plot, vase plot, violin plot and bean plot.  
Within each plot, the distributions from left to right are: standard normal (n), right-skewed (s), leptikurtic (k), and bimodal (mm). A normal kernel and bandwidth of 0.2 are used in all plots for all groups.}
  \label{fig:density}
\end{figure}

A more sophisticated display is the sectioned density plot \citep{cohen:2006}, which uses both colour and space to stack a density estimate into a smaller area, hopefully without losing any information (not formally verified with a perceptual study). The sectioned density plot is similar in spirit to horizon graphs for time series \citep{reijner:2008}, which have been found to be just as readable as regular line graphs despite taking up much less space \citep{heer:2009}. The density strips of \citet{jackson:2008} provide a similar compact display that uses colour instead of width to display density. These methods are shown in Figure~\ref{fig:density-display}. 

The summary plot \citep{potter:2010} is a similar idea. It combines a minimal boxplot with glyphs representing the first five moments (mean, standard deviation, skewness, kurtosis and tailings), and a sectioned density plot crossed with a violin plot (both colour and width are mapped to estimated density), and an overlay of a reference distribution. It is a rather busy display.

\begin{figure}[htbp]
  \centering
  \includegraphics[width = 0.25\linewidth]{four-sectioned}%
  \includegraphics[width = 0.25\linewidth]{four-denstrip}
  \caption{(Left) sectioned density plot and (right) density strips, same four distributions as Figure~\ref{fig:density}. A normal kernel and bandwidth of 0.2 are used in all plots for all groups.}
  \label{fig:density-display}
\end{figure}

The highest density region (HDR) boxplot \citep{hyndman:1996a} is a compromise between a boxplot and a density boxplot. It uses a density estimate but shows only two regions of highest density: the top 50\% and 99\%. These regions do not need to be contiguous and make it easy to spot multi-modality. The disadvantage of HDR boxplots is a less-sophisticated definition of extremes, making the outliers less useful for non-normal data. Figure~\ref{fig:hdr} shows the HDR boxplot for the four distributions previously described.

\begin{figure}[htbp]
  \centering
    \includegraphics[width = 0.25\linewidth]{four-hdr}%    
  \caption{The highest density region boxplot for the same four distributions as Figure~\ref{fig:density}. The multimodality in the fourth distribution is easy to spot.}
  \label{fig:hdr}
\end{figure}

Each author has suggested a different density estimate to use in conjunction with their new display, but there is no reason not to use any desired estimator. This is the price of density boxplots: the explosion of choices. Which density estimate should you use? Which choice of bandwidth or bin width is best? Bandwidth estimation is particularly challenging: if multiple groups are displayed, should each group get its own bandwidth, or should one bandwidth be used for all? \citet{kampstra:2008} suggests using the average of the per-group bandwidth estimates. The following two methods attempt a richer display of density without the cost of additional tuning parameters.

The box-percentile plot \citep{esty:2003} displays a modified empirical cumulative density function (ECDF). The width of each box is proportional to the percentile, up to the 50th percentile, after which the width is proportional to one minus the the percentile. Lines mark the median and upper and lower quartiles. While this display of the ECDF contains all information about the distribution, it is not always easy to parse this data into an informative mental model. This is illustrated in Figure~\ref{fig:bpp}: without training, it is very difficult to tell that the fourth distribution is bimodal.

Mountain plot: \citep{monti:1995}

\begin{figure}[htbp]
  \centering
  \includegraphics[width = 0.25\linewidth]{four-bpp}
  \caption{Box-percentile plots for the same four distributions used in Figure~\ref{fig:density}. The multimodality in the fourth distribution is hard to spot.}
  \label{fig:bpp}
\end{figure}

The letter-value boxplot \citep{hofmann:2006a} was designed to overcome the shortcomings of the boxplot for large data. For large datasets ($n \gtrsim 10,000$), the boxplot displays many outliers, and doesn't take advantage of the more reliable estimates of tail behaviour. The letter-value boxplot extends the boxplot with additional letter-values apart from the median (M) and fourths (F): eigths (E), sixteenths (D), ..., until the estimation error becomes too large. Each additional letter-value is displayed with a slightly smaller box, as shown in Figure~\ref{fig:letter-value}. Because this display just adds extra letter values, it suffers from the same problems as the original boxplot, and multimodality is almost impossible to spot.

\begin{figure}[htbp]
  \centering
  \includegraphics[width = 0.25\linewidth]{letter-value}
  \caption{Letter value plots of 100, 1000, 10,000, and 100,000 points drawn from a standard normal distribution.  The letter-value boxplot automatically displays additional letter values as the sample size increases.}
  \label{fig:letter-value}
\end{figure}

\section{Extensions to 2d}
\label{sec:2d}

Extending the boxplot to work in 2d is challenging because of the difficulty of defining order statistics, depth and quantiles in 2d. There is no unique definition of rank in 2d dimensions, and hence the 2d analogues of medians, fourths and extremes becomes more complex (mathematically and computationally). Perhaps due to the increased complexity of creating just a single plot, there has been little development of effective methods for comparing multiple groups. Compared to 1d, it is less obvious that 2d boxplots provide significant advantages over contours of density estimates or heatmaps of binned counts.

The first attempt at generalising the boxplot to 2d was the \emph{rangefinder plot} \citep{becketti:1987}. This generalisation is simple, treating the two variables as independent and drawing two independent 1d boxplots. The \emph{relplot} \citep{goldberg:1992} relaxes the assumption of independence by robustly fitting a bivariate Gaussian to the data, and drawing 50\% (corresponding to the box) and 99\% (corresponding to the whiskers) confidence ellipses. The \emph{quelplot}\citep{goldberg:1992} relaxes the assumption of normality by adding two degrees of asymmetry, accounting for residuals on both the major and minor axes of the ellipse. Another variation is the 2d boxplot of \citet{tongkumchum:2005}, which effectively constructs a 1d boxplot parallel to a (robust) line of best fit to the data.

The \emph{bagplot}, \citep{rousseuw:1999}, is a 2d analog of the boxplot with a \emph{bag} and \emph{fence} (cf. box) containing the middle 50\% of the data and a \emph{loop} (cf. whiskers) to separate outliers. The definitions of these regions follows naturally from the boxplot (although the numerical constants differ somewhat, adjusted by simulation to get behaviour that matches the boxplot), using the half-space depth in lieu of rank. The half-space depth, the smallest number of points contained within a half-plane of any orientation at a point, was informally defined early by \citet{tukey:1975}, but an efficient algorithm for computing it was not available until over 20 years later \citep{rousseeuw:1996}. A similar approach is the robust bivariate boxplot \citep{zani:1998}, which uses convex hull peeling to find the central regions, and then displays the hulls smoothed with a b-spline.

Figure~\ref{fig:2d} shows the range finder plot, relplot, quelplot and bagplot. The data is generated from a mixture of two highly correlated bivariate normals.

\begin{figure}[htbp]
  \centering
  \includegraphics[width = .25\linewidth]{2d-rangefinder}%
  \includegraphics[width = .25\linewidth]{2d-relplot}%
  \includegraphics[width = .25\linewidth]{2d-quelplot}%
  \includegraphics[width = .25\linewidth]{2d-bagplot}
  \caption{From left to right: a range finder plot, a relplot, a quelplot and a bagplot.}
  \label{fig:2d}
\end{figure}

The 1d HDR boxplot extends in a straightforward manner to 2d \citep{hyndman:1996a}. Unlike the methods described above, this 2d boxplot relies on a density estimate, and hence the selection of bandwidth is critically important. Figure~\ref{fig:2d-hdr}, shows the 2d HDR region boxplots with the 50\% and 95\% highest density regions, for three different bandwidths.

\begin{figure}[htbp]
  \centering
  \includegraphics[width = .25\linewidth]{2d-hdr-5}%
  \includegraphics[width = .25\linewidth]{2d-hdr-2-5}%
  \includegraphics[width = .25\linewidth]{2d-hdr-1}
  \caption{The 2d HDR boxplot, with both bandwidths set to (from left to right), 5, 2.5 and 1. Choice of bandwidth makes a critical difference in the appearance of the plot.}
  \label{fig:2d-hdr}
\end{figure}

Two other approaches take a circular approach. The rotational boxplot of \citet{muth:2000} bins the data into multiple overlapping circular sections around the centroid then computes the boxplot summary statistics for each section. The whiskers, hinges and medians are connected with lines. The clockwise bivariate boxplots of \citet{corbellini:2002} are a similar idea, but based on projections rather than sections: the full dataset is projected onto multiple lines through the origin and the boxplot statistics computed for each. Both plots are shown in Figure~\ref{fig:2d-circular}. These graphics must be read carefully because while they look similar to topological maps, but there is no guarantee that density is highest in the central region. They are most useful for displaying the distribution of distances from a location.

\begin{figure}[htbp]
  \centering
  \includegraphics[width = .25\linewidth]{2d-clockwise}%
  \includegraphics[width = .25\linewidth]{2d-rotational}%

  \caption{(Left) The clockwise bivariate boxplot displays hinges and whiskers, and (right) the rotational boxplot displays hinges, median (thicker) and whiskers (in grey). The rotational boxplot has a bin width of $90\,^{\circ}$. }

  \label{fig:2d-circular}
\end{figure}

Functional boxplots \citep{hyndman:2010,sun:2011} take the ideas of 2d boxplots and extend them to the infinite-dimensional functional case. The extensions are similar in spirit to the bagplot, first finding a definition of rank in a functional space, and then extending the definition of the boxplot to use that definition, and finally finding an effective display technique.

\section{Discussion}
\label{sec:discussion}

Boxplots were created to provide a succinct distributional summary that could easily be created by hand, and supported comparison across groups. As computers have become more prevalent and more powerful, it has become easier  to produce compact summaries that display more data. This has led to an explosion of boxplot variations that stay true to the original goals to various extents, while supporting much richer display of the underlying distributions.

There are a few places where the existing literature is weak. We wonder if the proposed boxplot variations do actually make it easier to compare distributions. Cleveland's hierarchy \citep{cleveland:1984} provides some guidance suggesting that techniques that use position, rather than colour, should be more effective, but there has been little empirical testing. Additionally, despite the number of 2d boxplots that have been developed, none attack one of the original motivating problems of boxplots: comparing the distributions of multiple groups.

Boxplots have a rich history, and a promising future. We hope that this paper has illuminated the past, and suggested new avenues for the future.

\section{Acknowledgements}
\label{sec:acknowledgements}

One of the biggest challenges of this paper was recreating all of the plots with consistent formatting. We're indebted to the authors who have provided implementations in R packages: {\tt hdrcde} \citep{hdrcde}, {\tt UsingR} \citep{UsingR}, {\tt beanplot} \citep{kampstra:2008}, {\tt denstrip} \citep{jackson:2008}, {\tt aplpack} \citep{aplpack} and {\tt lvplot}. Code for the relplot came from the website of Rand Wilcox. We wrote our own implementations of the vase plot and the quelplot (but only using non-robust estimators). Tufte's mid gap plot and Carr's colourful boxplot were drawn by hand.

% bibtool -x boxplots.aux -c > references.bib 
\bibliography{references}

\end{document}