2  Data & Visualization

This chapter classifies the kinds of data you will work with and gives you the basic tools to visualize and summarize a single variable. We start with data types (cardinal vs factor), then look at the structures that hold datasets in R (lists and data.frames), and finish with the standard ways to depict a distribution: proportions, histograms, ECDFs, quantiles, and boxplots. Which depiction fits depends on the data type: proportions suit discrete data, histograms suit continuous data, and the ECDF, quantiles, and boxplot work for either.

2.1 Data Types

Before any analysis, it helps to know whether the values in our dataset support arithmetic (their differences carry consistent meaning) or whether they are merely labels.

ImportantKey Definition

Cardinal (or numeric) data are values where the difference between elements always means the same thing.

Factor data are values where the difference between elements does not always mean the same thing.

The type of a variable governs what summaries are meaningful: cardinal data can be averaged, scattered, and used in regressions, while factor data is best counted and tabulated. Trying to compute the mean of factor data is a category error rather than a meaningful summary, and most statistical software will refuse to do so.

Basic Types

Cardinal data come in two flavors (discrete or continuous) and factor data come in two flavors (ordered or unordered), giving four sub-types in all:

  • Cardinal (Numeric)
    • Discrete: E.g. \(\{ 1,2,3\}\) and notice that \(2-1=3-2\).
    • Continuous: E.g., \(\{1.4348, 2.4348, 2.9, 3.9 \}\) and notice that \(2.9-1.4348=3.9-2.4348\)
  • Factor
    • Ordered: E.g., \(\{1^{st}, 2^{nd}, 3^{rd}\}\) place in a race and notice that \(1^{st}\) - \(2^{nd}\) place does not equal \(2^{nd}\) - \(3^{rd}\) place for a very competitive person who cares only about winning.
    • Unordered (categorical): E.g., \(\{Amanda, Bert, Charlie\}\) and notice that \(Amanda - Bert\) never makes sense.

Note that for theoretical analysis, the types are sometimes grouped differently as

  • Discrete (discrete cardinal, ordered factor, and unordered factor data). You can count the potential values. E.g., the set \(\{A,B,C\}\) has three potential values.
  • Continuous (continuous cardinal data). There are uncountably infinite potential values. E.g., Try counting the numbers between \(0\) and \(1\) including decimal points, and notice that any two numbers have another potential number between them.

Here are some examples

Code 
dat_card1 <- c(1, 2, 3) # Cardinal data (Discrete)
dat_card1
## [1] 1 2 3

dat_card2 <- c(1.1, 2/3, 3) # Cardinal data (Continuous)
dat_card2
## [1] 1.1000000 0.6666667 3.0000000

dat_fact1 <- factor( c('A', 'B', 'C'), ordered=TRUE) # Factor data (Ordinal)
dat_fact1
## [1] A B C
## Levels: A < B < C

dat_fact2 <- factor( c('Leipzig', 'Los Angeles', 'Logan'), ordered=FALSE) # Factor data (Categorical)
dat_fact2
## [1] Leipzig     Los Angeles Logan      
## Levels: Leipzig Logan Los Angeles

dat_fact3 <- factor( c(TRUE, FALSE), ordered=FALSE) # Factor data (Categorical)
dat_fact3
## [1] TRUE  FALSE
## Levels: FALSE TRUE

# Explicitly check the data types:
#class(dat_card1)
#class(dat_card2)

Strings

Note that R allows for unstructured plain text, called character strings, which we can then format as factors

Code 
c('A', 'B', 'C')  # character strings
## [1] "A" "B" "C"
c('Leipzig', 'Los Angeles', 'Logan')  # character strings
## [1] "Leipzig"     "Los Angeles" "Logan"

Also note that strings are encounter in a variety of settings, and you often have to format them after reading them into R (the RStudio regex cheatsheet is a handy reference).1

Code 
# Strings
paste( 'hi', 'mom')
## [1] "hi mom"
paste( c('hi', 'mom'), collapse='--')
## [1] "hi--mom"

kingText <- 'The king infringes the law on playing curling.'
gsub(pattern='ing', replacement='', kingText)
## [1] "The k infres the law on play curl."
# advanced usage
#gsub('[aeiouy]', '_', kingText)
#gsub('([[:alpha:]]{3})ing\\b', '\\1', kingText)

2.2 Datasets

Datasets can be stored in a variety of formats on your computer. But they can be analyzed in R in three basic ways.

Lists

Lists are probably the most basic type

Code 
x0_hat <- seq(1, 10)                       # vector of 1..10
y_hat <- 2*x0_hat                              # element-wise doubling
list(x0_hat, y_hat)                            # list holds two vectors of equal length
## [[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## [[2]]
##  [1]  2  4  6  8 10 12 14 16 18 20

x_mat1 <- matrix( seq(2, 7), 2, 3)    # 2x3 matrix
x_mat2 <- matrix( seq(4, -1), 2, 3)
list(x_mat1, x_mat2)                  # list can also hold matrices
## [[1]]
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7
## 
## [[2]]
##      [,1] [,2] [,3]
## [1,]    4    2    0
## [2,]    3    1   -1

Lists are useful for storing unstructured data

Code 
list(list(x_mat1), list(x_mat2))  # list of lists
## [[1]]
## [[1]][[1]]
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7
## 
## 
## [[2]]
## [[2]][[1]]
##      [,1] [,2] [,3]
## [1,]    4    2    0
## [2,]    3    1   -1

list(x_mat1, list(x_mat1, x_mat2)) # list of different objects
## [[1]]
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7
## 
## [[2]]
## [[2]][[1]]
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7
## 
## [[2]][[2]]
##      [,1] [,2] [,3]
## [1,]    4    2    0
## [2,]    3    1   -1

# ...inception...
list(x_mat1,
    list(x_mat1, x_mat2), 
    list(x_mat1, list(x_mat2)
    )) 
## [[1]]
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7
## 
## [[2]]
## [[2]][[1]]
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7
## 
## [[2]][[2]]
##      [,1] [,2] [,3]
## [1,]    4    2    0
## [2,]    3    1   -1
## 
## 
## [[3]]
## [[3]][[1]]
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7
## 
## [[3]][[2]]
## [[3]][[2]][[1]]
##      [,1] [,2] [,3]
## [1,]    4    2    0
## [2,]    3    1   -1

Data.frames

A data.frame looks like a matrix but each column can be a different type of vector (or even list). This allows you to combine different data types into a single object for analysis, which is why it might be your most common object.

Code 
# data.frames: your most common data type
    # matrix of different data-types
    # well-ordered lists
data.frame(x=x0_hat, y=y_hat)  # list of vectors
##     x  y
## 1   1  2
## 2   2  4
## 3   3  6
## 4   4  8
## 5   5 10
## 6   6 12
## 7   7 14
## 8   8 16
## 9   9 18
## 10 10 20

Create a data.frame storing two different types of data. Then show print only the second column

Code 
dat <- data.frame(x=dat_fact2, y=dat_card2)
dat
##             x         y
## 1     Leipzig 1.1000000
## 2 Los Angeles 0.6666667
## 3       Logan 3.0000000

dat[, 'y']
## [1] 1.1000000 0.6666667 3.0000000

Data.frames stack two dimensions, rows and columns. The generalization to three or more dimensions is the array.

Initial Data Inspection

Regardless of the data types you have, you typically begin by inspecting your data by examining the first few observations.

Consider, for example, historical data on crime in the US.

Code 
head(USArrests) # Actual Data
##            Murder Assault UrbanPop Rape
## Alabama      13.2     236       58 21.2
## Alaska       10.0     263       48 44.5
## Arizona       8.1     294       80 31.0
## Arkansas      8.8     190       50 19.5
## California    9.0     276       91 40.6
## Colorado      7.9     204       78 38.7

# Check NA values
x0_na <- c(3, 3.1, NA, 0.02) # Small dataset with a missing value
sum( is.na(x0_na) )
## [1] 1

Packages

R already comes with many datasets, such as USArrests. They live in a package called datasets, which R loads every time it starts. A package is a bundle of functions and datasets that someone has written and shared, so you can use them without writing them yourself.

Many more datasets live in packages that are not part of R when you first install it. You install such a package once, and then load its data in each script that uses it. For example, the Ecdat package contains economics datasets, such as Wages1 on the hourly wages of \(3294\) workers.

Code 
install.packages('Ecdat') # only once, typed in the console
Code 
data(Wages1, package='Ecdat') # every session, saved in your script
head(Wages1)
##   exper    sex school     wage
## 1     9 female     13 6.315296
## 2    12 female     12 5.479770
## 3    11 female     11 3.642170
## 4     9 female     14 4.593337
## 5     8 female     14 2.418157
## 6     9 female     14 2.094058

It is best to keep install.packages() out of your script, so that re-running the script does not reinstall the package every time. To list every dataset in a package, use data(package='Ecdat'). To read about one dataset, use ?Ecdat::Wages1.

Restart R and run head(Wages1) before loading the data. What happens?

R reports that it cannot find Wages1, because a new session starts with only the packages R loads by default. The package is still installed, so running data(Wages1, package='Ecdat') again fixes this without reinstalling anything.

2.3 Proportions

In what follows, we will often work with data as a vector, where there are \(n\) observations and \(\hat{x}_{i}\) is the value of the \(i\)th one. The hat marks a value computed from observed data. We often analyze observations in comparison to some value \(x\).

Code 
x0_hat <- c(3, 3.1, 0.02) # Small dataset we will use in numerical examples
x <- 2 # particular value: 'little x'
sum(x0_hat <= x)
## [1] 1
sum(x0_hat == x)
## [1] 0

Empirical Mass Function

If you have factor data, or discrete cardinal data, you can directly plot the proportions. For each unique outcome \(x\) we compute \[ \hat{p}(x)=\sum_{i=1}^{n}\mathbf{1}\left(\hat{x}_{i}=x\right)/n. \] where \(n\) is the number of observations and \(\sum_{i=1}^{n}\mathbf{1}\left(\hat{x}_{i}=x\right)\) counts the number of observations equal to \(x\).

For example, consider the dataset \(\{1,2,1,3,1\}\) with \(n=5\) observations. The value \(x=1\) appears three times, so \(\hat{p}(1)=3/5\). The value \(x=2\) appears once, so \(\hat{p}(2)=1/5\). The value \(x=3\) appears once, so \(\hat{p}(3)=1/5\). The proportions sum to one: \(3/5 + 1/5 + 1/5 = 1\).

Code 
x0_hat <- c(1, 2, 1, 3, 1)
table(x0_hat)
## x0_hat
## 1 2 3 
## 3 1 1
table(x0_hat)/length(x0_hat)
## x0_hat
##   1   2   3 
## 0.6 0.2 0.2
Code 
# Discretized data
x_hat <- USArrests[, 'Murder'] # lab data: murder arrests per 100,000
x_hat[1:3]
## [1] 13.2 10.0  8.1
x_rounded <- floor(x_hat) #rounded down, to be discrete
x_rounded[1:3]
## [1] 13 10  8

p_hat <- table(x_rounded)/length(x_rounded)
plot(p_hat, col=grey(0, .5),
    xlab='Murder Arrests (Discretized)',
    ylab='Proportion of States with each value')

There are several stylistic variations: a bar plot uses thick bars for each outcome, a lollipop plot above uses skinny lines with a circle on top. E.g., add points(names(p_hat), p_hat, pch=16) The most important point is that the height of each line equals the proportion of data with a specific value.

Histogram Density

For continuous data, no single value is likely to repeat, so we summarize the distribution by counting how many observations land in each of several equal-width intervals.

ImportantKey Definition

A histogram density divides the range of the data into bins of equal half-width \(h\). For an exclusive bin \((x-h, x+h]\) with midpoint \(x\), it reports the rescaled count \[\hat{f}(x) = \frac{\sum_{i=1}^{n} \mathbf{1}\left( \hat{x}_{i}~\text{in}~(x-h, x+h] \right)}{n \cdot 2h},\] where \(\mathbf{1}(\cdot)\) is an indicator function that equals \(1\) if the expression inside is TRUE and \(0\) otherwise.

The histogram is the workhorse picture for the shape of a continuous distribution: where observations are concentrated, whether there are multiple peaks, whether one tail is longer than the other, and whether there are gaps with no observations. These features are hard to spot in a list of numbers but jump out at a glance from a histogram. This rescaling makes the total area of all bars equal to one, so bar heights have units of probability per unit of \(x\). It also makes histogram heights comparable across samples of different sizes and lets us interpret the bar height as a density: the proportion of data per unit of \(x\). The half-width \(h\) acts as a bandwidth: a small \(h\) produces a jagged picture sensitive to single observations, while a large \(h\) smooths over real features. Note that \(h\) is measured from the bin midpoint out to one edge, so a bin is \(2h\) wide.

E.g., suppose \(\hat{x}_{i}=3.8\) and \(h=1/2\).

For \(x=1\) we have \(\mathbf{1}\left( \hat{x}_{i}~\text{in}~\left(1-h, 1+h \right] \right)=\mathbf{1}\left( 3.8~\text{in}~\left(0.5, 1.5\right] \right)=0\).

For \(x=4\) \(\mathbf{1}\left( \hat{x}_{i}~\text{in}~\left(4-h, 4+h \right] \right)=\mathbf{1}\left( 3.8~\text{in}~\left(3.5, 4.5\right] \right)=1\).

Note that rectangle area equals the proportion of data in the bin. Recalling that the area of the rectangle is “base x height”: \[\begin{aligned} 2h \times \hat{f}(x) &= \frac{\sum_{i=1}^{n} \mathbf{1}\left( \hat{x}_{i}~\text{in}~\left(x-h, x+h \right] \right)}{n}. \end{aligned}\] We compute \(\hat{f}(x)\) for each bin midpoint \(x\), and the area of all rectangles sums to one. 2

For example, consider the dataset \(\{3,3.1,0.02\}\) and use bins \((0,1], (1,2], (2,3], (3,4]\). In this case, the midpoints are \(x=(0.5,1.5,2.5,3.5)\) and \(h=1/2\). Then the counts at each midpoints are \((1,0,1,1)\). Since \(\frac{1}{n 2h}=\frac{1}{3\times1}=\frac{1}{3}\), we can rescale the counts to compute the density as \(\hat{f}(x)=(1,0,1,1) \frac{1}{3}=(1/3,0,1/3,1/3)\).

Code 
# Intuitive Examples
x0_hat <- c(3, 3.1, 0.02)
f_hat <- hist(x0_hat, breaks=c(0, 1, 2, 3, 4), plot=FALSE, freq=FALSE)
f_hat
## $breaks
## [1] 0 1 2 3 4
## 
## $counts
## [1] 1 0 1 1
## 
## $density
## [1] 0.3333333 0.0000000 0.3333333 0.3333333
## 
## $mids
## [1] 0.5 1.5 2.5 3.5
## 
## $xname
## [1] "x0_hat"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

# base x height
base <- 1
height <- f_hat[['density']]
sum(base*height)
## [1] 1

For another example, use the bins \((0,2]\) and \((2,4]\). So the midpoints are \(1\) and \(3\), and the bin half-width is \(1\). Only one observation, \(0.02\), falls in the bin \((0,2]\). The other two observations, \(3\) and \(3.1\), fall into the bin \((2,4]\). The scaling factor is \(\frac{1}{n 2h}=\frac{1}{3\times 2 \times 1}=\frac{1}{6}\). So the first bin has density \(\hat{f}(1)=1\frac{1}{6}=1/6\) and the second bin has density \(\hat{f}(3)=2\frac{1}{6}=2/6\). The area of the first bin’s rectangle is \(2 \times \hat{f}(1)=2/6=1/3\) and the area of the second rectangle is \(2 \times \hat{f}(3)=4/6=2/3\).

Now intuitively work through an example with three bins instead of four. Compute the areas

Code 
f_hat <- hist(x0_hat, breaks=c(0, 4/3, 8/3, 4), plot=FALSE, freq=FALSE)
base <- 4/3
height <- f_hat[['density']]
sum(base*height)
## [1] 1

# as a default, R uses bins (, ] instead of [, )
# but you can change that with 'right=F'
# hist(x0_hat, breaks=c(0, 4/3, 8/3, 4), plot=F, right=F)
Code 
# Practical Example
hist(x_hat,
    breaks=seq(0, 20, by=1), # bin width =1, so the half-width is h=1/2
    freq=FALSE,
    border=NA,
    main=NA,
    xlab='Murder Arrests')
title('Murder Arrests', font.main=1)
# Raw Observations
rug(x_hat, col=grey(0, .5))

Code 

# Since 2h=1, the density equals the proportion of states in each bin
# Redo this example with h=1 (i.e., breaks spaced 2h=2 apart)

Note that histograms qualitatively depict proportions over multiple bins (since each bar represents the proportion of data in its bins), but you can also directly compute proportions.

E.g., the proportion of murder rates in \((2,4)\) is

Code 
mean( x_hat > 2 & x_hat <4)
## [1] 0.22

2.4 Distributions

Empirical Cumulative Distribution Function

Another way to summarize a distribution asks, for each candidate value \(x\), what proportion of the sample lies at or below it. It is a step function (see this ECDF walkthrough for visual examples) and contains exactly the same information as the sorted dataset, just packaged as a function of \(x\) rather than a list of values.

ImportantKey Definition

The empirical cumulative distribution function (ECDF) is the proportion of observations in the sample whose values are less than or equal to \(x\), \[\hat{F}(x) = \frac{1}{n}\sum_{i=1}^{n}\mathbf{1}(\hat{x}_{i} \leq x).\] It jumps up by \(1/n\) at each data point and equals \(1\) once \(x\) passes the largest observation.

For example, consider the dataset \(\{3,3.1,0.02 \}\). We reorder the observations as \(\{0.02, 3, 3.1 \}\), so that there are discrete jumps of \(1/n=1/3\) at each value. Consider the points \(x~\text{in}~\{0.5,2.5,3.5\}\). At \(x=0.5\), \(\hat{F}(0.5)\) measures the proportion of the data \(\leq 0.5\). Since only one observation, \(0.02\), of three is \(\leq 0.5\), we can compute \(\hat{F}(0.5)=1/3\). Similarly, since only one observation, \(0.02\), of three is \(\leq 2.5\), we can compute \(\hat{F}(2.5)=1/3\). Since all observations are \(\leq 3.5\), we can compute \(\hat{F}(3.5)=1\).

Code 
F_hat <- ecdf(x0_hat)

# Visualized
plot(F_hat)

Code 

# Evaluated at each data point
x <- x0_hat
F_hat(x0_hat)
## [1] 0.6666667 1.0000000 0.3333333
#sum(x0_hat<=3)/length(x0_hat)
#sum(x0_hat<=3.1)/length(x0_hat)
#sum(x0_hat<=0.02)/length(x0_hat)

# Evaluated at other points
x <- c(0.5, 2.5, 3.5)
F_hat(x)
## [1] 0.3333333 0.3333333 1.0000000
#sum(x0_hat<=0.5)/length(x0_hat)
#sum(x0_hat<=2.5)/length(x0_hat)
#sum(x0_hat<=3.5)/length(x0_hat)

# Evaluated at the histogram bin edges
x <- c(0, 1, 2, 3, 4)
F_hat(x)
## [1] 0.0000000 0.3333333 0.3333333 0.6666667 1.0000000
diff(F_hat(x)) # proportion of data in each bin, as in the histogram
## [1] 0.3333333 0.0000000 0.3333333 0.3333333
Code 
F_hat <- ecdf(x_hat)
# proportion of murders <= 10
F_hat(10)
## [1] 0.7
# proportion of murders <= x, for all x
plot(F_hat, main=NA,
    xlab='Murder Arrests (x)',
    ylab='Proportion of States with Murder Arrests <= x',
    pch=16, col=grey(0, .5))
title('Murder Arrests ECDF', font.main=1)
rug(x_hat)

The ECDF is useful when we want precise numbers rather than a picture of shape: “what fraction of states have at most \(10\) murders per \(100\)k?” reads off directly without choosing bin widths or smoothing parameters. The ECDF also applies to more data types than the two previous depictions. The mass function \(\hat{p}\) needs values that repeat, so it suits discrete data, and the histogram \(\hat{f}\) needs a bin width, so it suits continuous data. The ECDF needs neither, since its only operation is checking whether \(\hat{x}_{i} \leq x\), so it works for discrete and continuous data alike. The ECDF does not work for unordered factors, where “\(\leq\)” has no meaning, but it does work for ordered factors such as letter grades.

Still, the ECDF and the histogram are two views of the same counts. The area of a histogram bar is the proportion of data in the bin \((x-h, x+h]\), which is also the rise in the ECDF across that bin. \[\begin{aligned} 2h \times \hat{f}(x) &= \hat{F}(x+h) - \hat{F}(x-h). \end{aligned}\] Conversely, adding up bar areas from left to right recovers \(\hat{F}\) at each bin edge.3 For a continuous example, the shaded bars on the histogram cover the states with at most \(5\) murder arrests, and their total area equals the height of the colored point on the ECDF, \(\hat{F}(5)\).

Code 
F_hat <- ecdf(x_hat)
breaks <- seq(0, 20, by=1) # bin width 2h=1
f_hat <- hist(x_hat, breaks=breaks, plot=FALSE)

# Same data, stacked: the shaded area on top equals the height of the point below
x <- 5 # particular value
below <- f_hat[['mids']] < x # bins at or below x
col_x <- rgb(1, 0, 0, .5)
par(mfrow=c(2, 1), mar=c(4, 4, 1, 1))

# Numerically verify shaded area equals point height
sum(f_hat[['density']][below]*diff(breaks)[below])
## [1] 0.32
F_hat(x)
## [1] 0.32

# Histogram: shade the bins at or below x
hist(x_hat, breaks=breaks, freq=FALSE, border=NA,
    col=ifelse(below, col_x, grey(.8)),
    main=NA, xlab=NA, ylab='Density')

# ECDF: mark the point (x, F_hat(x))
plot(F_hat, xlim=range(breaks), pch=16, col=grey(0, .5),
    main=NA, xlab='Murder Arrests', ylab='Proportion <= x')
segments(x, F_hat(x), -10, F_hat(x), col=col_x)
segments(x, F_hat(x), x, -1, col=col_x)
points(x, F_hat(x), pch=16, col=rgb(1, 0, 0, .8), cex=1.5)

For a concrete example, consider the dataset \(\{3,3.1,0.02 \}\). At the bin edges \(x~\text{in}~\{0,1,2,3,4\}\), the ECDF is \(\hat{F}=(0,1/3,1/3,2/3,1)\), so it rises by \((1/3, 0, 1/3, 1/3)\) across the bins \((0,1], (1,2], (2,3], (3,4]\). Those rises are exactly the bin proportions that the histogram rescaled by \(1/(2h)\) to get \(\hat{f}(x)=(1/3,0,1/3,1/3)\). The relationship between \(\hat{f}\) and \(\hat{F}\) is easy to check with real data.

Code 
# Histogram density from the ECDF, and back again
F_hat <- ecdf(x_hat)
breaks <- seq(0, 20, by=1) # bin width 2h=1
f_hat <- hist(x_hat, breaks=breaks, plot=FALSE)

# Rise in F_hat across each bin, divided by the bin width
f_from_F <- diff(F_hat(breaks))/diff(breaks)
all.equal(f_from_F, f_hat[['density']])
## [1] TRUE

# Bar areas added up from the left, evaluated at the upper bin edges
F_from_f <- cumsum(f_hat[['density']]*diff(breaks))
all.equal(F_from_f, F_hat(breaks[-1]))
## [1] TRUE

To see this generally, start from the definition of \(\hat{f}\) and write the indicator for the bin as the difference of two “at or below” indicators: \[\begin{aligned} 2h \times \hat{f}(x) &= 2h \times \frac{\sum_{i=1}^{n} \mathbf{1}\left( \hat{x}_{i}~\text{in}~(x-h, x+h] \right)}{n \cdot 2h} \\ &= \frac{1}{n}\sum_{i=1}^{n} \mathbf{1}\left( x-h < \hat{x}_{i} \leq x+h \right) \\ &= \frac{1}{n}\sum_{i=1}^{n} \left[ \mathbf{1}\left( \hat{x}_{i} \leq x+h \right) - \mathbf{1}\left( \hat{x}_{i} \leq x-h \right) \right] \\ &= \frac{1}{n}\sum_{i=1}^{n} \mathbf{1}\left( \hat{x}_{i} \leq x+h \right) - \frac{1}{n}\sum_{i=1}^{n} \mathbf{1}\left( \hat{x}_{i} \leq x-h \right) \\ &= \hat{F}(x+h) - \hat{F}(x-h). \end{aligned}\] The step that splits the indicator holds observation by observation: an \(\hat{x}_{i}\) at or below \(x-h\) contributes \(1-1=0\), one inside the bin contributes \(1-0=1\), and one above \(x+h\) contributes \(0-0=0\). Dividing by the bin width, \(\hat{f}(x)\) is the slope of \(\hat{F}\) across the bin (“rise over run”): where \(\hat{F}\) climbs steeply the histogram is tall, and where \(\hat{F}\) is flat the histogram is near zero.

Quantiles

Reading the ECDF in reverse asks: at what value \(x\) does a given proportion \(p\) of the data lie below?

ImportantKey Definition

The \(p\)th quantile of a sample is the value \(x\) where \(p\) percent of the data fall below \(x\) and \((1-p)\) percent fall above. The empirical quantile function \(\hat{q}(p) = \hat{F}^{-1}(p)\) inverts the ECDF, returning the quantile at any probability \(p~\text{in}~[0,1]\).

Quantiles are useful for precise numerical anchors: “the top \(10\%\) of cities”, “half of households earn less than…”, or “the middle \(50\%\) of test scores fall between \(x_1\) and \(x_2\)”. They also describe spread (range, \(IQR\)) and shape (the gap between median and mean) without committing to a smooth distributional model.

Five quantiles get their own names because they summarize the distribution at a glance:

  • The min (the \(p=0\) quantile) is the smallest value, where \(0\%\) of the data have lower values.
  • The lower quartile (the \(p=0.25\) quantile) is where \(25\%\) of the data have lower values.
  • The median (the \(p=0.5\) quantile) is the middle value, where \(50\%\) of the data have lower values.
  • The upper quartile (the \(p=0.75\) quantile) is where \(75\%\) of the data have lower values.
  • The max (the \(p=1\) quantile) is the largest value, where \(100\%\) of the data have lower values.

We also speak of deciles at \(p=0.1, 0.2, \ldots, 0.9\) when we want a finer split.

For example, that dataset \(\{0,0,0.02,3,5\}\) has a median of \(0.02\), with the lower quartile \(0\) and the upper quartile \(3\). (The number \(0\) is also special: the most frequent observation is called the mode.)

Code 
quantile(x0_hat, probs=c(0, .5, 1))
##   0%  50% 100% 
## 0.02 3.00 3.10

Now work through an intuitive example with observations \(\{1,2,...,13\}\). Hint: split the ordered observations into four groups.

Code 
# common quantiles
quantile(x_hat)
##     0%    25%    50%    75%   100% 
##  0.800  4.075  7.250 11.250 17.400

# All deciles are quantiles
quantile(x_hat, probs=seq(0, 1, by=.1))
##    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
##  0.80  2.56  3.38  4.75  6.00  7.25  8.62 10.12 12.12 13.32 17.40

# Visualized: Inverting the Empirical Distribution
F_hat <- ecdf(x_hat)
plot(F_hat, lwd=2, xlim=c(0, 20),
    pch=16, col=grey(0, .5), main=NA)
title('Inverting the ECDF', font.main=1)
# Two Examples of Quantiles
p <- c(.25, .9) # Lower Quartile, Upper Decile
cols <- c(rgb(1, 0, 0, .8), rgb(0, 0, 1, .8))
q_hat <- quantile(x_hat, p, type=1)
q_hat
##  25%  90% 
##  4.0 13.2
segments(q_hat, p, -10, p, col=cols)
segments(q_hat, p, q_hat, 0, col=cols)
mtext( round(q_hat, 2), 1, at=q_hat, col=cols)

There are some issues with quantiles with smaller datasets. E.g., to compute the median of \(\{3,3.1,0,1\}\), we need some ways to break ties (of which there are many options). The most common way is perhaps a simple average of the tied values: \((3+1)/2=2\).

To calculate quantiles, the computer defaults to sorting the observations from smallest to largest as \(\hat{x}_{(1)}, \hat{x}_{(2)},... \hat{x}_{(n)}\), and then computes quantiles as \(\hat{x}_{ (p \times n) }\), where \(p~\text{in}~[0,1]\) refers to probability and \((p\times n)\) is rounded. There are different ways to break ties and the min (\(p=0\)) is a special case.

Code 
x_sorted_hat <- sort(x_hat)
x_sorted_hat
##  [1]  0.8  2.1  2.1  2.2  2.2  2.6  2.6  2.7  3.2  3.3  3.4  3.8  4.0  4.3  4.4
## [16]  4.9  5.3  5.7  5.9  6.0  6.0  6.3  6.6  6.8  7.2  7.3  7.4  7.9  8.1  8.5
## [31]  8.8  9.0  9.0  9.7 10.0 10.4 11.1 11.3 11.4 12.1 12.2 12.7 13.0 13.2 13.2
## [46] 14.4 15.4 15.4 16.1 17.4

# median
x_sorted_hat[ceiling(length(x_sorted_hat)*.5)]
## [1] 7.2
quantile(x_hat, probs=.5, type=1) #tie break rule #1
## 50% 
## 7.2

## arbitrary quantiles
p <- 0.11
x_sorted_hat[ceiling(length(x_sorted_hat)*p)]
## [1] 2.6
quantile(x_hat, probs=p, type=1)
## 11% 
## 2.6

# min is a special case
min(x_sorted_hat)
## [1] 0.8
x_sorted_hat[ceiling(length(x_sorted_hat)*0)]
## numeric(0)
quantile(x_hat, probs=0, type=1)
##  0% 
## 0.8
x_sorted_hat[1]
## [1] 0.8

Boxplots

Boxplots also summarize the distribution. The boxplot shows the median (solid black line) and interquartile range (\(IQR=\) upper quartile \(-\) lower quartile; filled box).4 As a default, whiskers are shown as \(1.5\times IQR\) and values beyond that are highlighted as outliers (so whiskers do not typically show the data range). No state sits more than \(1.5\times IQR\) beyond a quartile, so the boxplot below draws no outlier circles. Like the ECDF, the boxplot works for discrete and continuous data alike, since it is built from quantiles rather than from bins or repeated values. It does need cardinal data: the whiskers extend \(1.5\times IQR\) beyond the box, and that distance is only meaningful when differences are.

Code 
# Default: median, box (IQR), whiskers at 1.5*IQR, outliers as circles
boxplot(x_hat,
    main=NA, ylab='Murder Arrests')
title('Default', font.main=1)

The default boxplot compresses \(n\) observations into five numbers, which discards information about shape. You can alternatively show all the raw data points instead of whisker+outliers. This exposes the gaps, clusters, and repeated values that a five-number summary hides (Weissgerber et al. 2015).

Code 
# Alternative: every observation, in place of whiskers and outliers
boxplot(x_hat,
    main=NA, ylab='Murder Arrests',
    whisklty=0, staplelty=0, outline=FALSE)
title('With Raw Data', font.main=1)
# Raw Observations
stripchart(x_hat,
    pch='-', col=grey(0, .5), cex=2,
    vert=TRUE, add=TRUE)

Each depiction in this chapter suits some data types and not others. Boxplots are generally not well suited for factor data but better suited for comparing groups (discussed later in Bivariate Data). For a single data vector, the mass function needs repeating values but works for unordered factors whereas the ECDF works for discrete and continuous data requires an ordering.5 Here is a table to help you recall

Visual Unordered factor Ordered factor Discrete cardinal Continuous cardinal
Mass Function \(\hat{p}(x)\) Yes Yes Yes No
Histogram \(\hat{f}(x)\) No No Yes Yes
ECDF \(\hat{F}(x)\) No Yes Yes Yes
Boxplot No No Yes Yes

2.5 Exercises

  1. Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.

  2. Classify each of the following as cardinal (discrete or continuous) or factor (ordered or unordered): (a) temperature in Celsius, (b) letter grade on a transcript, (c) number of siblings, (d) province of residence. Briefly explain each choice.

  3. The built-in dataset mtcars has a column mpg. Compute the histogram density by hand for the bins \((10,15]\), \((15,20]\), \((20,25]\), \((25,30]\), \((30,35]\). That is, for each bin count the number of observations, then divide by \(n \times 2h\) where \(2h=5\). Verify your answer matches hist(mtcars$mpg, breaks=seq(10,35,by=5), freq=FALSE).

  4. Load USArrests in R. Compute the ECDF of the Assault variable and use it to find the proportion of states with 200 or fewer assault arrests. Then find the upper quartile of Assault using quantile.

Further Reading

Recall

This chapter sorted observations into cardinal and factor types and then showed several complementary ways to picture the distribution of a single variable: bar plots and lollipops for proportions, the histogram density, the empirical CDF, quantiles (min, median, max, quartiles), and the boxplot. Proportions suit discrete data and histograms suit continuous data, while the ECDF, quantiles, and boxplot work for either. The running three-point dataset \(\{3, 3.1, 0.02\}\) ties them together: the ECDF jumps by \(1/3\) at each of \(0.02\), \(3\), and \(3.1\) before flattening at \(1\), its rises across the unit bins \((1/3, 0, 1/3, 1/3)\) are the histogram bar areas, and the median reads off as \(3\). Each picture answers a different question about the same data, and switching between them is a habit worth building.

Frigge, Michael, David C. Hoaglin, and Boris Iglewicz. 1989. “Some Implementations of the Boxplot.” The American Statistician 43 (1): 50–54. https://doi.org/10.1080/00031305.1989.10475612.
Weissgerber, Tracey L., Natasa M. Milic, Stacey J. Winham, and Vesna D. Garovic. 2015. “Beyond Bar and Line Graphs: Time for a New Data Presentation Paradigm.” PLOS Biology 13 (4): e1002128. https://doi.org/10.1371/journal.pbio.1002128.

  1. We will not cover the statistical analysis of text in this book, but strings are amenable to statistical analysis.↩︎

  2. If \(L\) distinct bins exactly cover the range, then \(2h=[\text{max}(\hat{x}_{i}) - \text{min}(\hat{x}_{i})]/L\) and we can write the design points as \(x~\text{in}~\left\{ \text{min}(\hat{x}_{i}) + (2\ell-1)h \right\}_{\ell=1}^{L}\), since consecutive midpoints sit \(2h\) apart.↩︎

  3. Summing the identity above over the bin midpoints \(x'\) at or below \(x\) telescopes, since the upper edge of each bin is the lower edge of the next: \(\sum_{x' \leq x} 2h \times \hat{f}(x') = \sum_{x' \leq x} \left[ \hat{F}(x'+h) - \hat{F}(x'-h) \right] = \hat{F}(x+h) - \hat{F}(x_{1}-h)\) where \(x_{1}\) is the first midpoint. The lowest bin edge \(x_{1}-h\) lies below every observation, so \(\hat{F}(x_{1}-h)=0\) and the bar areas up to \(x\) add up to \(\hat{F}(x+h)\). For discrete data the same relationship holds without bins: \(\hat{F}(x)=\sum_{x' \leq x} \hat{p}(x')\) accumulates the proportions, and the size of each jump in \(\hat{F}\) is \(\hat{p}(x)\).↩︎

  4. Technically, the upper and lower hinges use two different versions of the first and third quartile. Statistical packages differ in how they compute both the quartiles and the fences that flag outliers (Frigge et al. 1989). See also https://stackoverflow.com/questions/40634693/lower-and-upper-quartiles-in-boxplot-in-r.↩︎

  5. In R, ecdf and quantile(..., type=1) accept ordered factors, but boxplot silently treats the levels as the numbers \(1, 2, 3, \ldots\).↩︎