Earlier chapters summarized one variable at a time. We now study two variables together. This chapter builds the two basic ways to summarize them: the joint distribution, which records how often each combination of values occurs, and the conditional distribution, which records the distribution of one variable within fixed values of the other. We do each for discrete data and then for continuous data. Throughout, the data for each observation are grouped together as a vector \((\hat{x}_{i}, \hat{y}_{i})\).
# USArrests has 50 rows (one per US state) and four crime columns;
# pick two of those columns to form a bivariate dataset.
xy_hat <- USArrests[, c('UrbanPop', 'Murder')]
colnames(xy_hat) <- c('x_hat', 'y_hat') # x_hat=urban population share, y_hat=murder arrests per 100K
xy_hat[1, ] # first state's two values
## x_hat y_hat
## Alabama 58 13.2When two variables are recorded together for each observation, the first question is how often each combination of values appears.
The joint distribution is also called a frequency table (or cross-tabulation) and is useful for summarizing how two categorical variables co-occur: it is the starting point for every other two-variable statistic in this chapter. Dividing the count by \(n\) puts the cells on a probability scale, so each cell reports a proportion of the sample and the whole table sums to one.
Sometimes we only need to summarize one of the two variables and ignore the other.
The marginal (so called because it lives in the row and column margins of the joint table) is useful for recovering the univariate picture once the joint table is built, and for supplying the base rates (\(\hat{p}(x)\), \(\hat{p}(y)\)) that conditional distributions divide by. The same numbers can also be computed directly from the univariate data without ever forming the joint table (see here).
For example, consider the Wages1 data on US workers. Let \(\hat{x}_{i}\) be years of schooling and \(\hat{y}_{i}\) be sex. We keep only workers with at least \(10\) years of schooling, which drops the few workers with less schooling and leaves \(n=3001\).
library(Ecdat)
# Workers with at least 10 years of schooling: schooling and sex
xy2_hat <- Wages1[, c('school', 'sex')]
xy2_hat <- xy2_hat[xy2_hat[, 'school'] >= 10, ]
# Frequency table and joint distribution
tab <- table(xy2_hat) # raw counts: rows=school, cols=sex
tab
## sex
## school female male
## 10 166 233
## 11 307 354
## 12 594 594
## 13 189 150
## 14 140 124
## 15 68 66
## 16 11 5
prop <- tab / sum(tab) # joint sample distribution (sums to 1)
round(prop, 3)
## sex
## school female male
## 10 0.055 0.078
## 11 0.102 0.118
## 12 0.198 0.198
## 13 0.063 0.050
## 14 0.047 0.041
## 15 0.023 0.022
## 16 0.004 0.002
# addmargins() appends row and column totals (i.e., the marginals)
prop_m <- addmargins(prop)
round(prop_m, 3)
## sex
## school female male Sum
## 10 0.055 0.078 0.133
## 11 0.102 0.118 0.220
## 12 0.198 0.198 0.396
## 13 0.063 0.050 0.113
## 14 0.047 0.041 0.088
## 15 0.023 0.022 0.045
## 16 0.004 0.002 0.005
## Sum 0.492 0.508 1.000Each cell of the joint table is the share of all workers with that schooling and sex. For example, \(\hat{p}(x=12, y=\text{female}) \approx 0.198\), so about \(20\%\) of workers are women with \(12\) years of schooling. The row totals are the marginal distribution of schooling: \(\hat{p}(x=12) \approx 0.396\). The column totals are the marginal distribution of sex: \(\hat{p}(y=\text{female}) \approx 0.49\).
Note that you can compute various univariate statistics using the marginal distributions. E.g., the mean of \(\hat{x}_{i}\) is \(\hat{m}_{X}=\sum_{x} x \hat{p}(x)\). In the above example, that equals about \(12.0\) years of schooling. Try computing the standard deviation \(\hat{s}_{X}\) yourself.
# marginals and univariate statistics
p_x <- margin.table(prop, 1)
sum(p_x * as.numeric(names(p_x)))
## [1] 11.95801
mean(xy2_hat[, 'school'])
## [1] 11.95801
# p_y <- margin.table(prop, 2)Scatterplots are used frequently to summarize the joint distribution of continuous data. They can be enhanced in several ways. As a default, use semi-transparent points so as not to hide any points (and perhaps see if your observations are concentrated anywhere). You can also add other features that help summarize the relationship.
plot(y_hat ~ x_hat, xy_hat, pch=16, col=grey(0, .5),
main=NA, xlab='Urban Population', ylab='Murder Arrests')
You can also show the marginal distributions of each variable along each axis (see this margin and layout cheatsheet for how layout(), mar, and oma work together).
# Setup Plot
layout( matrix(c(2, 0, 1, 3), ncol=2, byrow=TRUE),
widths=c(9/10, 1/10), heights=c(1/10, 9/10))
# Scatterplot
par(mar=c(4, 4, 1, 1))
plot(y_hat ~ x_hat, xy_hat, pch=16, col=grey(0, .5),
main=NA, xlab='Urban Population', ylab='Murder Arrests')
# Add Marginals
par(mar=c(0, 4, 1, 1))
f_x_hat <- hist(xy_hat[, 'x_hat'], plot=FALSE)
barplot(f_x_hat[['counts']], axes=FALSE, space=0, border=NA)
par(mar=c(4, 0, 1, 1))
f_y_hat <- hist(xy_hat[, 'y_hat'], plot=FALSE)
barplot(f_y_hat[['counts']], axes=FALSE, space=0, horiz=TRUE, border=NA)
Often we want to know how the distribution of one variable changes when we restrict attention to a subgroup defined by another.
Conditional distributions are useful when one variable splits the data into meaningful subgroups (sex, education level, treatment status) and we want to compare the outcome distribution across those subgroups. The renormalization is what changes the unit of analysis: row sums equal to one means each row is itself a probability distribution over \(y\) for the given \(x\), rather than a share of the whole sample.
For example, return to the Wages1 workers and their joint distribution table. To compute the conditional distribution of \(\hat{y}_{i}\) given \(\hat{x}_{i}=x\), we rescale the numbers in each row so that each row sums to 1.
# Conditional distribution of sex given schooling: each row sums to 1
round(prop.table(tab, 1), 3)
## sex
## school female male
## 10 0.416 0.584
## 11 0.464 0.536
## 12 0.500 0.500
## 13 0.558 0.442
## 14 0.530 0.470
## 15 0.507 0.493
## 16 0.688 0.312In this example, we say that
We can also compute the conditional distribution of \(\hat{x}_{i}\) given \(\hat{y}_{i}=y\), \(\hat{p}(x \mid y) = \hat{p}(x,y) / \hat{p}(y)\), just as we did above. In this case, we rescale the numbers in each column of the joint distribution table so that each column sums to 1.
# Conditional distribution of schooling given sex: each column sums to 1
round(prop.table(tab, 2), 3)
## sex
## school female male
## 10 0.113 0.153
## 11 0.208 0.232
## 12 0.403 0.389
## 13 0.128 0.098
## 14 0.095 0.081
## 15 0.046 0.043
## 16 0.007 0.003For example, among female workers, about \(40\%\) have \(12\) years of schooling.
Different conditional distributions change the unit of analysis. Swapping which variable is in the denominator like this generally answers a different question.
Confusing the different distributions is so common that they have names.
The base-rate fallacy often shows up as reading a large count as a high rate. A common category can have the most cases even when its rate per case is the lowest.
How far a signal moves \(\hat{p}(y\mid x)\) away from \(\hat{p}(y)\) depends on the base rate \(\hat{p}(y)\) itself. The fraud-detection example below applies one screen to a rare outcome and then to a common one.
These describe the relationship between \(\hat{y}_{i}\) and \(\hat{x}_{i}\). We show how \(Y\) changes according to \(X\) using a histogram or bar plot. When \(X\) is continuous, as it often is, we split it into distinct bins and convert it to a factor variable. E.g.,
# Split Data by Urban Population above/below mean
m_x_hat <- mean(xy_hat[, 'x_hat'])
low_urban <- xy_hat[, 'x_hat'] <= m_x_hat
y1_hat <- xy_hat[low_urban, 'y_hat'] # murder arrests, low urban population
y2_hat <- xy_hat[!low_urban, 'y_hat'] # murder arrests, high urban population
cols <- c(low=rgb(0, 0, 1, .75), high=rgb(1, 0, 0, .75))
# Common Histogram
ylim <- c(0, .25)
xbks <- seq(min(xy_hat[, 'y_hat'])-1, max(xy_hat[, 'y_hat'])+1, by=1)
par(mfrow=c(1, 2))
hist(y1_hat,
breaks=xbks, col=cols[1],
main=NA,
xlab='Murder Arrests', freq=FALSE,
border=NA, ylim=ylim)
title('Urban Pop <= Mean', font.main=1)
hist(y2_hat,
breaks=xbks, col=cols[2],
main=NA,
xlab='Murder Arrests', freq=FALSE,
border=NA, ylim=ylim)
title('Urban Pop > Mean', font.main=1)
It is sometimes preferable to show the ECDF instead. And you can glue various combinations together to convey more information all at once
layout( t(c(1, 2, 2)))
# Full Sample Density
hist(xy_hat[, 'y_hat'],
main=NA,
xlab='Murder Arrests',
breaks=xbks, freq=FALSE, border=NA)
title('Full Sample Density', font.main=1)
# Split Sample Distribution Comparison
F_lowpop_hat <- ecdf(y1_hat)
plot(F_lowpop_hat, col=cols[1],
pch=16, xlab='Murder Arrests',
main=NA, bty='n')
title('Split Sample Distributions', font.main=1)
F_highpop_hat <- ecdf(y2_hat)
plot(F_highpop_hat, add=TRUE, col=cols[2], pch=16)
legend('bottomright', col=cols,
pch=16, bty='n', inset=c(0, .1),
title='% Urban Pop.',
legend=c('Low (<= Mean)', 'High (> Mean)'))
You can also split data into more than two groups. For more than three groups, boxplots are often more effective than histograms or ECDF’s.
# K Groups with even spacing (not even counts)
K <- 4
x_cut_hat <- cut(xy_hat[, 'x_hat'], K)
table(x_cut_hat)
## x_cut_hat
## (31.9,46.8] (46.8,61.5] (61.5,76.2] (76.2,91.1]
## 6 13 17 14
# Boxplots for each group
Kcols <- hcl.colors(K, alpha=.5)
boxplot(xy_hat[, 'y_hat'] ~ x_cut_hat,
main=NA, col=Kcols,
whisklty=0, staplelty=0, outline=FALSE,
varwidth=TRUE, #show number of obs. per group
xlab='Urban Population', ylab='Murder Arrests')
# 4 Groups with equal numbers of observations
#Qcuts <- c(
# '0%'=min(xy_hat[, 'x_hat'])-10*.Machine[['double.eps']],
# quantile(xy_hat[, 'x_hat'], probs=c(.25, .5, .75, 1)))
#x_qcut_hat <- cut(xy_hat[, 'x_hat'], Qcuts)
#boxplot(xy_hat[, 'y_hat'] ~ x_qcut_hat, col=hcl.colors(4, alpha=.5))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.
Explain the difference between a joint distribution, a marginal distribution, and a conditional distribution. Using the schooling-and-sex example from this chapter, describe in words what \(\hat{p}(y=\text{female} \mid x=13)\) tells you and how it differs from \(\hat{p}(x=13 \mid y=\text{female})\).
Suppose you observe a sample of \(n=20\) people with two discrete variables: \(\hat{x}_{i}\) takes values in \(\{A, B\}\) and \(\hat{y}_{i}\) takes values in \(\{1, 2\}\). The counts are: \((A,1)=6\), \((A,2)=4\), \((B,1)=3\), \((B,2)=7\). Compute the joint distribution \(\hat{p}(x,y)\), both marginal distributions, and the conditional distribution \(\hat{p}(y \mid x)\).
Using the USArrests dataset, split the variable Assault into two groups based on whether UrbanPop is above or below its median. Plot overlapping histograms of Assault for both groups with semi-transparent colors. Then plot the ECDF for each group on the same axes with a legend.
This chapter built up two-variable summaries from joint to marginal to conditional distributions, illustrated with the Wages1 schooling-and-sex table where \(\hat{p}(y=\text{female}\mid x=13)\approx 0.56\) is the share of female workers within the 13-years-of-school subgroup. We then saw how confusing the direction of conditioning produces the base-rate fallacy (the “most-stolen car” calculation, where Civics are the most-stolen model but the least likely to be stolen per car).