13  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; Association and Prediction Are Not Causation explains why they do not identify causal effects.

13.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 the building block for every other linear two-variable statistic: it is positive when \(\hat{X}_{i}\) and \(\hat{Y}_{i}\) tend to be above or below their means together, 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}\) 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 workhorse 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. A value close to \(-1\) suggests negative association, a value close to \(0\) suggests no linear association, and a value close to \(+1\) suggests positive association.

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
X <- c(0, 1, 2)
X
## [1] 0 1 2
Y <- c(0.1, 0.3, 0.2)
Y
## [1] 0.1 0.3 0.2

# Compute the Means
mX <- mean(X)
mY <- mean(Y)

# Compute the Deviances
dev_X <- X - mX
dev_Y <- Y - mY

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

# Compute the Standard Deviations
var_X <- sum(dev_X^2) / length(X)
sd_X <- sqrt(var_X)
var_Y <- sum(dev_Y^2) / length(Y)
sd_Y <- sqrt(var_Y)

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

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

Hypothesis Testing

You can conduct hypothesis tests for these statistics using the same procedures we learned for univariate data, ch.5. For example, by inverting a confidence interval.

Code
# Load the Data
xy <- USArrests[, c('Murder', 'UrbanPop')]
xy_cor <- cor(xy[, 1], xy[, 2])
#plot(xy, pch=16, col=grey(0, .25))
    
# Bootstrap Distribution of Correlation
n <- nrow(xy)
bootstrap_cor <- rep(NA, 9999)
for(b in seq(bootstrap_cor) ){
    xy_b <- xy[sample(n, replace=TRUE), ]
    xy_cor_b <- cor(xy_b[, 1], xy_b[, 2])
    bootstrap_cor[b] <- xy_cor_b
}
hist(bootstrap_cor, 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_ci <- quantile(bootstrap_cor, probs=c(0.025, 0.975))
abline(v=boot_ci)
abline(v=0, col=rgb(1, 0, 0, .8))

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

Code
xy <- USArrests[, c('Murder', 'UrbanPop')]
xy_cor <- cor(xy[, 1], xy[, 2])
#plot(xy, pch=16, col=grey(0, .25))
    
# Null Bootstrap Distribution of Correlation
n <- nrow(xy)
null_bootstrap_cor <- rep(NA, 9999)
for(b in seq(null_bootstrap_cor) ){
    xy_b <- xy
    xy_b[, 'UrbanPop'] <- xy[sample(n, replace=TRUE), 'UrbanPop'] ## Reshuffle X
    xy_cor_b <- cor(xy_b[, 1], xy_b[, 2])
    null_bootstrap_cor[b] <- xy_cor_b
}
hist(null_bootstrap_cor, breaks=100,
    border=NA, freq=FALSE,
    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)
abline(v=xy_cor, col=rgb(0, 0, 1, .8))

Because all dependence resides in the pairing, and permuting one margin destroys the pairing completely.

To construct a permutation test, we need to use replace=FALSE. Rework the above code to make a Permutation Null Distribution and conduct a permutation test.

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\) \(X,Y\): Deterministically leave-one-out observation Variance estimate (after rescaling) and Normal CI estimate
Bootstrap \(n\) \(B\) \(X,Y\): Random resample with replacement Percentile CI estimate
Null Bootstrap \(n\) \(B\) \(X,Y\): Random resample with replacement and shifted Percentile CI under imposed null, \(p\)-values
Permutation \(n\) \(B\) \(X\): Random resample without replacement Percentile CI under imposed null of no association, \(p\)-values

Falk Codeviance

When the data contain outliers, the mean-based covariance can be misleading and a robust alternative is needed.

The Falk codeviance \(\tilde{C}_{XY}\) 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\}.\]

The 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 (income, prices, response times), where one extreme observation can move \(\hat{C}_{XY}\) by a lot but \(\tilde{C}_{XY}\) by very little (the same robustness rationale as for the median, \(IQR\), and \(\text{MAD}\) from Part 1).1 Unlike the Pearson correlation, the median correlation \(\tilde{R}_{XY}\) typically lies in \([-1,1]\) but not always.

Code
codev <- function(xy) {
  # Compute medians for each column
  med <- apply(xy, 2, median)
  # Subtract the medians from each column
  xm <- sweep(xy, 2, med, '-')
  # Compute CoDev
  CoDev <- median(xm[, 1] * xm[, 2])
  # Compute the medians of absolute deviation
  MadProd <- prod( apply(abs(xm), 2, median) )
  # Return the robust correlation measure
  return( CoDev / MadProd)
}
xy_codev <- codev(xy)
xy_codev
## [1] 0.005707763

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 the median correlation is \(\tilde{R}_{XY}=2/(1\cdot1)=2\). Here \(\tilde{R}_{XY}\) exceeds \(1\): unlike Pearson’s correlation, the median correlation is not confined to \([-1,1]\).

Code
xy_small <- cbind(X=c(1,2,3,4,5), Y=c(2,1,4,3,6))
codev(xy_small)
## [1] 2

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

Code
xy <- USArrests[, c('Murder', 'UrbanPop')]
xy_cor <- cor(xy[, 1], xy[, 2])
#plot(xy, pch=16, col=grey(0, .25))
    
# Null Permutation Distribution of Codeviance
n <- nrow(xy)
null_permutation_codev <- rep(NA, 9999)
for(b in seq(null_permutation_codev) ){
    xy_b <- xy
    xy_b[, 'UrbanPop'] <- xy[sample(n, replace=FALSE), 'UrbanPop'] ## Reshuffle X
    xy_codev_b <- codev(xy_b)
    null_permutation_codev[b] <- xy_codev_b
}
hist(null_permutation_codev, breaks=100,
    border=NA, freq=FALSE,
    xlab='Codeviance',
    main=NA)
title('Null Permutation Distribution', font.main=1)

## Test whether correlation is statistically different from 0
abline(v=quantile(null_permutation_codev, probs=c(0.025, 0.975)))
abline(v=xy_codev, col=rgb(0, 0, 1, .8))

13.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{KT}\) is the share of variable-pairs that agree in direction (concordant) minus the share that disagree (discordant): \[\hat{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{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.

Code
xy <- USArrests[, c('Murder', 'UrbanPop')]
xy[, 1] <- rank(xy[, 1] )
xy[, 2] <- rank(xy[, 2] )
# plot(xy, pch=16, col=grey(0, .25))
KT <- cor(xy[, 1], xy[, 2], method = 'kendall')
round(KT, 3)
## [1] 0.074

Compute Kendall’s \(\hat{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.

Pair \(X\) moves \(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{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 as for Pearson’s Correlation statistic and Falk’s Codeviance statistic.

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

Code
xy <- USArrests[, c('Murder', 'UrbanPop')]
KT <- cor(xy[, 1], xy[, 2], 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 visual 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.

ImportantKey Definition

Cramer’s V \(\hat{CV}\) 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{CV} = \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{RF}_{k} \cdot \hat{CF}_{j} / n\) is the expected frequency under independence, and \(\hat{RF}_{k}=\sum_{j} \hat{O}_{kj}\), \(\hat{CF}_{j}=\sum_{k} \hat{O}_{kj}\) are the row and column totals.

Cramer’s V is useful for two unordered categorical variables (occupation, region, brand), because it treats every reordering of rows or columns the same and makes no use of any ordering. The chi-squared piece measures how far observed counts are from what we would expect if \(\hat{X}\) and \(\hat{Y}\) were independent, and dividing by \(n \cdot \min(J - 1, K - 1)\) rescales it into a bounded score so different tables are comparable: \(0\) suggests no association, and a value closer to \(1\) suggests a strong association.

Code
xy <- USArrests[, c('Murder', 'UrbanPop')]
xy[, 1] <- cut(xy[, 1], 3)
xy[, 2] <- cut(xy[, 2], 4)
table(xy)
##               UrbanPop
## Murder         (31.9,46.8] (46.8,61.5] (61.5,76.2] (76.2,91.1]
##   (0.783,6.33]           4           5           8           5
##   (6.33,11.9]            0           4           7           6
##   (11.9,17.4]            2           4           2           3

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)
## X-squared 
## 0.2307071

# DescTools::CramerV( table(xy) )

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{RF}_{k}\hat{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+3.33+3.33 = 16.67. \end{aligned}\] With \(K=2\) rows and \(J=2\) columns, \(\min(J-1,K-1)=1\), so \[ \hat{CV} = \sqrt{\frac{\hat{\chi}^2/n}{\min(J-1,K-1)}} = \sqrt{\frac{16.67/100}{1}} \approx 0.41. \]

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

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

# 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/n)/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 as the other statistics.

13.3 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 USArrests dataset, compute both the Pearson correlation and Kendall’s rank correlation between Murder and Assault. Then write a bootstrap with \(B = 9999\) iterations to construct a \(95\%\) confidence interval for the Pearson correlation. Does the interval contain zero?

Further Reading

Recall

This chapter summarized two-variable association with a small toolbox matched to data type: Pearson correlation \(\hat{R}_{XY}\) for cardinal data (worked out by hand on \(\{(0, 0.1), (1, 0.3), (2, 0.2)\}\), yielding \(\hat{R}_{XY}=1/2\)), Falk codeviance \(\tilde{C}_{XY}\) as the robust cousin for outlier-prone data, Kendall’s \(\hat{KT}\) for ordered factors (walked through the six pairs of \(\{(1, 1), (2, 3), (3, 2), (4, 4)\}\) to get \(\hat{KT}=2/3\)), and Cramer’s V \(\hat{CV}\) for unordered factors (the \(2\times 2\) degree-vs-employment table giving \(\hat{CV}\approx 0.41\)). These statistics describe association rather than causation, a distinction developed in Association and Prediction Are Not Causation. The next chapter takes the cardinal case and fits a line, which moves from “how strongly do \(X\) and \(Y\) move together” to “by how much does \(Y\) change per unit of \(X\)”.


  1. See Medians and Absolute Deviations. See also the Theil-Sen Estimator, which may be seen as a precursor.↩︎