Skip to main content

3  Summary Statistics

We often summarize distributions with statistics: functions of data. We will go through some specific examples mathematically below, but can intuitively understand that they are generally computed as statistic <- function(x){ .... }.

Note that functions can take functions as arguments, meaning we can also program statistics generally as

Code 
statistic <- function(x_hat, f){
    stat_hat <- f(x_hat)
    return(stat_hat)
}

x0_hat <- c(0, 1, 3, 10, 6) # Data
statistic(x0_hat, sum)
## [1] 20

This chapter asks four questions.

The most basic way to compute statistics is with summary, which reports multiple values that can all be calculated individually.

Code 
## MAIN EXAMPLE
# A sample of real data: 3294 US workers
data(Wages1, package='Ecdat')
x_hat <- Wages1[, 'wage'] # lab data: hourly wages in 1980 dollars
summary(x_hat)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.07656  3.62157  5.20578  5.75758  7.30451 39.80892

Together, the sample mean and variance statistics summarize the central tendency and dispersion of a distribution. In some special cases, such as with the normal distribution, they completely describe the distribution. Other distributions are better described with other statistics, either as an alternative or in addition to the mean and variance.

3.1 Mean and Variance

The mean and variance are the two most basic statistics that summarize the center and how spread apart the values are for data in your sample. As before, we represent data as vector \(\hat{x}=(\hat{x}_{1}, \hat{x}_{2}, ....\hat{x}_{n})\), where there are \(n\) observations and \(\hat{x}_{i}\) is the value of the \(i\)th one.

Mean

The first step in describing a sample is to reduce its center to a single number.

ImportantKey Definition

The sample mean is the sum of the observations divided by the number of observations, \[\hat{m} = \frac{1}{n} \sum_{i=1}^{n} \hat{x}_{i}.\] When data are discrete, we can also compute the mean using proportions, \[\hat{m}=\sum_{x} x \hat{p}(x).\]

The mean is also called the sample mean or empirical mean. It is the simplest and most commonly reported single-number summary of a sample. It is the natural measure of central tendency when the data are roughly symmetric. The mean is useful when totals matter: total income across a population equals the mean income times the number of people, a relation that does not hold for the median or other rank-based summaries. Because every observation contributes equally to the sum, a single extreme value can pull the mean substantially in its direction.

For example, a dataset of \(\{4,1,10\}\) has a mean of \([4+1+10]/3=5\).

Code 
x0_hat <- c(4, 1, 10)
sum(x0_hat)/length(x0_hat)
## [1] 5
mean(x0_hat)
## [1] 5

E.g., the dataset \(\{1,2,1,3,1\}\) has \(\hat{p}(1)=\frac{3}{5}\), \(\hat{p}(2)=\frac{1}{5}\), \(\hat{p}(3)=\frac{1}{5}\). Then, sorting the dataset as \(\{1,1,1,2,3\}\), we can see that \(\hat{m} = [1+1+1+2+3]/5 = [1+1+1]/5 + 2/5 + 3/5 = 1 \frac{3}{5} + 2 \frac{1}{5} + 3 \frac{1}{5} = 1 \hat{p}(1) + 2 \hat{p}(2) + 3 \hat{p}(3)\). In either case, we end up with \(\frac{3+2+3}{5}=8/5=1.6\)

Code 
x0_hat <- c(1, 2, 1, 3, 1)
mean(x0_hat)
## [1] 1.6

proportions <- table(x0_hat)/length(x0_hat)
proportions
## x0_hat
##   1   2   3 
## 0.6 0.2 0.2
vals <- sort(unique(x0_hat))
sum(vals*proportions)
## [1] 1.6

Variance

Once we know where the data are centered, the next question is how widely they vary around that center.

ImportantKey Definition

The sample variance is the average squared deviation from the mean, \[\hat{v} = \frac{1}{n} \sum_{i=1}^{n} (\hat{x}_{i} - \hat{m})^2.\] When data are discrete, we can also compute the variance using proportions, \[\hat{v}=\sum_{x} (x-\hat{m})^2 \hat{p}(x).\] The sample standard deviation is its square root, \(\hat{s} = \sqrt{\hat{v}}\), and is in the same units as the data.

The variance and standard deviation are the second-most-common single-number summaries of a sample. They tell us how spread out the data are around the mean. A small variance means the observations are close to the mean, and a large variance means many are far from it. The standard deviation in particular is useful because it sits in the original units of the data, making sentences like “monthly returns averaged \(0.6\%\) with a standard deviation of \(4\%\)” interpretable at a glance. Squaring the deviations also means a few extreme values can dominate the result, which is one reason the median absolute deviation (below) is often preferred for heavy-tailed data.

For example, a dataset of \(\{4,1,10\}\) has a mean of \([4+1+10]/3=5\). The variance is \([(4-5)^2+(1-5)^2+(10-5)^2]/3=[1+16+25]/3=42/3=14\) and the standard deviation is \(\sqrt{14}\).

Code 
x0_hat <- c(4, 1, 10)
m_hat <- mean(x0_hat)
v_hat <- mean( (x0_hat - m_hat)^2 )
sqrt(v_hat)
## [1] 3.741657

E.g., the dataset \(\{1,2,1,3,1\}\) from above has \(\hat{m}=\frac{8}{5}\) and \(\hat{p}(1)=\frac{3}{5}\), \(\hat{p}(2)=\frac{1}{5}\), \(\hat{p}(3)=\frac{1}{5}\). The squared deviations of the distinct values are \((1-\frac{8}{5})^2=\frac{9}{25}\), \((2-\frac{8}{5})^2=\frac{4}{25}\), and \((3-\frac{8}{5})^2=\frac{49}{25}\). Then, sorting the dataset as \(\{1,1,1,2,3\}\), we can see that \(\hat{v} = [\frac{9}{25}+\frac{9}{25}+\frac{9}{25}+\frac{4}{25}+\frac{49}{25}]/5 = \frac{9}{25} \times \frac{3}{5} + \frac{4}{25} \times \frac{1}{5} + \frac{49}{25} \times \frac{1}{5}\). In either case, we end up with \(\frac{27+4+49}{125}=\frac{80}{125}=\frac{16}{25}=0.64\), so the standard deviation is \(\hat{s}=\sqrt{16/25}=4/5=0.8\).

Code 
x0_hat <- c(1, 2, 1, 3, 1)
m_hat <- mean(x0_hat)
mean( (x0_hat - m_hat)^2 )
## [1] 0.64

proportions <- table(x0_hat)/length(x0_hat)
vals <- sort(unique(x0_hat))
sum( (vals - m_hat)^2 * proportions )
## [1] 0.64

When data are binary, taking only the values \(0\) and \(1\), both formulas simplify further. Because \(\hat{p}(0)=1-\hat{p}(1)\), the mean is \(\hat{m}=0 \cdot \hat{p}(0) + 1 \cdot \hat{p}(1)=\hat{p}(1)\), the proportion of ones. Show that the variance is \(\hat{v}=\hat{p}(1)[1-\hat{p}(1)]\).

Substituting \(\hat{m}=\hat{p}(1)\) into the variance formula, \[\begin{aligned} \hat{v} &= [0-\hat{p}(1)]^2 \hat{p}(0) + [1-\hat{p}(1)]^2 \hat{p}(1) \\ &= \hat{p}(1)^2 [1-\hat{p}(1)] + [1-\hat{p}(1)]^2 \hat{p}(1) \\ &= \hat{p}(1)[1-\hat{p}(1)] \left( \hat{p}(1) + 1-\hat{p}(1) \right) \\ &= \hat{p}(1)[1-\hat{p}(1)]. \end{aligned}\] The variance is largest, \(1/4\), when half the observations are ones, and is zero when every observation is the same.

Check this on the lab data.

Code 
# 1 if a worker earns more than 10 dollars an hour, 0 otherwise
x1_hat <- as.numeric(x_hat > 10)
sum(x1_hat) # 290 of the 3294 workers
## [1] 290
p1_hat <- mean(x1_hat) # proportion of ones, which is also the mean
p1_hat*(1 - p1_hat)
## [1] 0.08028802

# compare to the definition
mean( (x1_hat - p1_hat)^2 )
## [1] 0.08028802
Code 
## MAIN EXAMPLE
# compute the mean and standard deviation of the lab data
m_hat <- mean(x_hat)
s1_hat <- sqrt(mean((x_hat - m_hat)^2)) # divides by n, as defined above; sd() divides by n-1
m1_sd_bounds <- c(m_hat - s1_hat,  m_hat + s1_hat)
m2_sd_bounds <- c(m_hat - 2*s1_hat,  m_hat + 2*s1_hat)
m2_sd_bounds
## [1] -0.779794 12.294964
range(x_hat) # smallest and largest observation
## [1]  0.07655561 39.80891720
# The upper bound lands inside the data.
# The lower bound, about -0.78, is below the smallest observation.
# It is also below zero, which no wage can be.

# visualize on a histogram, widened to include both sets of bounds
# The 18 wages above 20 dollars are not drawn
hist(x_hat[x_hat <= 20], breaks=seq(0, 20, by=1), border=NA, freq=FALSE, main=NA,
    xlim=c(m2_sd_bounds[1], 20),
    xlab='Hourly wage (1980 dollars)')
rug(x_hat[x_hat <= 20], col=grey(0, .1))
abline(v=m_hat, col=rgb(1, 0, 0, .8), lwd=2) # mean
abline(v=m1_sd_bounds, col=rgb(0, 0, 0, .8), lwd=2) # mean +- sd
abline(v=m2_sd_bounds, col=rgb(0, 0, 0, .8), lwd=2, lty=2) # mean +- 2 sd
mtext(expression(hat(m)-2*hat(s), hat(m)-hat(s), hat(m), hat(m)+hat(s), hat(m)+2*hat(s)),
    side=3, at=c(m2_sd_bounds[1], m1_sd_bounds[1], m_hat, m1_sd_bounds[2], m2_sd_bounds[2]),
    col=c(rgb(0, 0, 0, .8), rgb(0, 0, 0, .8), rgb(1, 0, 0, .8), rgb(0, 0, 0, .8), rgb(0, 0, 0, .8)))

Code 
# The bounds are symmetric, but wages are not.
# They stop at zero on the left and stretch much further to the right.

The bounds are symmetric around the mean by construction, but the data need not be. So a bound can land outside the data, or even on a value the variable cannot take.

Repeat this figure with the murder arrest rates in USArrests.

Code 
x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000, one value per state

# compute the mean and standard deviation
m2_hat <- mean(x2_hat)
s2_hat <- sqrt(mean((x2_hat - m2_hat)^2)) # divides by n, as defined above; sd() divides by n-1
m1_sd_bounds2 <- c(m2_hat - s2_hat,  m2_hat + s2_hat)
m2_sd_bounds2 <- c(m2_hat - 2*s2_hat,  m2_hat + 2*s2_hat)
m2_sd_bounds2
## [1] -0.8354694 16.4114694
range(x2_hat) # smallest and largest observation
## [1]  0.8 17.4
# The upper bound lands inside the data.
# The lower bound, about -0.8, is below the smallest observation.
# It is also below zero, which no arrest rate can be.

# visualize on a histogram, widened to include both sets of bounds
hist(x2_hat, breaks=15, border=NA, freq=FALSE, main=NA,
    xlim=range(c(x2_hat, m2_sd_bounds2)),
    xlab='Murder arrests (per 100k)')
rug(x2_hat, col=grey(0, .5))
abline(v=m2_hat, col=rgb(1, 0, 0, .8), lwd=2) # mean
abline(v=m1_sd_bounds2, col=rgb(0, 0, 0, .8), lwd=2) # mean +- sd
abline(v=m2_sd_bounds2, col=rgb(0, 0, 0, .8), lwd=2, lty=2) # mean +- 2 sd
mtext(expression(hat(m)-2*hat(s), hat(m)-hat(s), hat(m), hat(m)+hat(s), hat(m)+2*hat(s)),
    side=3, at=c(m2_sd_bounds2[1], m1_sd_bounds2[1], m2_hat, m1_sd_bounds2[2], m2_sd_bounds2[2]),
    col=c(rgb(0, 0, 0, .8), rgb(0, 0, 0, .8), rgb(1, 0, 0, .8), rgb(0, 0, 0, .8), rgb(0, 0, 0, .8)))

Code 
# The bounds are symmetric, but murder arrests are not.
# They stop at zero on the left and stretch further to the right.

Note that an adjusted version of the empirical variance, also called the “unbiased version”, is used by R and many statisticians: \(\hat{v}_{\text{adj}} =\frac{\sum_{i=1}^{n} [\hat{x}_{i} - \hat{m}]^2}{n-1}\) and \(\hat{s}_{\text{adj}} = \sqrt{\hat{v}_{\text{adj}}}\). In this book, we use the version defined previously because we do not yet know about “bias/unbiased” statistics. Do not be concerned, as there is hardly any difference when \(n\) is large: e.g., \(\frac{1}{n}\approx \frac{1}{n-1}\) for \(n=100,000\).

Code 
x0_hat <- c(4, 1, 10)

m_hat <- mean(x0_hat)
v_hat <- mean( (x0_hat - m_hat)^2 )
v_hat
## [1] 14

# Adjusted Version
n <- length(x0_hat)
v_adj_hat <- sum( (x0_hat - m_hat)^2 )/(n-1)
v_adj_hat
## [1] 21

var(x0_hat) # R-Version
## [1] 21

We can see how changing the mean or the variance alters a distribution. The figure below compares two normal distributions with the same standard deviation, \(1\), but different means, \(0\) and \(2\). Increasing the mean slides the whole distribution to the right without changing its shape. On the histogram, the bars keep the same heights and move over by \(2\). On the ECDF, the curve shifts right by \(2\) at every height.

The next figure compares two normal distributions with the same mean, \(0\), but different standard deviations, \(1\) and \(2\) (variances \(1\) and \(4\)). Increasing the variance spreads the distribution out on both sides of the mean. On the histogram, the distribution becomes wider and flatter, since the same total area covers a wider range. On the ECDF, the curve rises more slowly, but both curves still cross \(0.5\) at the mean.

Note that the mean and standard deviation affect a histogram or an ECDF, but cannot be read directly from them. We have to compute them separately. This is unlike the median and \(\hat{\text{IQR}}\) below, which can be read directly from the ECDF.

3.2 Other Center/Spread Statistics

A general rule of applied statistics is that there are multiple ways to measure something. Mean and Variance are measurements of Center and Spread, but there are others that have different theoretical properties and may be better suited for your dataset.

Center/Spread based on Quantiles

When the data are skewed or contain outliers, statistics built from quantiles often give a more representative summary of center and spread than the mean and variance. Every quantile lies between the smallest and largest observation, so a quantile-based range can never produce an impossible value, as \(\hat{m}-2\hat{s}\) can.

ImportantKey Definition

The sample median \(\tilde{m}\) is the value where one half of the data fall below and one half above; it is the \(0.5\) quantile.

The interquartile range (\(\hat{\text{IQR}}\)) is the upper quartile minus the lower quartile, the spread of the middle half of the data.

The median absolute deviation (\(\hat{\text{MAD}}\)) is the median of the absolute deviations from the median, \(\text{Med}(|\hat{x}_{i} - \tilde{m}|)\) (Hampel 1974; Rousseeuw and Croux 1993).

These three statistics are built from ranks, so extreme values have little influence on them. A single outlier can pull the mean by an arbitrary amount, but it moves the median by at most one position. The \(\hat{\text{IQR}}\) and \(\hat{\text{MAD}}\) are robust for the same reason, since they are built from medians and other \(p\)-th quantiles. Use these robust summaries when the distribution is skewed or has values you cannot trust to be typical. For example, income and house prices have a few very large values that can pull the mean far from a typical observation.

Code 
## MAIN EXAMPLE
median(x_hat)
## [1] 5.205781
quantile(x_hat, prob=0.5)
##      50% 
## 5.205781
# The median, about 5.21, is below the mean, about 5.76.
# The long right tail pulls the mean up, but not the median.

Examine robustness to an extreme value

Code 
x1_extreme_hat <- c(x_hat, 1000) # add one extreme value
#par(mfrow=c(1, 2)) # visualize side-by-side
#hist(x_hat)
#hist(x1_extreme_hat)

# Which measures of central tendency are robust
# to a single extreme value?
mean( x_hat)
## [1] 5.757585
mean( x1_extreme_hat )
## [1] 6.059328

quantile(x_hat, prob=0.5)
##      50% 
## 5.205781
quantile(x1_extreme_hat, prob=0.5)
##      50% 
## 5.205781

The \(\hat{\text{IQR}}\) corresponds to the box in the boxplot (Tukey 1977), and the \(\hat{\text{MAD}}\) uses the same robust idea but measures distances from the median directly rather than reading off quartile cut-offs.

Code 
## MAIN EXAMPLE
IQR(x_hat)
## [1] 3.682936
mad(x_hat, constant=1) # constant=1 gives the MAD as defined above
## [1] 1.837335

Unlike the mean and standard deviation, the median and \(\hat{\text{IQR}}\) can be read directly from the ECDF. The median is the value where \(\hat{F}\) crosses \(0.5\). The \(\hat{\text{IQR}}\) is the horizontal distance between where \(\hat{F}\) crosses \(0.25\) and where it crosses \(0.75\).

Code 
# median and quartiles
p <- c(.25, .5, .75)
q_hat <- quantile(x_hat, probs=p)

# visualize on the ECDF
# The axis stops just past 20 dollars, so the 16 largest wages are not drawn
F_hat <- ecdf(x_hat)
plot(F_hat, verticals=TRUE, do.points=FALSE, col.01line=NA, lwd=2, xlim=c(0, 20), main=NA,
    xlab='Hourly wage (1980 dollars)', ylab='Proportion <= x')
cols <- c(rgb(0, 0, 1, .8), rgb(1, 0, 0, .8), rgb(0, 0, 1, .8))
segments(q_hat, p, -10, p, col=cols, lwd=c(1, 2, 1)) # read across from p
segments(q_hat, p, q_hat, -1, col=cols, lwd=c(1, 2, 1)) # read down to the quantile
mtext(round(q_hat, 2), side=1, at=q_hat, col=cols)

# IQR is the distance between the quartiles
arrows(q_hat[1], .9, q_hat[3], .9, code=3, length=.1, col=rgb(0, 0, 1, .8))
text(mean(q_hat[c(1, 3)]), .9, 'IQR', pos=3, col=rgb(0, 0, 1, .8))

Repeat these statistics with the murder arrest rates in USArrests, and read the quartiles off the ECDF.

Code 
x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000, one value per state
median(x2_hat)
## [1] 7.25
IQR(x2_hat)
## [1] 7.175
mad(x2_hat, constant=1)
## [1] 3.65

# median and quartiles
q2_hat <- quantile(x2_hat, probs=p)

# visualize on the ECDF
F2_hat <- ecdf(x2_hat)
plot(F2_hat, pch=16, col=grey(0, .5), main=NA,
    xlab='Murder arrests (per 100k)', ylab='Proportion <= x')
segments(q2_hat, p, -10, p, col=cols, lwd=c(1, 2, 1)) # read across from p
segments(q2_hat, p, q2_hat, -1, col=cols, lwd=c(1, 2, 1)) # read down to the quantile
mtext(round(q2_hat, 2), side=1, at=q2_hat, col=cols)

# IQR is the distance between the quartiles
arrows(q2_hat[1], .9, q2_hat[3], .9, code=3, length=.1, col=rgb(0, 0, 1, .8))
text(mean(q2_hat[c(1, 3)]), .9, 'IQR', pos=3, col=rgb(0, 0, 1, .8))

Compute the \(\hat{\text{IQR}}\) statistic for the dataset \(\{-100,1,4,10,10\}\).

\(\hat{\text{IQR}} =\) Upper quartile \(-\) Lower quartile \(= 10 - 1 = 9\).

Code 
x0_hat <- c(-100, 1, 4, 10, 10)
# An intuitive alternative to sd(x0_hat), used in the boxplot
quants <- quantile(x0_hat, probs=c(.25, .75))
quants[2]-quants[1]
## 75% 
##   9
IQR(x0_hat)
## [1] 9

Compute the \(\hat{\text{MAD}}\) statistic for the same dataset \(\{-100,1,4,10,10\}\). First compute the median, \(\text{Med}(-100,1,4,10,10)=4\). Then compute \(\text{Med}(|-100-4|,~ |1-4|,~ |4-4|,~ |10-4|,~ |10-4|) = \text{Med}(104, 3, 0, 6, 6) = 6\).

Code 
#Another alternative to sd(x0_hat)
mad(x0_hat, constant=1)
## [1] 6

# Computationally equivalent
# med_hat <- quantile(x0_hat, probs=.5)
# absdev_hat <- abs(x0_hat -med_hat)
# quantile(absdev_hat, probs=.5)

Compare the robustness of various “spread” metrics to an extreme value

Code 
# Standard deviation, which divides by n, unlike sd()
sqrt(mean((x_hat - mean(x_hat))^2))
## [1] 3.26869
sqrt(mean((x1_extreme_hat - mean(x1_extreme_hat))^2))
## [1] 17.62372

IQR(x_hat)
## [1] 3.682936
IQR(x1_extreme_hat)
## [1] 3.682427

mad(x_hat, constant=1)
## [1] 1.837335
mad(x1_extreme_hat, constant=1)
## [1] 1.837335

Note that there are other “absolute deviation” statistics

Code 
# sometimes seen elsewhere
mean( abs(x_hat - mean(x_hat)) )
mean( abs(x_hat - median(x_hat)) )

Weighted Statistics

The mean generalizes to a weighted mean: an average where different values contribute to the final result with varying levels of importance. For each outcome \(x\) we have a weight \(W(x)\) and compute \[\begin{aligned} \hat{m} &= \frac{\sum_{x} x W(x)}{\sum_{x} W(x)} = \sum_{x} x w(x), \end{aligned}\] where \(w(x)=\frac{W(x)}{\sum_{x'} W(x')}\) is the normalized version of \(W(x)\), so that the weights of all unique values \(x'\) sum to unity, \(\sum_{x'}w(x')=1\). You have already seen the special case of “probability weights” when computing the mean above, with \(w(x)=\hat{p}(x)=\sum_{i=1}^{n}\mathbf{1}\left(\hat{x}_{i}=x\right)/n\).

Weights can also be given to individual observations rather than to unique values. For each observation \(i\) we have a weight \(W_{i}\) and compute \[\begin{aligned} \hat{m} &= \frac{\sum_{i=1}^{n} \hat{x}_{i} W_{i}}{\sum_{i=1}^{n} W_{i}} = \sum_{i=1}^{n} \hat{x}_{i} w_{i}, \end{aligned}\] where \(w_{i}=\frac{W_{i}}{\sum_{i'=1}^{n} W_{i'}}\) is the normalized weight, so that \(\sum_{i=1}^{n} w_{i}=1\). Equal weights, \(w_{i}=1/n\), give the ordinary mean. Adding up the weights of all observations with the same value gives the weight for that value, \(w(x)=\sum_{i=1}^{n} w_{i}\mathbf{1}\left(\hat{x}_{i}=x\right)\).

For another example, suppose a student has these scores

Code 
Homework1 <- c(score=88, weight=0.25)
Homework2 <- c(score=92, weight=0.25)
Exam1 <- c(score=67, weight=0.2)
Exam2 <- c(score=90, weight=0.3)

Grades <- rbind(Homework1, Homework2, Exam1, Exam2)
Grades
##           score weight
## Homework1    88   0.25
## Homework2    92   0.25
## Exam1        67   0.20
## Exam2        90   0.30

Each assessment is an observation \(i\) with its own weight \(w_{i}\), and the weights already sum to \(1\). We can compute the final grade as a weighted mean

Code 
# Manual Way
88*0.25 + 92*0.25 + 67*0.2 + 90*0.3
## [1] 85.4

# Computerized Way
Values <- Grades[, 'score'] * Grades[, 'weight']
FinalGrade <- sum(Values)
FinalGrade
## [1] 85.4

# Built-in function, which also normalizes the weights
weighted.mean(Grades[, 'score'], w=Grades[, 'weight'])
## [1] 85.4

In principle, we can also compute other weighted statistics such as a weighted median.

See that we can also compute weighted quantiles

Code 
weighted.quantile <- function(x_hat, w, probs){
    #See spatstat.univar::weighted.quantile
    oo <- order(x_hat)
    x_hat <- x_hat[oo]
    w <- w[oo]
    F_hat <- cumsum(w)/sum(w)
    quantile_id <- max(which(F_hat <= probs))+1
    q_hat <- x_hat[quantile_id]
    return(q_hat)
}

## Unweighted
quantile(x0_hat, probs=.5)
## 50% 
##   4
weights <- rep(1, length(x0_hat))
weighted.quantile(x_hat=x0_hat, w=weights, probs=.5)
## [1] 4

## Weighted
weights <- seq(x0_hat)
weights <- weights/sum(weights) # normalize
weighted.quantile(x_hat=x0_hat, w=weights, probs=.5)
## [1] 10

Standardized Values

Whether an observation is unusual depends on the center and spread of the data, not on its raw value alone. The standardized value of an observation is its deviation from the mean measured in standard deviations, \[\hat{z}_{i} = \frac{\hat{x}_{i} - \hat{m}}{\hat{s}}.\]

Standardized values are also called z-scores. A standardized value of \(2\) means the observation lies two standard deviations above the mean, and a value of \(-1\) means it lies one standard deviation below. The units of \(\hat{x}_{i}-\hat{m}\) and \(\hat{s}\) cancel, so standardized values have no units. Standardized values are useful for flagging observations that are far from the center of their own distribution, and for comparing observations across variables measured on different scales, such as a worker’s wage and years of schooling.

Compute the standardized values of the dataset \(\{1,2,1,3,1\}\) from above.

The dataset has \(\hat{m}=1.6\) and \(\hat{s}=0.8\). The standardized values are \((1-1.6)/0.8=-0.75\) for each \(1\), \((2-1.6)/0.8=0.5\), and \((3-1.6)/0.8=1.75\). The observation \(3\) is the most unusual, lying \(1.75\) standard deviations above the mean.

Code 
x0_hat <- c(1, 2, 1, 3, 1)
m_hat <- mean(x0_hat)
s_hat <- sqrt(mean((x0_hat - m_hat)^2))
z0_hat <- (x0_hat - m_hat)/s_hat
z0_hat
## [1] -0.75  0.50 -0.75  1.75 -0.75
Code 
# standardize the lab data, using the standard deviation that divides by n
m_hat <- mean(x_hat)
s_hat <- sqrt(mean((x_hat - m_hat)^2))
z_hat <- (x_hat - m_hat)/s_hat
# scale(x_hat) is similar, but divides by sd(), the n-1 version

# mean 0 and standard deviation 1, up to rounding error
mean(z_hat)
## [1] 1.24078e-16
sqrt(mean((z_hat - mean(z_hat))^2))
## [1] 1

# same shape, new units: standardize the bin edges too
# The 18 wages above 20 dollars are not drawn in either panel
drawn <- x_hat <= 20
f_hat <- hist(x_hat[drawn], breaks=seq(0, 20, by=1), plot=FALSE)
par(mfrow=c(1, 2), mar=c(4, 4, 1, 1))
hist(x_hat[drawn], breaks=f_hat[['breaks']], border=NA, freq=FALSE, main=NA,
    xlab='Hourly wage (1980 dollars)')
hist(z_hat[drawn], breaks=(f_hat[['breaks']] - m_hat)/s_hat, border=NA, freq=FALSE, main=NA,
    xlab='Standardized hourly wage')

Standardizing moves the center to \(0\) and changes the units on the horizontal axis, but the two histograms have the same shape.

Chebyshev’s inequality (Chebyshev 1867) says that, for any dataset, the proportion of observations at least \(k\) standard deviations from the mean is at most \(1/k^2\), \[\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}\left(|\hat{z}_{i}| \geq k\right) \leq \frac{1}{k^2}.\] For example, at most a quarter of the observations can be two or more standard deviations from the mean, whatever the shape of the distribution. To see why, note that the standardized values satisfy \(\frac{1}{n}\sum_{i=1}^{n}\hat{z}_{i}^2=1\), and every observation with \(|\hat{z}_{i}| \geq k\) adds at least \(k^2/n\) to that sum.

Check the bound on the lab data for \(k=2\) and \(k=3\).

Code 
# proportion of workers at least k standard deviations from the mean
mean(abs(z_hat) >= 2)
## [1] 0.0343048
mean(abs(z_hat) >= 3)
## [1] 0.01305404
# About 3.4% of workers are at least 2 standard deviations from the mean.
# About 1.3% are at least 3.

# Chebyshev's upper bound
1/c(2, 3)^2
## [1] 0.2500000 0.1111111
# Both proportions are far below what the bound allows.

The bound is loose because it must hold for every dataset, including ones with much heavier tails.

Standardizing also lets us compare variables measured on different scales. Their raw numbers cannot be compared directly, but their standardized values can. Below, we standardize two variables and compare them for two observations.

Code 
## MAIN EXAMPLE
# a function that standardizes any variable
standardize <- function(x_hat){
    m_hat <- mean(x_hat)
    s_hat <- sqrt(mean((x_hat - m_hat)^2))
    z_hat <- (x_hat - m_hat)/s_hat
    return(z_hat)
}

# standardize two variables measured on different scales
# Wages are in dollars and schooling is in years.
z_wages_hat <- data.frame(
    wage=standardize(Wages1[, 'wage']),
    school=standardize(Wages1[, 'school']))

# raw and standardized values for two workers
# which.max() returns the row number of the first largest value
workers <- c(which.max(Wages1[, 'wage']), which.max(Wages1[, 'school']))
Wages1[workers, c('wage', 'school')]
##           wage school
## 1813 39.808917     11
## 50    5.447893     16
round(z_wages_hat[workers, ], 2)
##       wage school
## 1813 10.42  -0.38
## 50   -0.09   2.64
# Worker 1813 has the highest wage, 10.42 standard deviations above the mean.
# That worker's schooling is 0.38 below the mean.
# Worker 50 has the most schooling, 16 years, which is 2.64 above the mean.
# That worker's wage is 0.09 below the mean.

Repeat the comparison with the murder and assault arrest rates in USArrests, for the state with the highest rate of each.

Code 
# standardize two variables measured on different scales
# Both rates are per 100,000 people, but assaults are far more common.
# So the raw numbers cannot be compared directly.
z2_hat <- data.frame(
    Murder=standardize(USArrests[, 'Murder']),
    Assault=standardize(USArrests[, 'Assault']),
    row.names=rownames(USArrests))

# raw and standardized values for two states
# Georgia has the highest murder arrest rate of any state.
# North Carolina has the highest assault arrest rate.
states <- c('Georgia', 'North Carolina')
USArrests[states, c('Murder', 'Assault')]
##                Murder Assault
## Georgia          17.4     211
## North Carolina   13.0     337
round(z2_hat[states, ], 2)
##                Murder Assault
## Georgia          2.23    0.49
## North Carolina   1.21    2.02
# Georgia's murder rate is 2.23 standard deviations above the mean.
# Its assault rate is only 0.49 above.
# North Carolina is the reverse, with an assault rate 2.02 above the mean.
# Its murder rate is 1.21 above.

Replacing the mean and standard deviation in the standardized value with the median and \(\hat{\text{MAD}}\) gives a robust z-score, \[\tilde{z}_{i} = \frac{\hat{x}_{i} - \tilde{m}}{\hat{\text{MAD}}},\] which a single extreme value cannot distort by inflating the spread.

Take the worker with the highest hourly wage in the lab data. Compute that worker’s standard and robust z-scores, with and without the extreme value in x1_extreme_hat. Which one still flags this worker as unusual?

Code 
# robust version of the standardize function above
standardize_robust <- function(x_hat){
    med_hat <- median(x_hat)
    mad_hat <- mad(x_hat, constant=1)
    z_hat <- (x_hat - med_hat)/mad_hat
    return(z_hat)
}

# the worker with the highest hourly wage, about 39.81 dollars
id <- which.max(x_hat)

# standard z-scores, without and with the extreme value
standardize(x_hat)[id]
## [1] 10.41743
standardize(x1_extreme_hat)[id]
## [1] 1.915009
# The extreme value inflates the standard deviation from about 3.3 to 17.6.
# So this worker's z-score falls from about 10.42 to about 1.92.
# The highest wage in the lab data no longer looks very unusual.

# robust z-scores, without and with the extreme value
standardize_robust(x_hat)[id]
## [1] 18.83333
standardize_robust(x1_extreme_hat)[id]
## [1] 18.83333
# The robust z-score stays at about 18.83.
# The MAD, about 1.84, is smaller than the standard deviation.
# So the robust z-scores are larger than the standard ones.

A new extreme value can inflate the standard deviation enough to hide an observation that was unusual before. The median and \(\hat{\text{MAD}}\) hardly change, so the robust z-score still flags that observation. R’s default mad(x_hat) multiplies by \(1.4826\) so that the \(\hat{\text{MAD}}\) and the standard deviation are on the same scale for Normal data.

3.3 Shape Statistics

Central tendency and dispersion are often insufficient to describe a distribution. To further describe shape, we can compute sample skew and kurtosis to measure asymmetry and extreme values.

There are many other statistics we could compute on an ad-hoc basis. However, shape is often best understood with graphical descriptions: histogram, ECDF, Boxplot. These should be made before numerical descriptions: skewness and kurtosis statistics.

Skewness

Center and spread do not tell us whether one tail of the distribution is longer than the other.

ImportantKey Definition

The sample skewness measures asymmetry: the average cubed deviation from the mean, scaled by the cube of the standard deviation, \[\hat{\text{skew}} = \frac{\sum_{i=1}^{n}(\hat{x}_{i}-\hat{m})^3 / n}{\hat{s}^3} = \sum_{i=1}^{n}\hat{z}_{i}^3 / n,\] where \(\hat{z}_{i}\) is the standardized value of observation \(i\). Zero skew is symmetric and non-zero skew is asymmetric (positive skew means a long right tail, negative skew a long left tail).

Notice that cubing a deviation keeps its sign. The numerator is therefore close to zero when the data are roughly balanced around the mean, and large in absolute value when one tail is much longer than the other. Dividing by the cubed standard deviation makes skewness dimensionless, so the values are comparable across datasets measured in different units (Joanes and Gill 1998). Skewness is useful for checking whether the mean and median will be close. With low skew they are close, and with high skew they can be far apart. Skewness also helps you decide whether to take the log or square root of the data before further analysis.

The figure below compares two distributions with the same mean and standard deviation. One is symmetric, and the other has a long right tail. On the histogram, the right-skewed distribution piles up just below the mean and stretches far to the right. On the ECDF, it rises steeply at first and then approaches \(1\) slowly.

Consider the right-skewed dataset \(\{0,1,2,3,9\}\) with mean \(\hat{m}=15/5=3\). The cubed deviations from the mean are \((-3)^3, (-2)^3, (-1)^3, 0^3, 6^3\), which equal \(-27, -8, -1, 0, 216\). A cube keeps the sign of its deviation, so the single large positive deviation outweighs the small negative ones. Their average is the numerator: \(\sum_{i=1}^{n}[\hat{x}_{i}-\hat{m}]^3/n = 180/5 = 36\). The variance is \([9+4+1+0+36]/5=10\), so \(\hat{s}=\sqrt{10}\), and dividing by \(\hat{s}^3 \approx 31.6\) gives a skewness of \(36/31.6 \approx 1.14\). The positive value reflects the long right tail created by the value \(9\).

Code 
x_demo_hat <- c(0, 1, 2, 3, 9)
m3_hat <- mean( (x_demo_hat - mean(x_demo_hat))^3 ) # numerator
m3_hat
## [1] 36
s_hat <- sqrt(mean( (x_demo_hat - mean(x_demo_hat))^2 )) # divides by n, unlike sd()
m3_hat / s_hat^3 # skewness
## [1] 1.13842
Code 
## MAIN EXAMPLE
skewness <- function(x_hat){
    z_hat <- standardize(x_hat)
    skew <- mean(z_hat^3) # average cubed standardized value
    return(skew)
}

# skewness of hourly wages
skewness(x_hat)
## [1] 1.970904
# Wages have a large positive skew, about 1.97.
# So their right tail is much longer than their left.
# The mean, about 5.76, sits above the median, about 5.21.

Kurtosis

A second shape question, independent of asymmetry, asks how much of the variability comes from a few far-from-center observations rather than from the typical bulk.

The sample kurtosis measures how heavy the tails are: the average fourth power of deviations from the mean, scaled by \(\hat{s}^4\), \[\hat{\text{kurt}} = \frac{\sum_{i=1}^{n}(\hat{x}_{i}-\hat{m})^4 / n}{\hat{s}^4}.\] The Normal distribution has kurtosis \(3\), so subtracting \(3\) gives a comparison-to-Normal version (excess kurtosis).

Where skewness uses the cube (and so can be negative), kurtosis uses the fourth power, which forces every term to be positive and disproportionately amplifies the largest deviations. Equivalently, kurtosis is the average fourth power of the standardized values, \(\sum_{i=1}^{n}\hat{z}_{i}^4/n\). A few values far from the mean can raise kurtosis a lot, even when most of the distribution looks ordinary. So kurtosis is often described as a “tail-heaviness” measure. Kurtosis is useful for distinguishing distributions that look similar in their bulk but differ in how often extreme values occur: daily stock returns, for example, have moderate variance but kurtosis well above \(3\), reflecting occasional crashes that a Normal model would treat as vanishingly rare.

The figure below compares two symmetric distributions with the same median and \(\hat{\text{IQR}}\). The high-kurtosis distribution has more observations near the center and in the far tails, and fewer in between. On the ECDF, it approaches \(0\) and \(1\) much more slowly in the tails, and is still well above \(0\) and below \(1\) at the edges of the plot.

Using the same dataset \(\{0,1,2,3,9\}\) with mean \(\hat{m}=3\), the quartic deviations are \((-3)^4, (-2)^4, (-1)^4, 0^4, 6^4\), which equal \(81, 16, 1, 0, 1296\). Because the power is even, every term is positive, and the large deviation from the value \(9\) dominates the sum. Their average is the numerator: \(\sum_{i=1}^{n}[\hat{x}_{i}-\hat{m}]^4/n = 1394/5 = 278.8\). Dividing by \(\hat{s}^4 = 10^2 = 100\) gives a kurtosis of \(278.8/100 \approx 2.79\).

Code 
x_demo_hat <- c(0, 1, 2, 3, 9)
m4_hat <- mean( (x_demo_hat - mean(x_demo_hat))^4 ) # numerator
m4_hat
## [1] 278.8
s_hat <- sqrt(mean( (x_demo_hat - mean(x_demo_hat))^2 )) # divides by n, unlike sd()
m4_hat / s_hat^4 # kurtosis
## [1] 2.788

Some authors further subtract \(3\) to explicitly compare against the normal distribution (the normal distribution has a kurtosis of \(3\)).

Boxplot whiskers are a great way to examine kurtosis, with more circle “outliers” indicating more kurtosis. You can also see skew in the boxplot when one quartile is much further from the median than the other.

Code 
## MAIN EXAMPLE
kurtosis <- function(x_hat){
    z_hat <- standardize(x_hat)
    kurt <- mean(z_hat^4) # average fourth power of standardized values
    return(kurt)
    # use instead to compare against normal
    # excess_kurt <- kurt - 3
}

# kurtosis of hourly wages
kurtosis(x_hat)
## [1] 12.63408
# Wages have a kurtosis of about 12.63, well above 3.
# So their tails are much heavier than normal.

Compute both shape statistics for the murder arrest rates in USArrests. How do they compare with the values for the lab data?

Code 
x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000, one value per state

# skewness of murder arrests
skewness(x2_hat)
## [1] 0.3820378
# Murder arrests have a small positive skew.
# So their right tail is slightly longer than their left.

# kurtosis of murder arrests
kurtosis(x2_hat)
## [1] 2.135329
# Murder arrests have a kurtosis below 3.
# So their tails are lighter than normal.

Clusters/Gaps

You can also describe distributions in terms of how clustered the values are, including the number of modes (peaks), bunching, and many other statistics. A gap is a range with few or no observations, and a cluster is a group of values bunched together. A distribution with one peak is unimodal, and one with two peaks is bimodal. The first figure below is bimodal: two clusters of values separated by a gap. The second has a single cluster, with a long flat shelf on its left that ends abruptly at \(-8\). Each pairs a histogram with the ECDF of the same values, where a cluster shows up as a steep rise and a gap as a flatter stretch. The mean and variance alone would miss this structure, since two very different distributions can share the same mean and variance. So remember that “a picture is worth a thousand words”: plot your data before summarizing it numerically.

Use the built-in LifeCycleSavings data. Plot a histogram of pop15 with breaks=15. How many clusters do you see, and where is the gap? Then compute the mean and mark it on the histogram with abline(v=...). How many countries have a value close to the mean?

Code 
# percent of the population under age 15, fifty countries
# Each country's value is averaged over 1960 to 1970.
pop15 <- LifeCycleSavings[, 'pop15']

3.4 Factor Data

The statistics above all add, subtract, or rank the observations, so they require cardinal data. For factor data, we instead start from the proportions \(\hat{p}(x)\) from the empirical mass function and summarize those. Which summaries make sense depends on whether the factor is ordered.

Ordered Factors

Ordered factors, such as letter grades or survey responses from strongly disagree to strongly agree, have a ranking but no meaningful distance between levels. The mode and Herfindahl index defined below also apply, since they only use the proportions. Because the levels can be sorted, the median and other quantiles also apply: the median is the middle level, and the quartiles bracket the middle half of the observations. Quantiles are useful for ordered factors because they only use the ordering: “half the class earned a \(B\) or better” is a statement the data can support. The mean and variance do not apply, because there is nothing to add. A grade point average gets around this by coding the levels as numbers, treating the gap between \(C\) and \(B\) as equal to the gap between \(B\) and \(A\), but that is an assumption imposed by the grader rather than a property of the data. Surveys often record a cardinal variable only in bands, such as income brackets, which turns it into an ordered factor.

For example, the letter grades \(\{B, A, C, B, F, B, A\}\) of \(n=7\) students sort from lowest to highest as \(F, C, B, B, B, A, A\). The median is the fourth sorted value, \(B\), the lower quartile is the second, \(C\), and the upper quartile is the sixth, \(A\). The mode is also \(B\).

Code 
grade <- factor(c('B', 'A', 'C', 'B', 'F', 'B', 'A'),
    levels=c('F', 'D', 'C', 'B', 'A'), ordered=TRUE)
sort(grade)
## [1] F C B B B A A
## Levels: F < D < C < B < A

# median and quartiles (type=1 is required for ordered factors)
quantile(grade, probs=c(.25, .5, .75), type=1)
## 25% 50% 75% 
##   C   B   A 
## Levels: F < D < C < B < A

Compute the median of the original USArrests[, 'UrbanPop'] and check which band it falls in. Then move the inner breaks from \(50\) and \(75\) to \(40\) and \(60\) and recompute the median band. Does the median band always contain the cardinal median?

Below, a cardinal variable is recorded in three bands.

Code 
# urban share of each state's population in three bands, an ordered factor
urban_band <- cut(USArrests[, 'UrbanPop'], breaks=c(0, 50, 75, 100),
    labels=c('mostly rural', 'mixed', 'mostly urban'), ordered_result=TRUE)
p_hat <- table(urban_band)/length(urban_band)
plot(p_hat, col=grey(0, .5),
    xlab='Urban Population', ylab='Proportion of States')

Code 

# mode and Herfindahl index
names(p_hat)[p_hat==max(p_hat)]
## [1] "mixed"
sum(p_hat^2)
## [1] 0.4024

# median and quartiles
quantile(urban_band, probs=c(.25, .5, .75), type=1)
##          25%          50%          75% 
##        mixed        mixed mostly urban 
## Levels: mostly rural < mixed < mostly urban

Unordered Factors

For categories with no natural order, such as regions or brands, the center of a distribution is simply the most common category. Categories have no distance between them, so for factor data spread instead means how evenly the observations divide across the categories.

ImportantKey Definition

The sample mode is the value that occurs most often: the \(x\) with the largest proportion \(\hat{p}(x)\).

The Herfindahl index (Hirschman 1945; Herfindahl 1950) is the sum of the squared proportions across the distinct values, \(\hat{H} = \sum_{x} \hat{p}(x)^2.\)

The mode is useful when you want to know which category is typical. For example, the mode can be the most common region, the best-selling brand, or the most frequent answer on a survey. It is the only measure of center that applies to unordered factors, since the categories can be neither averaged nor ranked. If several categories tie for the highest proportion, the data have more than one mode.

The Herfindahl index ranges from \(1/K\), when every category has the same proportion, up to \(1\), when every observation falls in a single category. It is useful for asking how dominated a distribution is by its largest categories. Ecologists use a close relative, Simpson’s index: the chance that two observations drawn at random without replacement fall in the same category (Simpson 1949). The Herfindahl index is that same chance when the draws are made with replacement, so the two differ little unless \(n\) is small. Economists use the same statistic, with market shares in place of proportions, to measure how concentrated an industry is, and antitrust agencies report it with shares in percentage points so that the scale runs to \(10{,}000\). The complement \(1-\hat{H}\) is sometimes reported instead so that larger values mean a more even spread, and after rescaling by \(K/(K-1)\) it is called the index of qualitative variation.

For example, the dataset \(\{A, B, A, C, C, A\}\) has \(n=6\) observations in \(K=3\) categories. \(A\) appears three times, \(C\) twice, and \(B\) once, so \(\hat{p}(A)=1/2\), \(\hat{p}(B)=1/6\), \(\hat{p}(C)=1/3\), and the mode is \(A\). Notice the mode is not the alphabetically middle letter, \(B\). The Herfindahl index is \((1/2)^2+(1/6)^2+(1/3)^2=[9+1+4]/36=14/36\approx 0.39\), above the minimum of \(1/3\) that would occur if all three letters were equally common.

Code 
x0_hat <- c('A', 'B', 'A', 'C', 'C', 'A')
p_hat <- table(x0_hat)/length(x0_hat)
p_hat
## x0_hat
##         A         B         C 
## 0.5000000 0.1666667 0.3333333

# mode(s)
names(p_hat)[p_hat==max(p_hat)]
## [1] "A"

# Herfindahl index
sum(p_hat^2)
## [1] 0.3888889
Code 
## MAIN EXAMPLE
# census region of each state, an unordered factor
region <- state.region
p_hat <- table(region)/length(region)
plot(p_hat, col=grey(0, .5),
    xlab='Census Region', ylab='Proportion of States')

Code 

# mode
names(p_hat)[p_hat==max(p_hat)]
## [1] "South"

# Herfindahl index, compared to its minimum 1/K
K <- length(p_hat)
c(H_hat=sum(p_hat^2), H_min=1/K)
## H_hat H_min 
##  0.26  0.25

Suppose \(100\) customers each buy from one of four firms: \(60\) from firm \(A\), \(20\) from \(B\), \(10\) from \(C\), and \(10\) from \(D\). Compute each firm’s market share, the mode, and the Herfindahl index. Then suppose firms \(C\) and \(D\) merge into a single firm and recompute. Does the merger raise or lower concentration?

Code 
x0_hat <- rep(c('A', 'B', 'C', 'D'), times=c(60, 20, 10, 10))
p_hat <- table(x0_hat)/length(x0_hat)
sum(p_hat^2)
## [1] 0.42

# after the merger, C and D are one firm
x_merged_hat <- x0_hat
x_merged_hat[x_merged_hat=='D'] <- 'C'

Each statistic in this chapter suits some data types and not others. Here is a table to help you recall which apply. Note that the mode and Herfindahl index need repeated values, so for continuous cardinal data we can describe proportions with a histogram instead (seeing the highest peak as the mode, for example).

Summary Statistics by Data Type
Statistic Factor Cardinal
Mean, variance, skewness, kurtosis, \(\hat{\text{IQR}}\) No Yes
Median, quantiles Ordered only (not Unordered) Yes
Mode, Herfindahl index Yes Discrete only (not Continuous)

3.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. Explain in your own words why the median is more robust to extreme values than the mean. Give a concrete example with a small dataset where adding one outlier changes the mean substantially but barely moves the median.

  3. Using the dataset \(\{2, 5, 5, 8, 12\}\), compute by hand: (a) the sample mean \(\hat{m}\), (b) the sample variance \(\hat{v}\), (c) the sample standard deviation \(\hat{s}\), and (d) the \(\hat{\text{IQR}}\). Show your work for each step.

  4. Load USArrests in R. Compute the mean, standard deviation, and MAD (with constant=1) of the Assault variable. Then apply a log transformation with log(USArrests[,'Assault']) and compute the skewness of both the original and transformed data using the skewness function defined in the chapter. Which version is less skewed?

  5. A survey asks \(12\) customers to rate a product on a five-level scale from poor to excellent, with responses \(\{\text{good}, \text{fair}, \text{excellent}, \text{good}, \text{poor}, \text{good}, \text{fair}, \text{excellent}, \text{good}, \text{very good}, \text{fair}, \text{good}\}\). Store the responses as an ordered factor in R, then report the proportions, the mode, the median, and the Herfindahl index. Explain why the median is a meaningful summary here but the mean of the numeric codes \(1\) to \(5\) requires an extra assumption.

Further Reading

Recall

This chapter summarized a distribution with single numbers. The mean, variance, and standard deviation describe center and spread. Standardized values measure distance from the mean in standard deviations. The median, IQR, and MAD are robust alternatives, and the skewness and kurtosis describe shape. For factor data, nothing can be added, so we used the mode and the Herfindahl index \(\hat{H}\). The four census regions in state.region have mode South and \(\hat{H}=0.26\). That is barely above the \(0.25\) floor for four equally common categories. Ordered factors, such as letter grades, also have a median and quartiles. We used the five-point dataset \(\{0, 1, 2, 3, 9\}\) for the shape statistics. Its mean is \(3\), and its cubed deviations average to \(36\), giving a skewness of about \(1.14\) and a long right tail. Its quartic deviations give a kurtosis of about \(2.79\). The single outlier at \(9\) drives both numbers. For real data such as Wages1[, 'wage'], a histogram with the mean and standard deviation already shows most of the shape. The median, MAD, skewness, and kurtosis add detail when the data are not symmetric.

Chebyshev, P. L. 1867. “Des valeurs moyennes.” Journal de Mathématiques Pures et Appliquées, 2e série 12: 177–84. https://www.numdam.org/item/JMPA_1867_2_12__177_0/.
Hampel, Frank R. 1974. “The Influence Curve and Its Role in Robust Estimation.” Journal of the American Statistical Association 69 (346): 383–93. https://doi.org/10.1080/01621459.1974.10482962.
Herfindahl, Orris C. 1950. “Concentration in the Steel Industry.” PhD thesis, Columbia University.
Hirschman, Albert O. 1945. National Power and the Structure of Foreign Trade. University of California Press. https://doi.org/10.1525/9780520378179.
Joanes, D. N., and C. A. Gill. 1998. “Comparing Measures of Sample Skewness and Kurtosis.” Journal of the Royal Statistical Society: Series D (The Statistician) 47 (1): 183–89. https://doi.org/10.1111/1467-9884.00122.
Rousseeuw, Peter J., and Christophe Croux. 1993. “Alternatives to the Median Absolute Deviation.” Journal of the American Statistical Association 88 (424): 1273–83. https://doi.org/10.1080/01621459.1993.10476408.
Simpson, E. H. 1949. “Measurement of Diversity.” Nature 163 (4148): 688. https://doi.org/10.1038/163688a0.
Tukey, John Wilder. 1977. Exploratory Data Analysis. Vol. 2. Springer.