Skip to main content

14  Statistics of Association


There are several ways to statistically assess the relationship between two variables. The major differences surround whether the data are cardinal or an ordered/unordered factor. In all cases, these statistics measure association, not causation.

This chapter asks three questions.

14.1 Cardinal Data

Pearson (Linear) Correlation

When two variables are cardinal, the foundational measure of how they move together is the covariance.

ImportantKey Definition

The covariance \(\hat{c}_{XY}\) is the average product of paired deviations from each variable’s mean (where \(\hat{m}_{X}\) and \(\hat{m}_{Y}\) are the sample means): \[\hat{c}_{XY} = \sum_{i=1}^{n} [\hat{x}_{i} - \hat{m}_{X}] [\hat{y}_{i} - \hat{m}_{Y}] / n.\]

The covariance is useful as a building block for other linear two-variable statistics. It is positive when \(\hat{x}_{i}\) and \(\hat{y}_{i}\) tend to be above or below their means together. It is negative when they tend to move apart, and zero when there is no linear co-movement. Its units are the units of \(\hat{x}_{i}\) times the units of \(\hat{y}_{i}\), which makes raw values hard to compare across datasets. Note that the covariance of \(\hat{x}\) with itself is just its variance, \(\hat{c}_{XX}=\hat{v}_{X}\), and the standard deviation is \(\hat{s}_{X}=\sqrt{\hat{v}_{X}}\).

For an interpretable, unit-free measure that is comparable across datasets, we rescale the covariance.

ImportantKey Definition

The Pearson correlation \(\hat{r}_{XY}\) (Pearson 1896) rescales the covariance by the product of the two standard deviations: \[\hat{r}_{XY} = \frac{ \hat{c}_{XY} }{ \hat{s}_{X} \hat{s}_{Y}}.\]

The Pearson correlation is useful as the standard unit-free summary of linear association. It is constrained to \([-1, 1]\), with \(+1\) for a perfect rising line, \(-1\) for a perfect falling line, and \(0\) for no linear association. Because the bounds are the same in every dataset, it lets us compare the strength of two-variable associations across studies.

What is the correlation for the dataset \(\{ (0,0.1) , (1, 0.3), (2, 0.2) \}\)? Find the answer both mathematically and computationally.

Mathematically, there are five steps.

Step 1: Compute the means \[\begin{aligned} \hat{m}_{X} &= \frac{0+1+2}{3} = 1 \\ \hat{m}_{Y} &= \frac{0.1+0.3+0.2}{3} = 0.2 \end{aligned}\]

Step 2: Compute the deviances \[ \begin{array}{c|rrrr} \hat{x}_{i} & 0 & 1 & 2 \\ \hat{x}_{i}-\hat{m}_{X} & -1 & 0 & 1 \\ \hat{y}_{i} & 0.1 & 0.3 & 0.2 \\ \hat{y}_{i}-\hat{m}_{Y} & -0.1 & 0.1 & 0 \end{array} \]

Step 3: Compute the Covariance \[\begin{aligned} \hat{c}_{XY} &= \sum (\hat{x}_{i}-\hat{m}_{X})(\hat{y}_{i}-\hat{m}_{Y})/n = \left[ (-1)(-0.1) + 0(0.1) + 1(0) \right] \frac{1}{3} = (0.1) \frac{1}{3} = 1/30 \end{aligned}\]

Step 4: Compute Standard Deviations \[\begin{aligned} \hat{v}_{X} &= \sum_{i=1}^n \left(\hat{x}_{i}-\hat{m}_{X}\right)^2 / n = \left[(-1)^2+0^2+1^2 \right]/3 = 2/3 \\ \hat{s}_{X} &= \sqrt{2/3} \\ \hat{v}_{Y} &= \sum_{i=1}^n \left( \hat{y}_{i}-\hat{m}_{Y} \right)^2 / n = \left[ (-0.1)^2+(0.1)^2+0^2 \right]/3 = \left[0.01+0.01\right]/3 = \frac{2}{100} \frac{1}{3} = 2/300 \\ \hat{s}_{Y} &= \sqrt{2/300} \end{aligned}\]

Step 5: Compute the Correlation \[\begin{aligned} \frac{\hat{c}_{XY}}{\hat{s}_{X} \hat{s}_{Y}} &= \frac{1/30}{ \sqrt{2/3} \sqrt{2/300}} = \frac{1/30}{ 2 /\sqrt{900}} = \frac{1/30}{2/30} = 1/2 \end{aligned}\]

Note that this value suggests a positive relationship between the variables.

Computationally, we do the same steps

Code
# Create the Data
x0_hat <- c(0, 1, 2)
x0_hat
## [1] 0 1 2
y0_hat <- c(0.1, 0.3, 0.2)
y0_hat
## [1] 0.1 0.3 0.2

# Compute the Means
m_x_hat <- mean(x0_hat)
m_y_hat <- mean(y0_hat)

# Compute the Deviances
dev_X <- x0_hat - m_x_hat
dev_Y <- y0_hat - m_y_hat

# Compute the Covariance
cov_manual <-  sum(dev_X * dev_Y) / length(x0_hat)

# Compute the Standard Deviations
v_x_hat <- sum(dev_X^2) / length(x0_hat)
s_x_hat <- sqrt(v_x_hat)
v_y_hat <- sum(dev_Y^2) / length(y0_hat)
s_y_hat <- sqrt(v_y_hat)

# Compute the Correlation
cor_manual <- cov_manual / (s_x_hat * s_y_hat)
cor_manual
## [1] 0.5

# Verify with the built-in function
cor(x0_hat, y0_hat)
## [1] 0.5

For example, the code below draws a scatterplot of two cardinal variables and computes their covariance and correlation.

Code
## MAIN EXAMPLE
# Load the Data: 3294 US workers, with hourly wages in 1980 dollars
data(Wages1, package='Ecdat')
xy_hat <- Wages1[, c('school', 'wage')]
colnames(xy_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=hourly wage

# Scatterplot
plot(y_hat ~ x_hat, data=xy_hat, pch=16, col=grey(0, .1),
    main=NA, xlab='Years of schooling', ylab='Hourly wage (1980 dollars)')

Code
# Workers with more schooling tend to earn more.
# Wages vary a lot at every level of schooling.

# Covariance: cov() divides by n-1 rather than n, which barely matters with n=3294
cov(xy_hat[, 'x_hat'], xy_hat[, 'y_hat'])
## [1] 1.530776
# The covariance is positive.
# Its units (years times dollars per hour) are hard to judge.

# Correlation
xy_cor <- cor(xy_hat[, 'x_hat'], xy_hat[, 'y_hat'])
xy_cor
## [1] 0.2824925
# The correlation of about 0.28 suggests a positive but weak linear association.

Hypothesis Testing

You can conduct hypothesis tests for these statistics using the same procedures we learned for univariate data, Intervals. For example, you can do so by inverting a confidence interval. The histogram below shows the bootstrap distribution of the correlation, with the interval and \(0\) marked.

Code
## MAIN EXAMPLE
# xy_hat and xy_cor defined above
# Bootstrap Distribution of Correlation
n <- nrow(xy_hat)
bootstrap_cor <- rep(NA, 9999)
for(b in seq(bootstrap_cor) ){
    xy_boot <- xy_hat[sample(n, replace=TRUE), ]
    xy_cor_b <- cor(xy_boot[, 'x_hat'], xy_boot[, 'y_hat'])
    bootstrap_cor[b] <- xy_cor_b
}
hist(bootstrap_cor, breaks=100,
    border=NA, freq=FALSE,
    xlim=range(c(0, bootstrap_cor)),
    xlab='Correlation',
    main=NA)
title('Bootstrap Distribution', font.main=1)

## Test whether correlation is statistically different from 0
boot_ci <- quantile(bootstrap_cor, probs=c(0.025, 0.975))
abline(v=boot_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=0, col=rgb(0, 0, 1, .8), lty=2)

Code
boot_ci
##      2.5%     97.5% 
## 0.2493675 0.3154663
# The interval, from about 0.25 to 0.32, does not contain 0.

Importantly, we can also impose the null hypothesis of no association by reshuffling the data. If we resample without replacement, this is known as a permutation test.

Code
## MAIN EXAMPLE
# xy_hat and xy_cor defined above
# Null Permutation Distribution of Correlation
n <- nrow(xy_hat)
null_permutation_cor <- rep(NA, 9999)
for(b in seq(null_permutation_cor) ){
    xy_perm <- xy_hat
    xy_perm[, 'x_hat'] <- xy_hat[sample(n, replace=FALSE), 'x_hat'] ## Reshuffle X
    xy_cor_b <- cor(xy_perm[, 'x_hat'], xy_perm[, 'y_hat'])
    null_permutation_cor[b] <- xy_cor_b
}
hist(null_permutation_cor, breaks=100,
    border=NA, freq=FALSE,
    xlim=range(c(null_permutation_cor, xy_cor)),
    xlab='Correlation',
    main=NA)
title('Null Permutation Distribution', font.main=1)

## Test whether correlation is statistically different from 0
perm_ci <- quantile(null_permutation_cor, probs=c(0.025, 0.975))
abline(v=perm_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=xy_cor, col=rgb(1, 0, 0, .8), lwd=2)

Code
perm_ci
##        2.5%       97.5% 
## -0.03371363  0.03448983
# The shuffled correlations stay within a few hundredths of 0.
# So the observed correlation of about 0.28 is far outside them.

Any association between \(\hat{x}\) and \(\hat{y}\) comes from how their values are paired. Shuffling \(\hat{x}\) keeps its values but breaks the pairing.

Edit the code above to use years of experience (exper) in place of schooling. The code below computes that correlation. Is it statistically different from \(0\)?

Code
# Wages1 loaded above
cor(Wages1[, 'exper'], Wages1[, 'wage'])
## [1] 0.05138051
# The correlation of about 0.05 is much smaller than the 0.28 for schooling.

Repeat both tests for the correlation between urban population share and murder arrests across US states. Use one random sample of \(35\) of the \(50\) states, as in Comparing Multiple Groups.

Code
# USArrests covers all 50 states, so it is the whole population of states.
# To make the bootstrap meaningful, pretend we observed only a random sample.
# Draw one random sample of 35 of the 50 states, and use only it from here on
state_id <- sample(nrow(USArrests), 35)

# Load the Data for the sampled states
xy2_hat <- USArrests[state_id, c('UrbanPop', 'Murder')]
colnames(xy2_hat) <- c('x_hat', 'y_hat') # x_hat=urban population share, y_hat=murder arrests per 100K
xy2_cor <- cor(xy2_hat[, 'x_hat'], xy2_hat[, 'y_hat'])
xy2_cor
## [1] 0.1618443
# The correlation is about 0.16, a weak association.

# Bootstrap Distribution of Correlation
n2 <- nrow(xy2_hat)
bootstrap_cor2 <- rep(NA, 9999)
for(b in seq(bootstrap_cor2) ){
    xy2_boot <- xy2_hat[sample(n2, replace=TRUE), ]
    xy2_cor_b <- cor(xy2_boot[, 'x_hat'], xy2_boot[, 'y_hat'])
    bootstrap_cor2[b] <- xy2_cor_b
}
hist(bootstrap_cor2, breaks=100,
    border=NA, freq=FALSE,
    xlab='Correlation',
    main=NA)
title('Bootstrap Distribution', font.main=1)

## Test whether correlation is statistically different from 0
boot_ci2 <- quantile(bootstrap_cor2, probs=c(0.025, 0.975))
abline(v=boot_ci2, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=0, col=rgb(0, 0, 1, .8), lty=2)

Code
boot_ci2
##       2.5%      97.5% 
## -0.1783681  0.4866896
# The interval, from about -0.18 to 0.49, contains 0.

# Null Permutation Distribution of Correlation
null_permutation_cor2 <- rep(NA, 9999)
for(b in seq(null_permutation_cor2) ){
    xy2_perm <- xy2_hat
    xy2_perm[, 'x_hat'] <- xy2_hat[sample(n2, replace=FALSE), 'x_hat'] ## Reshuffle X
    xy2_cor_b <- cor(xy2_perm[, 'x_hat'], xy2_perm[, 'y_hat'])
    null_permutation_cor2[b] <- xy2_cor_b
}
hist(null_permutation_cor2, breaks=100,
    border=NA, freq=FALSE,
    xlab='Correlation',
    main=NA)
title('Null Permutation Distribution', font.main=1)

## Test whether correlation is statistically different from 0
perm_ci2 <- quantile(null_permutation_cor2, probs=c(0.025, 0.975))
abline(v=perm_ci2, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=xy2_cor, col=rgb(1, 0, 0, .8), lwd=2)

Code
perm_ci2
##       2.5%      97.5% 
## -0.3375322  0.3356207
# The observed correlation is well inside the middle 95% of shuffled values.

So the bootstrap evaluates sampling variability in the world that generated your data. A permutation test constructs the distribution of the statistic in a world where the null is true. Altogether, we have

Types of resampling
Distribution Sample Size per Iteration Number of Iterations Mechanism Typical Purpose
Jackknife \(n-1\) \(n\) \(\hat{x},\hat{y}\): Deterministically leave-one-out observation Variance estimate (after rescaling) and Normal CI estimate
Bootstrap \(n\) \(B\) \(\hat{x},\hat{y}\): Random resample with replacement Percentile CI estimate
Null Bootstrap \(n\) \(B\) \(\hat{x},\hat{y}\): Random resample with replacement and shifted Percentile CI under imposed null, \(p\)-values
Permutation \(n\) \(B\) \(\hat{x}\): Random resample without replacement Percentile CI under imposed null of no association, \(p\)-values
Code
# xy_hat and xy_cor defined above
# Null Bootstrap Distribution of Correlation
n <- nrow(xy_hat)
null_bootstrap_cor <- rep(NA, 9999)
for(b in seq(null_bootstrap_cor) ){
    xy_boot <- xy_hat
    xy_boot[, 'x_hat'] <- xy_hat[sample(n, replace=TRUE), 'x_hat'] ## Reshuffle X
    xy_cor_b <- cor(xy_boot[, 'x_hat'], xy_boot[, 'y_hat'])
    null_bootstrap_cor[b] <- xy_cor_b
}
hist(null_bootstrap_cor, breaks=100,
    border=NA, freq=FALSE,
    xlim=range(c(null_bootstrap_cor, xy_cor)),
    xlab='Correlation',
    main=NA)
title('Null Bootstrap Distribution', font.main=1)

## Test whether correlation is statistically different from 0
boot_ci <- quantile(null_bootstrap_cor, probs=c(0.025, 0.975))
abline(v=boot_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=xy_cor, col=rgb(1, 0, 0, .8), lwd=2)

Code
# The observed correlation of about 0.28 is far outside this distribution too.

Drawing \(\hat{x}\) with replacement also breaks the pairing, which gives the Null Bootstrap row of the table. Compare this distribution to the permutation distribution above.

Robust Correlation

When the data contain outliers, the mean-based covariance can be misleading and a robust alternative is needed. The Falk codeviance \(\tilde{c}_{XY}\) (Falk 1997) replaces the mean-based average in the covariance with the median of paired deviations from each variable’s median (where \(\tilde{m}_{X}\) and \(\tilde{m}_{Y}\) are the sample medians): \[\tilde{c}_{XY} = \text{Med}\left\{ (\hat{x}_{i} - \tilde{m}_{X})(\hat{y}_{i} - \tilde{m}_{Y}) \right\}.\]

Falk’s median correlation \(\tilde{r}_{XY}\) rescales the codeviance by the product of the two median absolute deviations: \[\tilde{r}_{XY} = \frac{ \tilde{c}_{XY} }{ \hat{\text{MAD}}_{X} \hat{\text{MAD}}_{Y}}.\]

The Falk codeviance is useful for cardinal data with outliers or heavy tails, such as income and prices. One extreme observation can move \(\hat{c}_{XY}\) by a lot but \(\tilde{c}_{XY}\) by very little. The median, \(\hat{\text{IQR}}\), and \(\hat{\text{MAD}}\) from Part 1 are robust for the same reason.1 Falk’s median correlation \(\tilde{r}_{XY}\) typically lies in \([-1,1]\) but, unlike the Pearson correlation, not always. Ties can be a problem, however. Every observation tied at a median has a product of \(0\), so with enough ties the median product is also \(0\).

Code
## MAIN EXAMPLE
codev <- function(xy) {
  # Compute medians for each column
  med_hat <- apply(xy, 2, median)
  # Subtract the medians from each column
  xm <- sweep(xy, 2, med_hat, '-')
  # Return the median of the paired products
  CoDev <- median(xm[, 1] * xm[, 2])
  return(CoDev)
}
med_cor <- function(xy) {
  # Compute the medians of absolute deviation
  med_hat <- apply(xy, 2, median)
  xm <- sweep(xy, 2, med_hat, '-')
  MadProd <- prod( apply(abs(xm), 2, median) )
  # Return the robust correlation measure
  return( codev(xy) / MadProd)
}
# xy_hat defined above
xy_codev <- codev(xy_hat)
xy_codev
## [1] 0
xy_med_cor <- med_cor(xy_hat)
xy_med_cor
## [1] 0
# Both are exactly 0, although the Pearson correlation is about 0.28.
# The median worker has 12 years of schooling, and so do 36% of all workers.
# Their products are all 0, and the median product falls among them.

Compute the Codeviance for the dataset \(\{(1,2),(2,1),(3,4),(4,3),(5,6)\}\).

The medians are \(\tilde{m}_{X}=3\) and \(\tilde{m}_{Y}=3\), giving signed deviations and products

\[\begin{array}{c|rrrrr} \hat{x}_{i}-\tilde{m}_{X} & -2 & -1 & 0 & 1 & 2 \\ \hat{y}_{i}-\tilde{m}_{Y} & -1 & -2 & 1 & 0 & 3 \\ (\hat{x}_{i}-\tilde{m}_{X})(\hat{y}_{i}-\tilde{m}_{Y}) & 2 & 2 & 0 & 0 & 6 \end{array}\]

The Codeviance is the median of those products: \(\tilde{c}_{XY}=\text{Med}\{0,0,2,2,6\}=2\). The median absolute deviations are \(\hat{\text{MAD}}_{X}=\text{Med}\{0,1,1,2,2\}=1\) and \(\hat{\text{MAD}}_{Y}=\text{Med}\{0,1,1,2,3\}=1\), so Falk’s median correlation is \(\tilde{r}_{XY}=2/(1\cdot1)=2\). Here \(\tilde{r}_{XY}\) exceeds \(1\): unlike Pearson’s correlation, Falk’s median correlation is not confined to \([-1,1]\).

Code
xy0_hat <- cbind(x_hat=c(1,2,3,4,5), y_hat=c(2,1,4,3,6))
codev(xy0_hat)
## [1] 2
med_cor(xy0_hat)
## [1] 2

You construct sampling distributions and conduct hypothesis tests for Falk’s codeviance and median correlation in the same way you do for Pearson’s Correlation statistic.

Shuffling a variable keeps its ties, so with enough ties every shuffled value is \(0\). Add noise to break the ties, then test whether Falk’s median correlation is statistically different from \(0\).

Code
# xy_hat and med_cor defined above
# Many workers tie at the median years of schooling.
# So shuffling them gives a value of 0 every time.
# Years of schooling are recorded as whole numbers.
# Noise of up to half a year keeps each worker within their recorded year.
n <- nrow(xy_hat)
xy_noise_hat <- xy_hat
xy_noise_hat[, 'x_hat'] <- xy_hat[, 'x_hat'] + runif(n, -.5, .5)
xy_noise_med_cor <- med_cor(xy_noise_hat)
xy_noise_med_cor
## [1] 0.1286916

# Null Permutation Distribution of Falk's Median Correlation
null_permutation_med_cor <- rep(NA, 9999)
for(b in seq(null_permutation_med_cor) ){
    xy_perm <- xy_noise_hat
    xy_perm[, 'x_hat'] <- xy_noise_hat[sample(n, replace=FALSE), 'x_hat'] ## Reshuffle X
    xy_med_cor_b <- med_cor(xy_perm)
    null_permutation_med_cor[b] <- xy_med_cor_b
}
hist(null_permutation_med_cor, breaks=100,
    border=NA, freq=FALSE,
    xlim=range(c(null_permutation_med_cor, xy_noise_med_cor)),
    xlab="Falk's Median Correlation",
    main=NA)
title('Null Permutation Distribution', font.main=1)

## Test whether Falk's median correlation is statistically different from 0
abline(v=quantile(null_permutation_med_cor, probs=c(0.025, 0.975)), col=rgb(0, 0, 0, .8), lwd=2)
abline(v=xy_noise_med_cor, col=rgb(1, 0, 0, .8), lwd=2)

Code
# The observed value of about 0.13 is far outside the shuffled values.

Here are the same statistics and the same test for the sampled states from Hypothesis Testing.

Code
# xy2_hat and n2 defined above
xy2_codev <- codev(xy2_hat)
xy2_codev
## [1] 3.5
xy2_med_cor <- med_cor(xy2_hat)
xy2_med_cor
## [1] 0.1129032

# Null Permutation Distribution of Falk's Median Correlation
null_permutation_med_cor2 <- rep(NA, 9999)
for(b in seq(null_permutation_med_cor2) ){
    xy2_perm <- xy2_hat
    xy2_perm[, 'x_hat'] <- xy2_hat[sample(n2, replace=FALSE), 'x_hat'] ## Reshuffle X
    xy2_med_cor_b <- med_cor(xy2_perm)
    null_permutation_med_cor2[b] <- xy2_med_cor_b
}
hist(null_permutation_med_cor2, breaks=100,
    border=NA, freq=FALSE,
    xlab="Falk's Median Correlation",
    main=NA)
title('Null Permutation Distribution', font.main=1)

## Test whether Falk's median correlation is statistically different from 0
abline(v=quantile(null_permutation_med_cor2, probs=c(0.025, 0.975)), col=rgb(0, 0, 0, .8), lwd=2)
abline(v=xy2_med_cor, col=rgb(1, 0, 0, .8), lwd=2)

Code
# The observed value of about 0.11 is inside the middle 95% of shuffled values.

14.2 Factor Data

Two Ordered Factors

When the data are ordered (rankings or ordered categories), we can summarize association by counting which pairs of observations agree in direction.

ImportantKey Definition

Kendall’s rank correlation \(\hat{\text{KT}}\) (Kendall 1938) is the share of variable-pairs that agree in direction (concordant) minus the share that disagree (discordant): \[\hat{\text{KT}} = \frac{2}{n(n-1)} \sum_{i} \sum_{j > i} \text{sgn} \Bigl( (\hat{x}_{i} - \hat{x}_{j})(\hat{y}_{i} - \hat{y}_{j}) \Bigr),\] where the sign function is \[\text{sgn}(z) = \begin{cases} +1 & \text{if } z > 0\\ 0 & \text{if } z = 0 \\ -1 & \text{if } z < 0 \end{cases}.\]

Kendall’s \(\hat{\text{KT}}\) is useful for any monotone relationship: the values do not have to fall on a straight line, only have consistent ordering. Because it only uses signs, it works for any ordinal data (rankings, ordered categories) where arithmetic averages are not meaningful. The bounds \([-1, 1]\) make it directly comparable to the Pearson correlation: a value closer to \(+1\) suggests positive association in rankings, \(-1\) negative, and \(0\) no association in the ordering. For example, the table below cross-classifies two ordered variables, and \(\hat{\text{KT}}\) summarizes it.

Code
## MAIN EXAMPLE
# Load the Data
data(Wages1, package='Ecdat')
# Group hourly wages into four ordered brackets
wage_bracket <- cut(Wages1[, 'wage'], breaks=c(0, 4, 6, 8, Inf), ordered_result=TRUE)
xy_ord_hat <- data.frame(x_hat=Wages1[, 'school'], y_hat=wage_bracket) # x_hat=years of schooling, y_hat=wage bracket
table(xy_ord_hat)
##      y_hat
## x_hat (0,4] (4,6] (6,8] (8,Inf]
##    3      0     1     0       0
##    4      1     0     0       1
##    5      2     2     1       0
##    6      4     7     3       0
##    7     11     9     3       1
##    8     51    22     9       4
##    9     90    52    13       6
##    10   169   119    78      33
##    11   214   218   129     100
##    12   365   375   251     197
##    13    63    93    93      90
##    14    50    55    59     100
##    15    13    23    40      58
##    16     0     6     2       8
# Both are ordered categories: 12 years is more than 10.
# The top bracket pays more than the bottom one.

# Kendall's rank correlation uses only the order of the categories
KT <- cor(xy_ord_hat[, 'x_hat'], as.integer(xy_ord_hat[, 'y_hat']), method='kendall')
round(KT, 3)
## [1] 0.236
# The value is positive.
# So workers with more schooling tend to be in higher wage brackets.

Kendall’s \(\hat{\text{KT}}\) also applies to cardinal data, because it only uses the ordering of the values. Converting each variable to ranks first does not change the statistic.

Code
# xy_hat defined above
# Convert each variable to ranks
xy_rank_hat <- xy_hat
xy_rank_hat[, 'x_hat'] <- rank(xy_hat[, 'x_hat'])
xy_rank_hat[, 'y_hat'] <- rank(xy_hat[, 'y_hat'])
# plot(xy_rank_hat, pch=16, col=grey(0, .1))
KT_rank <- cor(xy_rank_hat[, 'x_hat'], xy_rank_hat[, 'y_hat'], method = 'kendall')
round(KT_rank, 3)
## [1] 0.219
# Same value using the raw data
round(cor(xy_hat[, 'x_hat'], xy_hat[, 'y_hat'], method='kendall'), 3)
## [1] 0.219

Compute Kendall’s \(\hat{\text{KT}}\) for the dataset \(\{(1,1),(2,3),(3,2),(4,4)\}\).

With \(n=4\) there are \(n(n-1)/2=6\) pairs. For each pair we take the sign of \((\hat{x}_{i}-\hat{x}_{j})(\hat{y}_{i}-\hat{y}_{j})\): a pair is concordant (\(+1\)) if both variables move the same way, discordant (\(-1\)) if they move oppositely.

Concordant and Discordant Pairs
Pair \(\hat{x}\) moves \(\hat{y}\) moves sign
(1,1),(2,3) up up \(+1\)
(1,1),(3,2) up up \(+1\)
(1,1),(4,4) up up \(+1\)
(2,3),(3,2) up down \(-1\)
(2,3),(4,4) up up \(+1\)
(3,2),(4,4) up up \(+1\)

There are \(5\) concordant and \(1\) discordant pair, so \[ \hat{\text{KT}} = \frac{2}{n(n-1)}\sum_{i}\sum_{j>i}\text{sgn}(\cdots) = \frac{2}{4\cdot3}(5-1) = \frac{8}{12} = \frac{2}{3} \approx 0.67. \]

Code
cor(c(1,2,3,4), c(1,3,2,4), method='kendall')
## [1] 0.6666667

You construct sampling distributions and conduct hypothesis tests for Kendall’s rank correlation statistic in the same way you do for Pearson’s Correlation statistic and Falk’s Codeviance statistic.

Test whether Kendall’s correlation statistic is statistically different from \(0\). Expand on the example below to use bootstrapping.

Code
# xy_ord_hat defined above
KT <- cor(xy_ord_hat[, 'x_hat'], as.integer(xy_ord_hat[, 'y_hat']), method='kendall')

Kendall’s rank correlation coefficient can also be used for non-linear relationships, where Pearson’s correlation coefficient often falls short. It almost always helps to visualize your data first before summarizing it with a statistic.

Two Unordered Factors

When neither variable has a meaningful order, we organize the data as a contingency table and summarize the strength of association with a single bounded score. Cramer’s V \(\hat{\text{CramV}}\) (Cramér 1946) rescales the chi-squared statistic \(\hat{\chi}^{2}\) into a \([0, 1]\) score for the strength of association between two categorical variables in a \(K\times J\) contingency table: \[\hat{\chi}^2 = \sum_{k=1}^{K} \sum_{j=1}^{J} \frac{(\hat{o}_{kj} - \hat{e}_{kj})^2}{\hat{e}_{kj}}, \qquad \hat{\text{CramV}} = \sqrt{\frac{\hat{\chi}^2 / n}{\min(J - 1, \, K - 1)}},\] where \(\hat{o}_{kj}\) is the observed frequency in cell \((k, j)\), \(\hat{e}_{kj} = \hat{\text{RF}}_{k} \cdot \hat{\text{CF}}_{j} / n\) is the expected frequency under independence, and \(\hat{\text{RF}}_{k}=\sum_{j} \hat{o}_{kj}\), \(\hat{\text{CF}}_{j}=\sum_{k} \hat{o}_{kj}\) are the row and column totals.

Cramer’s V is useful for two unordered categorical variables, such as occupation and region. Reordering the rows or columns of the table does not change it. The chi-squared piece measures how far observed counts are from what we would expect if \(\hat{x}\) and \(\hat{y}\) were independent. Dividing by \(n \cdot \min(J - 1, K - 1)\) rescales it into a bounded score, so different tables are comparable. A value of \(0\) suggests no association, and a value closer to \(1\) suggests a strong association.

For example, the table below cross-classifies two categorical variables, and the function below computes their Cramer’s V.

Code
## MAIN EXAMPLE
# Pair each worker's sex with a schooling group, as in Comparing Multiple Groups
# The schooling groups have an order, but Cramer's V ignores it.
# So we treat both variables as unordered.
school_g <- cut(Wages1[, 'school'], breaks=c(0, 11, 12, 16), labels=c('<12', '12', '>12'))
xy_cat_hat <- data.frame(x_hat=Wages1[, 'sex'], y_hat=school_g) # x_hat=sex, y_hat=schooling group
table(xy_cat_hat)
##         y_hat
## x_hat    <12  12 >12
##   female 567 594 408
##   male   786 594 345
# Men are more likely than women to have fewer than 12 years of schooling.

CV <- function(xy){
    # Create a contingency table from the categorical variables
    tbl <- table(xy)
    # Compute the chi-square statistic (without Yates' continuity correction)
    chi2 <- chisq.test(tbl, correct=FALSE)[['statistic']]
    # Total sample size
    n <- sum(tbl)
    # Compute the minimum degrees of freedom (min(rows-1, columns-1))
    df_min <- min(nrow(tbl) - 1, ncol(tbl) - 1)
    # Calculate Cramer's V
    V <- sqrt((chi2 / n) / df_min)
    return(V)
}
CV(xy_cat_hat)
## X-squared 
## 0.1007045
# Cramer's V of about 0.10 suggests a weak association.

# DescTools::CramerV( table(xy_cat_hat) )

Here is the same calculation for a second table, UCBAdmissions.

Code
# CV defined above
# Built-in: 4526 applications to six graduate departments at Berkeley in 1973
# One row per (Admit, Gender, Dept) cell, with the number of applicants in Freq
ucb_cells <- as.data.frame(UCBAdmissions)
# Repeat each row Freq times, so there is one row per applicant
ucb_hat <- ucb_cells[rep(seq_len(nrow(ucb_cells)), ucb_cells[, 'Freq']), ]
# Department and admission outcome are both unordered.
xy_cat2_hat <- data.frame(x_hat=ucb_hat[, 'Dept'], y_hat=ucb_hat[, 'Admit']) # x_hat=department, y_hat=admission outcome
table(xy_cat2_hat)
##      y_hat
## x_hat Admitted Rejected
##     A      601      332
##     B      370      215
##     C      322      596
##     D      269      523
##     E      147      437
##     F       46      668
# Department A admits about 64% of its applicants and department F about 6%.
CV(xy_cat2_hat)
## X-squared 
## 0.4148446
# Cramer's V is about 0.41.
# This reflects a strong association between department and admission.

Repeat the calculation with the census region of each sampled state from Hypothesis Testing and its murder arrest rate grouped into three bands.

Code
# state_id and CV defined above
# Regions have no natural order, so we treat both variables as unordered.
# state.region lists states in the same alphabetical order as USArrests.
# So state_id picks out the same states from both.
# Pair each sampled state's region with its murder band (low, middle, high)
murder_band <- cut(USArrests[state_id, 'Murder'], 3, labels=c('low', 'middle', 'high'))
xy_cat3_hat <- data.frame(x_hat=state.region[state_id], y_hat=murder_band) # x_hat=census region, y_hat=murder band
table(xy_cat3_hat)
##                y_hat
## x_hat           low middle high
##   Northeast       6      1    0
##   South           0      3    3
##   North Central   6      4    1
##   West            5      4    2
CV(xy_cat3_hat)
## X-squared 
## 0.4085406

You can also apply Cramer’s V to cardinal data after cutting each variable into bins. This throws away the ordering of the bins, so Kendall’s \(\hat{\text{KT}}\) or Pearson’s \(\hat{r}_{XY}\) usually make better use of the data.

Code
# xy_hat and CV defined above
# Cut each variable into bins: years of schooling into 3, hourly wages into 4
# Wages have a long right tail, so use the brackets above, not equal widths
xy_bin_hat <- xy_hat
xy_bin_hat[, 'x_hat'] <- cut(xy_hat[, 'x_hat'], 3)
xy_bin_hat[, 'y_hat'] <- cut(xy_hat[, 'y_hat'], breaks=c(0, 4, 6, 8, Inf))
table(xy_bin_hat)
##              y_hat
## x_hat         (0,4] (4,6] (6,8] (8,Inf]
##   (2.99,7.33]    18    19     7       2
##   (7.33,11.7]   524   411   229     143
##   (11.7,16]     491   552   445     453
CV(xy_bin_hat)
## X-squared 
## 0.1471187

Compute Cramer’s V for this \(2\times2\) table of \(100\) workers, cross-classifying a college degree against employment.

\[\begin{array}{c|cc|c} & \text{employed} & \text{unemployed} & \text{Row total}\\ \hline \text{degree} & 30 & 10 & 40\\ \text{no degree} & 20 & 40 & 60\\ \hline \text{Column total} & 50 & 50 & 100 \end{array}\]

If degree and employment were independent, the expected count in each cell is \(\hat{e}_{kj}=\hat{\text{RF}}_{k}\hat{\text{CF}}_{j}/n\). For the top-left cell, \(\hat{e}_{11}=40\cdot50/100=20\); all four expected counts are \((20,20,30,30)\). \[\begin{aligned} \hat{\chi}^2 &= \frac{(30-20)^2}{20}+\frac{(10-20)^2}{20}+\frac{(20-30)^2}{30}+\frac{(40-30)^2}{30} = 5+5+\frac{10}{3}+\frac{10}{3} = \frac{50}{3} \approx 16.67. \end{aligned}\] With \(K=2\) rows and \(J=2\) columns, \(\min(J-1,K-1)=1\), so \[ \hat{\text{CramV}} = \sqrt{\frac{\hat{\chi}^2/n}{\min(J-1,K-1)}} = \sqrt{\frac{(50/3)/100}{1}} = \sqrt{\frac{1}{6}} \approx 0.41. \]

Code
# Observed 2x2 table
O <- rbind(degree    = c(employed=30, unemployed=10),
           no_degree = c(20, 40))
n0 <- sum(O)

# Expected counts under independence
E <- outer(rowSums(O), colSums(O)) / n0

# Chi-square statistic and Cramer's V
chi2 <- sum((O - E)^2 / E)
df_min <- min(nrow(O)-1, ncol(O)-1)
c(chi2=chi2, CramerV=sqrt((chi2/n0)/df_min))
##       chi2    CramerV 
## 16.6666667  0.4082483

You construct sampling distributions and conduct hypothesis tests for Cramer’s V statistic in the same way you do for the other statistics.

14.3 Association is not Causation

None of the tools in this chapter alone identifies what would happen to \(Y_{i}\) if we intervened and changed \(X_{i}\) (Holland 1986). Sampling variability is one problem: relationships in the population may not appear in a sample, and apparent sample relationships may not exist in the population. Hypothesis tests help quantify that uncertainty, but statistical significance does not distinguish association from causation. Two other problems remain: real relationships can average out, and random data can acquire mechanically induced relationships. To make both problems concrete, we focus on the Pearson correlation from Cardinal Data and examine

  • Causation without correlation
  • Correlation without causation

Causation without correlation

Examples of the first problem include nonlinear effects and heterogeneous effects that average out. For example, wages tend to rise with age early in a career and then flatten or fall, and a minimum wage increase can raise pay for some workers while others lose their jobs.

Code
## MAIN EXAMPLE
n <- 10000

# X causes Y via Y = X^2 + noise
x_sim <- runif(n, min = -1, max = 1)
epsilon <- rnorm(n, mean = 0, sd = 0.1)
y_sim <- x_sim^2 + epsilon  # clear causal effect of X on Y
plot(x_sim, y_sim, pch=16, col=grey(0, .05),
    main=NA, xlab='X', ylab='Y')

# Correlation over the entire range
cor_sim <- round(cor(x_sim, y_sim), 2)
title(paste0('Cor: ', cor_sim), font.main=1)

Code
# Heterogeneous Effects
x_sim <- rnorm(n) # Randomized 'treatment'

# Heterogeneous effects based on group
group <- rbinom(n, size = 1, prob = 0.5)
epsilon <- rnorm(n, mean = 0, sd = 1)
y_sim <- ifelse(group == 1,
            x_sim + epsilon,   # positive effect
            -x_sim + epsilon)  # negative effect
plot(x_sim, y_sim, pch=16, col=grey(0, .05),
    main=NA, xlab='X', ylab='Y')

# Correlation in the pooled sample
cor_sim <- round(cor(x_sim, y_sim), 2)
title(paste0('Cor: ', cor_sim), font.main=1)

Correlation without causation

Examples of the second problem include shared denominators and selection bias that induce correlations.

Consider three completely random variables. We can induce a mechanical relationship between the first two variables by dividing them both by the third variable. Per-capita statistics are a common real case, such as crimes per resident and police officers per resident across cities, which share population as the denominator.

Code
## MAIN EXAMPLE
n <- 20000

# Independent components
A <- runif(n)
B <- runif(n)
C <- runif(n)
par(mfrow=c(1, 2))
plot(A, B, pch=16, col=grey(0, .05),
    main=NA, xlab='A', ylab='B')
cor_ab <- round(cor(A, B), 2)
title(paste0('Independent Variables\nCor: ', cor_ab), font.main=1)

# Ratios with a shared denominator
x_sim <- A / C
y_sim <- B / C
plot(x_sim, y_sim, pch=16, col=grey(0, .05),
    main=NA, xlab='X', ylab='Y')
cor_sim <- round(cor(x_sim, y_sim), 2)
title(paste0('With Common Divisor\nCor: ', cor_sim), font.main=1)

Code

# Correlation
cor(x_sim, y_sim)
## [1] 0.8094966

Consider an admissions rule into university: applicants are accepted if they have either high test scores or strong extracurriculars. Even if there is no general relationship between test scores and extracurriculars, you will see one amongst university students. Real datasets are often selected in a similar way, so their relationships describe the selected group rather than everyone.

Code

# Independent traits in the population
test_score        <- rnorm(n, mean = 0, sd = 1)
extracurriculars  <- rnorm(n, mean = 0, sd = 1)

# Selection above thresholds
threshold <- 1.0
admitted <- (test_score > threshold) | (extracurriculars > threshold)
mean(admitted)  # admission rate
## [1] 0.294

par(mfrow = c(1, 2))
# Full population
plot(test_score, extracurriculars,
     pch=16, col=grey(0, .05),
     main=NA, xlab='Test Score', ylab='Extracurriculars')
cor_all <- round(cor(test_score, extracurriculars), 2)
title(paste0('General Sample\nCor: ', cor_all), font.main=1)
# Admitted only
plot(test_score[admitted], extracurriculars[admitted],
     pch=16, col=grey(0, .05),
     main=NA, xlab='Test Score', ylab='Extracurriculars')
cor_adm <- round(cor(test_score[admitted], extracurriculars[admitted]), 2)
title(paste0('University Sample\nCor: ', cor_adm), font.main=1)

Code
# Correlation among admitted applicants only
cor(test_score[admitted], extracurriculars[admitted])
## [1] -0.5570494

Real data often mix several of these issues. Two variables can also move together because both respond to common factors. The scatterplot below shows one correlation that need not be causal.

Code
# xy_hat defined above: years of schooling and hourly wages of 3294 workers
plot(y_hat ~ x_hat, xy_hat, pch=16, col=grey(0, .1),
    main=NA, xlab='Years of schooling', ylab='Hourly wage (1980 dollars)')
cor_wages <- round(cor(xy_hat[, 'x_hat'], xy_hat[, 'y_hat']), 2)
title(paste0('Cor: ', cor_wages), font.main=1)

Code

# Correlation across workers
cor(xy_hat[, 'x_hat'], xy_hat[, 'y_hat'])
## [1] 0.2824925
# The correlation is about 0.28.
# This alone does not show that schooling raises wages.
# Both likely respond to common factors, such as ability or family background.
# Wages1 also only includes people with observed wages.
# So its relationships describe workers, not everyone who could have worked.

Repeat this for murder and assault arrests across US states in USArrests.

Code
# Arrests per 100,000 residents in each US state in 1973
plot(Murder ~ Assault, USArrests, pch=16, col=grey(0, .5),
    main=NA, xlab='Assault Arrests', ylab='Murder Arrests')
cor_arrests <- round(cor(USArrests$Murder, USArrests$Assault), 2)
title(paste0('Cor: ', cor_arrests), font.main=1)

Code
# States with more assault arrests also tend to have more murder arrests.

# Correlation across states
cor(USArrests$Murder, USArrests$Assault)
## [1] 0.8018733
# The correlation is about 0.8.
# This does not mean that assaults cause murders.
# Both variables are divided by the same state population.
# Both also likely respond to common factors.
# Examples are policing practices or local economic conditions.

Note that the examples above are not the only examples of “correlation does not mean causation” (Freedman 1991). Many real datasets have temporal and spatial interdependence that create additional issues. Many real datasets also have economic interdependence, which also creates additional issues (Adamson and Stephenson 2026).

14.4 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. For the dataset \(\{(1, 5),\ (2, 3),\ (3, 4),\ (4, 2),\ (5, 1)\}\), compute the Pearson correlation \(\hat{r}_{XY}\) by hand. Show all five steps: means, deviances, covariance \(\hat{c}_{XY}\), standard deviations \(\hat{s}_{X}\) and \(\hat{s}_{Y}\), and the final correlation.

  3. Using the \(35\) sampled states in state_id, compute both the Pearson correlation and Kendall’s rank correlation between Murder and Assault in USArrests. Then write a bootstrap with \(B = 9999\) iterations to construct a \(95\%\) confidence interval for the Pearson correlation. Does the interval contain zero?

  4. Give an example of a dataset where Pearson’s correlation \(\hat{r}_{XY}\) is close to zero but there is a clear causal relationship between \(X_{i}\) and \(Y_{i}\). What feature of the relationship causes the correlation to miss it, and which alternative statistic from this chapter might detect it?

Further Reading

Recall

The right statistic of association depends on the type of data. For cardinal data, we computed the Pearson correlation by hand for \(\{(0, 0.1), (1, 0.3), (2, 0.2)\}\) and found \(\hat{r}_{XY}=1/2\). For cardinal data with outliers, Falk’s median correlation \(\tilde{r}_{XY}\), built from the codeviance \(\tilde{c}_{XY}\), is a robust alternative. For ordered factors, the six pairs of \(\{(1, 1), (2, 3), (3, 2), (4, 4)\}\) gave Kendall’s \(\hat{\text{KT}}=2/3\). For unordered factors, the \(2\times 2\) table of degree and employment gave Cramer’s V \(\hat{\text{CramV}}=\sqrt{1/6}\approx 0.41\). The correlation between schooling and wages in Wages1 was about \(0.28\), with a \(95\%\) bootstrap confidence interval from about \(0.25\) to \(0.32\). In the permutation test, the shuffled correlations stayed within a few hundredths of \(0\), so \(0.28\) is far outside them. These statistics describe association rather than causation. Simulated examples showed causation without correlation, when effects average out, and correlation without causation, from shared denominators or selection.

Adamson, Jordan, and Daniel Stephenson. 2026. “The Endogenous Nature of War and Its Economic Consequences.” European Economic Review 189: 105430. https://doi.org/10.1016/j.euroecorev.2026.105430.
Cramér, Harald. 1946. Mathematical Methods of Statistics. Princeton University Press. https://doi.org/10.1515/9781400883868.
Falk, Michael. 1997. “On MAD and Comedians.” Annals of the Institute of Statistical Mathematics 49 (4): 615–44. https://doi.org/10.1023/A:1003258024248.
Freedman, David A. 1991. “Statistical Models and Shoe Leather.” Sociological Methodology 21: 291–313. https://doi.org/10.2307/270939.
Holland, Paul W. 1986. “Statistics and Causal Inference.” Journal of the American Statistical Association 81 (396): 945–60. https://doi.org/10.1080/01621459.1986.10478354.
Kendall, M. G. 1938. “A New Measure of Rank Correlation.” Biometrika 30 (1-2): 81–93. https://doi.org/10.1093/biomet/30.1-2.81.
Pearson, Karl. 1896. “Mathematical Contributions to the Theory of Evolution. III. Regression, Heredity, and Panmixia.” Philosophical Transactions of the Royal Society of London. Series A 187: 253–318. https://doi.org/10.1098/rsta.1896.0007.
Sen, Pranab Kumar. 1968. “Estimates of the Regression Coefficient Based on Kendall’s Tau.” Journal of the American Statistical Association 63 (324): 1379–89. https://doi.org/10.1080/01621459.1968.10480934.
Theil, Henri. 1950. “A Rank-Invariant Method of Linear and Polynomial Regression Analysis, I, II, III.” Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen 53: 386–92, 521–25, 1397–412. https://doi.org/10.1007/978-94-011-2546-8_20.

  1. See Center/Spread based on Quantiles. See also the Theil-Sen Estimator (Theil 1950; Sen 1968), which may be seen as a precursor.↩︎