2  Data & Visualization

---

This chapter asks three questions.

• What kind of data do we have?
  • Cardinal and factor data: differences between cardinal values always mean the same thing, while differences between factor values do not.
• How does R hold a dataset?
  • Lists and data.frames: a list can hold objects of different shapes, and a data.frame is like a matrix whose columns can have different types.
• What does the distribution of one variable look like?
  • Proportions and histograms: these show how often values occur, with the proportion \(\hat{p}(x)\) for discrete data and the histogram density \(\hat{f}(x)\) for continuous data.
  • ECDFs and quantiles: the ECDF \(\hat{F}(x)\) is the proportion of data at or below \(x\), and a quantile reverses that question.
  • Boxplots: these compress the distribution into the median, the interquartile range, and any outliers.

Which depiction to use depends on the data type.

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 determines which summaries make sense. You can average cardinal data, plot it in a scatterplot, and use it in regressions. 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 can be discrete or continuous, and factor data can be ordered or unordered. So there are 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^{\text{st}}, 2^{\text{nd}}, 3^{\text{rd}}\}\) place in a race and notice that \(1^{\text{st}}\) - \(2^{\text{nd}}\) place does not equal \(2^{\text{nd}}\) - \(3^{\text{rd}}\) place for a very competitive person who cares only about winning.
  • Unordered (categorical): E.g., \(\{\text{Amanda}, \text{Bert}, \text{Charlie}\}\) and notice that \(\text{Amanda} - \text{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 encountered in a variety of settings. You often have to format them after reading them into R, and 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

## MAIN EXAMPLE
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

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=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

TipTest Yourself

Create a data.frame storing two different types of data. Then 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

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.

The Ecdat package you installed in First Steps contains economics datasets, such as Wages1.

## MAIN EXAMPLE
# Each row is one of 3294 US workers.
# The wage is hourly, in 1980 dollars, and school and exper are in years.
data(Wages1, package='Ecdat') # every session, saved in your script
head(Wages1) # Actual Data
##   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

# 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

To list every dataset in a package, use data(package='Ecdat'). To read about one dataset, use ?Ecdat::Wages1.

R also comes with built-in datasets, which need no data() call.

# 50 US states in 1973: arrests per 100,000 residents, and percent urban.
head(USArrests)
##            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

TipTest Yourself

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

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

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

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

## MAIN EXAMPLE
# Discretized data
x_hat <- Wages1[, 'wage'] # lab data: hourly wages in 1980 dollars
x_hat[1:3]
## [1] 6.315296 5.479770 3.642170
x_rounded <- floor(x_hat) #rounded down, to be discrete
x_rounded[1:3]
## [1] 6 5 3

p_hat <- table(x_rounded)/length(x_rounded)
plot(p_hat, col=grey(0, .5),
    xlab='Hourly wage (1980 dollars, discretized)',
    ylab='Proportion of workers with each value')

# The tallest line says about 16% of workers earn from 4 up to 5 dollars.

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.

TipTest Yourself

Repeat the empirical mass function with murder arrest rates across US states.

# Discretized data
x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000
x2_hat[1:3]
## [1] 13.2 10.0  8.1
x2_rounded <- floor(x2_hat) #rounded down, to be discrete
x2_rounded[1:3]
## [1] 13 10  8

p2_hat <- table(x2_rounded)/length(x2_rounded)
plot(p2_hat, col=grey(0, .5),
    xlab='Murder Arrests (Discretized)',
    ylab='Proportion of States with each 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 shows the shape of a continuous distribution. It shows where observations are concentrated and whether there are multiple peaks. It also shows 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. Dividing each count by \(n \cdot 2h\) 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 (Scott 2015). Note that \(h\) is measured from the bin midpoint out to one edge, so a bin is \(2h\) wide.

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}~\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. ²

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

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=FALSE'
# hist(x0_hat, breaks=c(0, 4/3, 8/3, 4), plot=FALSE, right=FALSE)

## MAIN EXAMPLE
# Practical Example
hist(x_hat,
    breaks=seq(0, 40, by=1), # bin width =1, so the half-width is h=1/2
    freq=FALSE,
    border=NA,
    main=NA,
    xlab='Hourly wage (1980 dollars)')
title('Hourly Wages', font.main=1)
# Raw Observations
rug(x_hat, col=grey(0, .5))

# About three quarters of wages lie between 2 and 8 dollars.
# A long right tail reaches almost 40 dollars.

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

TipTest Yourself

Repeat the histogram with murder arrest rates.

x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000
hist(x2_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(x2_hat, 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 wages in \((2,4)\) is

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

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.

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

F_hat <- ecdf(x0_hat)

# Visualized
plot(F_hat)

# 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

## MAIN EXAMPLE
F_hat <- ecdf(x_hat)
# proportion of wages <= 10
F_hat(10)
## [1] 0.9119611
# About 91% of workers earn at most 10 dollars an hour.

# proportion of wages <= x, for all x
# The axis stops just past 20 dollars, so the 16 largest wages are not drawn.
plot(F_hat, main=NA, xlim=c(0, 20),
    xlab='Hourly wage (1980 dollars)',
    ylab='Proportion of workers with wage <= x',
    do.points=FALSE, verticals=TRUE)
title('Hourly Wage ECDF', font.main=1)
rug(x_hat[x_hat <= 20], col=grey(0, .5))

# With 3294 workers each step is tiny, so the ECDF looks like a smooth curve.

TipTest Yourself

Repeat the ECDF with murder arrest rates. What fraction of states have at most \(10\) murders per \(100\)k?

x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000
F2_hat <- ecdf(x2_hat)
# proportion of murders <= 10
F2_hat(10)
## [1] 0.7
# 70% of states have at most 10 murders per 100k.

# proportion of murders <= x, for all x
plot(F2_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(x2_hat)

The ECDF is useful when we want precise numbers rather than a picture of shape: “what fraction of the data is at most \(x\)?” 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 (Dvoretzky et al. 1956). 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.³ For a continuous example, the shaded bars on the histogram cover the observations at or below a particular value \(x\), and their total area equals the height of the colored point on the ECDF, \(\hat{F}(x)\).

F_hat <- ecdf(x_hat)
breaks <- seq(0, 40, 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.4708561
F_hat(x)
## [1] 0.4708561
# Both are about 0.47, the proportion of workers earning at most 5 dollars.

# Histogram: shade the bins at or below x
# Both panels stop just past 20 dollars, so the 16 largest wages are not drawn.
hist(x_hat, breaks=breaks, freq=FALSE, border=NA, xlim=c(0, 20),
    col=ifelse(below, col_x, grey(.8)),
    main=NA, xlab=NA, ylab='Density')
rug(x_hat[x_hat <= 20], col=grey(0, .5))

# ECDF: mark the point (x, F_hat(x))
plot(F_hat, xlim=c(0, 20), do.points=FALSE, verticals=TRUE,
    main=NA, xlab='Hourly wage (1980 dollars)', 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)

TipTest Yourself

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.

# Histogram density from the ECDF, and back again
F_hat <- ecdf(x_hat)
breaks <- seq(0, 40, 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 making precise statements about a distribution. For example, we can talk about “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}\)”. Quantiles also describe spread, such as the range or \(\text{IQR}\), and shape, such as the gap between the median and the mean. You do not need to assume a smooth model for the distribution to compute them.

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, the 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.)

x1_hat <- c(0, 0, 0.02, 3, 5)
quantile(x1_hat, probs=c(0.25, 0.5, 0.75)) # lower quartile, median, upper quartile
##  25%  50%  75% 
## 0.00 0.02 3.00

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

## MAIN EXAMPLE
# common quantiles
quantile(x_hat)
##          0%         25%         50%         75%        100% 
##  0.07655561  3.62157010  5.20578148  7.30450612 39.80891720

# All deciles are quantiles
quantile(x_hat, probs=seq(0, 1, by=.1))
##          0%         10%         20%         30%         40%         50% 
##  0.07655561  2.44228726  3.24556561  3.96669544  4.59333660  5.20578148 
##         60%         70%         80%         90%        100% 
##  5.89692780  6.76912762  7.74962940  9.74935904 39.80891720
# The top tenth of wages spans about 9.75 to 39.81 dollars, a long right tail.

# Visualized: Inverting the Empirical Distribution
# The axis stops just past 20 dollars, so the 16 largest wages are not drawn.
F_hat <- ecdf(x_hat)
plot(F_hat, lwd=2, xlim=c(0, 20),
    do.points=FALSE, verticals=TRUE, 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% 
## 3.618992 9.749460
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)

# A quarter of wages are at or below about 3.62 dollars.
# Nine tenths are at or below about 9.75 dollars.

TipTest Yourself

Repeat the quantiles with murder arrest rates.

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

# All deciles are quantiles
quantile(x2_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
F2_hat <- ecdf(x2_hat)
plot(F2_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))
q2_hat <- quantile(x2_hat, p, type=1)
q2_hat
##  25%  90% 
##  4.0 13.2
segments(q2_hat, p, -10, p, col=cols)
segments(q2_hat, p, q2_hat, 0, col=cols)
mtext( round(q2_hat, 2), 1, at=q2_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~\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.

x_sorted_hat <- sort(x_hat)
head(x_sorted_hat) # the smallest few of 3294 wages
## [1] 0.07655561 0.09799118 0.10289074 0.15348016 0.18462107 0.21231423

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

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

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

Boxplots

Boxplots also summarize the distribution. The boxplot shows the median (solid black line) and interquartile range (\(\text{IQR}=\) upper quartile \(-\) lower quartile; filled box).⁴ As a default, whiskers are shown as \(1.5\times \text{IQR}\) and values beyond that are highlighted as outliers (so whiskers do not typically show the data range). 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 \text{IQR}\) beyond the box, and that distance is only meaningful when differences are.

## MAIN EXAMPLE
# Default: median, box (IQR), whiskers at 1.5*IQR, outliers as circles
boxplot(x_hat,
    main=NA, ylab='Hourly wage (1980 dollars)')
title('Default', font.main=1)

# 96 wages sit more than 1.5*IQR above the upper quartile, drawn as circles.
# None sits that far below the lower quartile.

The default boxplot compresses \(n\) observations into five numbers and any outliers, 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
# Without ylim, outline=FALSE would leave the long right tail out of the figure.
boxplot(x_hat,
    main=NA, ylab='Hourly wage (1980 dollars)', ylim=range(x_hat),
    border=NA, medcol=grey(.4), # a darker shade of the default grey fill, for the median
    whisklty=0, staplelty=0, outline=FALSE)
title('With Raw Data', font.main=1)
# Raw Observations
stripchart(x_hat,
    pch='-', col=grey(0, .1), cex=2,
    vert=TRUE, add=TRUE)

TipTest Yourself

Repeat both boxplots with murder arrest rates.

x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000

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

# No state sits more than 1.5*IQR beyond a quartile.
# So the boxplot draws no outlier circles.

# Alternative: every observation, in place of whiskers and outliers
boxplot(x2_hat,
    main=NA, ylab='Murder Arrests',
    border=NA, medcol=grey(.4), # a darker shade of the default grey fill, for the median
    whisklty=0, staplelty=0, outline=FALSE)
title('With Raw Data', font.main=1)
# Raw Observations
stripchart(x2_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. 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 but requires an ordering.⁵ Here is a table to help you recall

Visuals by Data Type

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

• 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 stored datasets in R as lists and data.frames. It then showed several ways to picture the distribution of a single variable. These were bar plots and lollipops for proportions, the histogram density, the empirical CDF, quantiles, and the boxplot. The named quantiles were the min, median, quartiles, and max. Proportions suit discrete data and histograms suit continuous data, while the ECDF, quantiles, and boxplot work for either. The lab data were the 3294 hourly wages in Wages1. We used the three-point dataset \(\{3, 3.1, 0.02\}\) for the histogram and ECDF hand calculations. Its 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. Its median is \(3\). Each picture answers a different question about the same data, so practice switching between them.

Dvoretzky, A., J. Kiefer, and J. Wolfowitz. 1956. “Asymptotic Minimax Character of the Sample Distribution Function and of the Classical Multinomial Estimator.” The Annals of Mathematical Statistics 27 (3): 642–69. https://doi.org/10.1214/aoms/1177728174.

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.

Scott, D. W. 2015. Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley Series in Probability and Statistics. Wiley. https://books.google.de/books?id=pIAZBwAAQBAJ.

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=[\max(\hat{x}_{i}) - \min(\hat{x}_{i})]/L\) and we can write the design points as \(x~\text{in}~\left\{ \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\).