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, data.frames, arrays), and finish with the standard ways to depict a distribution: proportions, histograms, kernel densities, ECDFs, quantiles, and boxplots.

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.

TipTest Yourself

Here are some examples

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

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).¹

# 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

x <- seq(1, 10)                       # vector of 1..10
y <- 2*x                              # element-wise doubling
list(x, y)                            # 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

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.

# data.frames: your most common data type
    # matrix of different data-types
    # well-ordered lists
data.frame(x, y)  # 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

TipTest Yourself

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

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

Arrays.

Arrays are generalization of matrices to multiple dimensions. They are a very efficient way to store well-formatted numeric data, and are often used in spatial econometrics and time series (often in the form of “data cubes”).

# data square (matrix)
array(data = seq(1, 24), dim=c(3, 8))
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    1    4    7   10   13   16   19   22
## [2,]    2    5    8   11   14   17   20   23
## [3,]    3    6    9   12   15   18   21   24

# data cube
a <- array(data = seq(1, 24), dim=c(3, 2, 4))
a
## , , 1
## 
##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6
## 
## , , 2
## 
##      [,1] [,2]
## [1,]    7   10
## [2,]    8   11
## [3,]    9   12
## 
## , , 3
## 
##      [,1] [,2]
## [1,]   13   16
## [2,]   14   17
## [3,]   15   18
## 
## , , 4
## 
##      [,1] [,2]
## [1,]   19   22
## [2,]   20   23
## [3,]   21   24

a[1, , , drop = FALSE]  # Row 1
#a[, 1, , drop = FALSE]  # Column 1
#a[, , 1, drop = FALSE]  # Layer 1

a[ 1, 1,  ]  # Row 1, column 1
#a[ 1,  , 1]  # Row 1, 'layer' 1
#a[  , 1, 1]  # Column 1, 'layer' 1
a[1 , 1, 1]  # Row 1, column 1, 'layer' 1

Apply extends to arrays

apply(a, 1, mean)    # Row means
## [1] 11.5 12.5 13.5
apply(a, 2, mean)    # Column means
## [1] 11 14
apply(a, 3, mean)    # 'Layer' means
## [1]  3.5  9.5 15.5 21.5
apply(a, c(1, 2), mean)  # Row/Column combination 
##      [,1] [,2]
## [1,]   10   13
## [2,]   11   14
## [3,]   12   15

The outer function produces an array by multiplying every element of one vector with every element of another.

x <- c(1, 2, 3)
x_mat1 <- outer(x, x) # x %o% x
x_mat1
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    2    4    6
## [3,]    3    6    9
is.array(x_mat1) # Matrices are arrays
## [1] TRUE

x_mat2 <- matrix( seq(6, 1), 2, 3)
outer(x_mat2, x)
## , , 1
## 
##      [,1] [,2] [,3]
## [1,]    6    4    2
## [2,]    5    3    1
## 
## , , 2
## 
##      [,1] [,2] [,3]
## [1,]   12    8    4
## [2,]   10    6    2
## 
## , , 3
## 
##      [,1] [,2] [,3]
## [1,]   18   12    6
## [2,]   15    9    3
# outer(x_mat2, matrix(x))
# outer(x_mat2, t(x))
# outer(x_mat1, x_mat2)

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.

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
X <- c(3, 3.1, NA, 0.02) #Small dataset we will use in numerical examples
sum( is.na(X) )
## [1] 1

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 will contrast it with theoretical quantities in Population Statistics. We often analyze observations in comparison to some value \(x\).

X <- c(3, 3.1, 0.02) # Data: 'big X'
x <- 2 # particular value: 'little x'
sum(X <= x)
## [1] 1
sum(X == 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 \[\begin{eqnarray} \hat{p}_{x}=\sum_{i=1}^{n}\mathbf{1}\left(\hat{X}_{i}=x\right)/n. \end{eqnarray}\] 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\).

NoteMust Know

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\).

X <- c(1, 2, 1, 3, 1)
table(X)/length(X)
## X
##   1   2   3 
## 0.6 0.2 0.2

# Discretized data
X <- USArrests[, 'Murder']
Xr <- floor(X) #rounded down, to be discrete

#table(Xr)
proportions <- table(Xr)/length(Xr)
plot(proportions, 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 uses skinny lines with a circle on top). 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} \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. The same \(h\) reappears in the kernel density below, and again in Local Regression, where it sets how much data around a point a local model uses.

TipTest Yourself

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

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

For \(x=4\) \(\mathbf{1}\left( \hat{X}_{i} \in \left(4-h, 4+h \right] \right)=\mathbf{1}\left( 3.8 \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{eqnarray} 2h \times \hat{f}(x) &=& \frac{\sum_{i=1}^{n} \mathbf{1}\left( \hat{X}_{i} \in \left(x-h, x+h \right] \right)}{n}. \end{eqnarray}\] We compute \(\hat{f}(x)\) for each bin midpoint \(x\), and the area of all rectangles sums to one. ²

NoteMust Know

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)\).

# Intuitive Examples
X <- c(3, 3.1, 0.02)
Xhist <- hist(X, breaks=c(0, 1, 2, 3, 4), plot=FALSE, freq=FALSE)
Xhist
## $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] "X"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

# base x height
base <- 1
height <- Xhist[['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

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

# as a default, R uses bins (, ] instead of [, )
# but you can change that with 'right=F'
# hist(X, breaks=c(0, 4/3, 8/3, 4), plot=F, right=F)

# Practical Example
X <- USArrests[, 'Murder']
hist(X,
    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(USArrests[, 'Murder'], col=grey(0, .5))

# 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.

TipTest Yourself

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

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

Kernel Density.

When we want a smoother picture than the histogram’s jagged steps, we replace fixed bins with a moving window centered on each value of \(x\) we want to evaluate.

A kernel density estimates the density at \(x\) by averaging a smooth weight function (the kernel) centered on each data point, with a bandwidth \(h\) that controls how wide the weight spreads. The bandwidth is a half-width, exactly as in the histogram: each kernel below weights observations within \(h\) of \(x\) and gives zero weight beyond. Unlike the histogram, the bins around each \(x\) can overlap, which is why the result is smooth instead of stepped.

Kernel densities are useful when we want a continuous picture of shape: easier to read than a histogram and better suited to overlaying two or three distributions on the same axes, where stacked or overlapping bars become visually muddled. Different kernels weight nearby points in different ways, producing curves with visually distinct (but qualitatively similar) shapes. The simplest kernel function is the “uniform” (or “rectangular”) kernel, which places a weight of \(1/2\) on every point inside the interval \(\left[ x-h, x + h\right]\) and zero outside: \[\begin{eqnarray} k_{U}\left( \hat{X}_{i}, x, h \right) &=& \frac{\mathbf{1}\left(\frac{|\hat{X}_{i}-x|}{h} \leq 1\right)}{2} = \frac{\mathbf{1}\left( \hat{X}_{i} \in \left[ x-h, x + h\right] \right) }{2} \\ \hat{f}_{U}(x) &=& \frac{1}{nh} \sum_{i}^{n} k_{U}(\hat{X}_{i}, x, h) = \frac{ \sum_{i}^{n} \mathbf{1}\left( \hat{X}_{i} \in \left[ x-h, x + h\right] \right) }{n 2 h} \end{eqnarray}\]

For example, take the dataset \(\{3, 3.1, 0.02\}\) with bandwidth \(h=1/2\) and evaluate the density at \(x=3\). The interval \([x-h, x+h]=[2.5, 3.5]\) contains the two points \(3\) and \(3.1\), so \(\sum_{i}\mathbf{1}\left(\hat{X}_{i} \in [2.5,3.5]\right)=2\) and \(\hat{f}_{U}(3)=\frac{2}{n 2h}=\frac{2}{3\times 1}=2/3\).

Notice that the uniform kernel is essentially the histogram but without the restriction that \(x\) must be a midpoint of exclusive bins. Typically, the points \(x\) are chosen to be either the unique observations or some equidistant set of “design points” (e.g., at \(512\) evenly spaced values of \(x\) the dataset, not just the midpoints of exclusive bins).³

# Practical Example
X <- USArrests[, 'Murder']
hist(X,
    breaks=seq(0, 20, by=2), # bin width =2, so h=1
    freq=FALSE,
    border=NA,
    main=NA,
    xlab='Murder Arrests')
title('Murder Arrests', font.main=1)
# Density Estimate with half-width h=3
# 'density' reads bw as the kernel's standard deviation, so pass bw=h/sqrt(3)
h <- 3
lines( density(X, bw=h/sqrt(3), kernel='rectangular') )
# Raw Observations
rug(X, col=grey(0, .5))

We can also replace the uniform kernel with other kernel functions to create even smoother lines.

TipTest Yourself

There are many kernels, but these are the most intuitive and commonly used.

# Kernel Density Functions

x <- 0 # Design point
X <- seq(-2, 2, length.out=1001) # Where to compute weight
h <- 1 # Half-width: all three kernels are zero beyond |X-x| > h


plot.new()
plot.window(xlim=c(-1.2, 1.2), ylim=c(0, 1))

## Uniform/rectangular
d_unif <- function(X, x, h){
    u <- abs(X-x)/h
    fu <- 1/2*(u <= 1)
    return(fu)
}
kU <- d_unif(X, x, h)
kern_cols <- hcl.colors(3, alpha=.45)
lines( kU ~ X, col=kern_cols[1], lty=1)

## Epanechnikov
d_epan <- function(X, x, h){
    u <- abs(X-x)/h
    fu <- 3/4*(1-u^2)*(u <= 1)
    return(fu)
}
kE <- d_epan(X, x, h)
lines( kE ~ X, col=kern_cols[2], lty=1)
## Try others using the 'density' function
# 'density' reads bw as the kernel's standard deviation, not the half-width h
#kE <-  density(x=0, bw=h/sqrt(5), kernel='epanechnikov')

## Tricubic
d_tricub <- function(X, x, h){
    u <- abs(X-x)/h
    fu <- 70/81*(1-u^3)^3*(u <= 1)
    return(fu)
}
kT <- d_tricub(X, x, h)
lines( kT ~ X, col=kern_cols[3], lty=1)


rug(0, lwd=2)
axis(1)
axis(2)

legend('topright', lty=1, col=kern_cols,
    legend=c('Uniform (h=1)', 'Epanechnikov (h=1)', 'Tricubic (h=1)'))

Each of the three functions integrates to one over \(u=(\hat{X}_{i}-x)/h\), so the same \(\frac{1}{nh}\) rescaling above turns any of them into a density.

Once we have picked a kernel (which particular one is not particularly important) we can use it to compute density estimates at each design point.

# Practical Example
X <- USArrests[, 'Murder']
hist(X,
    breaks=seq(0, 20, by=2), # bin width 2h=2, so h=1
    freq=FALSE,
    border=NA,
    main=NA,
    xlab='Murder Arrests')
title('Murder Arrests', font.main=1)
# Same half-width h=3 as the rectangular kernel above, but a smoother weight
# For the Epanechnikov kernel, bw=h/sqrt(5)
h <- 3
lines( density(X, bw=h/sqrt(5), kernel='epanechnikov') )
# Raw Observations
rug(X, col=grey(0, .5))

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.

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.

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. 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. Inverting the picture horizontally (asking “at what value \(x\) does the proportion reach \(p\)?”) gives the quantiles introduced next.

NoteMust Know

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 \in \{0.5,2.5,3.5\}\). At \(x=0.5\), \(F(0.5)\) measures the proportion of the data \(\leq 0.5\). Since only one observations, \(0.02\), of three is \(\leq 0.5\), we can compute \(F(0.5)=1/3\). Similarly, since only one observations, \(0.02\), of three is \(\leq 2.5\), we can compute \(F(2.5)=1/3\). Since all observations are \(\leq 3.5\), we can compute \(F(3.5)=1\).

X <- c(3, 3.1, 0.02)
Fhat <- ecdf(X)

# Visualized
plot(Fhat)

# Evaluated at the data
x <- X
Fhat(X)
## [1] 0.6666667 1.0000000 0.3333333
#sum(X<=3)/length(X)
#sum(X<=3.1)/length(X)
#sum(X<=0.02)/length(X)

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

F_murder <- ecdf(USArrests[, 'Murder'])
# proportion of murders <= 10
F_murder(10)
## [1] 0.7
# proportion of murders <= x, for all x
plot(F_murder, 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(USArrests[, 'Murder'])

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 \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.

NoteMust Know

For example, if \(X=(0,0,0.02,3,5)\) then the median is \(0.02\), the lower quartile is \(0\), and the upper quartile is \(3\). (The number \(0\) is also special: the most frequent observation is called the mode, which we revisit in Numerical Statistics.)

X <-  c(3.1, 3, 0.02)
quantile(X, 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.

# common quantiles
X <- USArrests[, 'Murder']
quantile(X)
##     0%    25%    50%    75%   100% 
##  0.800  4.075  7.250 11.250 17.400

# All deciles are quantiles
quantile(X, 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
FX_hat <- ecdf(X)
plot(FX_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))
QX_hat <- quantile(X, p, type=1)
QX_hat
##  25%  90% 
##  4.0 13.2
segments(QX_hat, p, -10, p, col=cols)
segments(QX_hat, p, QX_hat, 0, col=cols)
mtext( round(QX_hat, 2), 1, at=QX_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\).

TipTest Yourself

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 \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.

X <- USArrests[, 'Murder']
Xo <- sort(X)
Xo
##  [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
Xo[ceiling(length(Xo)*.5)]
## [1] 7.2
quantile(X, probs=.5, type=1) #tie break rule #1
## 50% 
## 7.2

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

# min is a special case
min(Xo)
## [1] 0.8
Xo[ceiling(length(Xo)*0)]
## numeric(0)
quantile(X, probs=0, type=1)
##  0% 
## 0.8
Xo[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).⁴ 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. They show up in Numerical Statistics, where heavier-tailed data push observations past the whiskers.

# Default: median, box (IQR), whiskers at 1.5*IQR, outliers as circles
boxplot(USArrests[, 'Murder'],
    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).

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

These univariate visualization tools extend to examining relationships between two variables. See Bivariate Distributions for how histograms, ECDFs, and boxplots adapt to bivariate data.

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.

• https://library.virginia.edu/data/articles/understanding-empirical-cumulative-distribution-functions – Clear walkthrough of the ECDF with visual examples.
• https://raw.githubusercontent.com/rstudio/cheatsheets/main/regex.pdf – RStudio cheat sheet for regular expressions and string handling in R.
• https://www.econometrics-with-r.org/2.1-random-variables-and-probability-distributions.html – Review of data types and distributions in an econometrics context.

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 smoother kernel density, the empirical CDF, quantiles (min, median, max, quartiles), and the boxplot. The running three-point dataset \(\{3, 3.1, 0.02\}\) ties three of them together: the uniform kernel density at \(x=3\) with half-width \(h=1/2\) equals \(2/3\) (two points fall inside the interval \([2.5, 3.5]\)); the ECDF jumps by \(1/3\) at each of \(0.02\), \(3\), and \(3.1\) before flattening at \(1\); 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. In the next chapter we summarize the same distributions with single numbers: mean, variance, median, MAD, skewness, and kurtosis.

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 \in \left\{ \text{min}(\hat{X}_{i}) + (2\ell-1)h \right\}_{\ell=1}^{L}\), since consecutive midpoints sit \(2h\) apart.

3. One detail matters when moving from the formula to R. The density function scales its bw argument to be the standard deviation of the kernel rather than the half-width, so bw and \(h\) are different numbers. Setting bw from \(h\), rather than reading \(h\) off bw, keeps the code matching the math. The conversion is a constant that depends on the kernel: the half-width is \(\sqrt{3}\) times bw for the rectangular kernel and \(\sqrt{5}\) times bw for the Epanechnikov kernel.

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.