We want to emphasize that most of our work relies heavily on covariance structures. Everything from correlation to regression to ANOVA to structural equation models is, at bottom, an argument about how variables vary together. So before we build models, we need to understand the raw material they are built from. Variance measured how a single variable spreads around its own mean; covariance measures how two variables spread around their means at the same time.
Science usually begins with a hunch about a cause: sleep affects mood, study time affects grades, a drug affects recovery. We cannot see causes directly; we can only watch variables move. If X truly drives Y, then when X is high we expect Y to be high (or low) systematically - the two should co-vary. Covariance is our first, humblest measurement of that co-movement. It cannot prove the cause (that famous warning is coming), but a cause that leaves no covariance behind is a cause with no evidence.
Recall variance: the average product of a deviation with itself, \(\frac{\sum(x-\bar x)(x-\bar x)}{n-1}\). Covariance makes one small change - it multiplies each variable’s deviation by the other’s:
\[\text{cov}(x, y) = \frac{\sum (x - \bar{x})(y - \bar{y})}{n - 1}\]
Think about the sign of that product. When a point is above the mean on both variables (or below on both), the two deviations share a sign and their product is positive. When it is high on one and low on the other, the product is negative. Sum them up: mostly-positive products mean the variables rise together; mostly-negative means one rises as the other falls.
set.seed(5)
study <- rnorm(100, 5, 2)
grade <- 60 + 4 * study + rnorm(100, 0, 5) # grade rises with study
# by hand:
cov_hand <- sum((study - mean(study)) * (grade - mean(grade))) / (length(study) - 1)
c(by_hand = cov_hand, r_cov = cov(study, grade)) # match by_hand r_cov
15.48643 15.48643 CORRELATIONS
/VARIABLES=study grade
/STATISTICS=XPROD.using Statistics
cov(df.study, df.grade)Positive, as expected: more study, more grade. Let’s see it, colouring each point by the sign of its deviation product:
sign_prod <- (study - mean(study)) * (grade - mean(grade))
plot(study, grade, pch = 19, col = ifelse(sign_prod > 0, "steelblue", "red"),
main = "Blue points push covariance up; red points push it down")
abline(v = mean(study), h = mean(grade), lty = 2, col = "grey50")
The blue points (agreeing deviations) dominate, so covariance is positive. That picture is covariance.
Covariance has a serious drawback for interpretation: its size depends on the units. Measure grade in points versus proportions and the covariance changes a lot, even though the relationship is identical. This is the same units problem the standard deviation had, and we fix it the same way, by standardizing. Divide covariance by both standard deviations and you get the correlation coefficient, \(r\):
\[r = \frac{\text{cov}(x, y)}{s_x \, s_y}\]
r_hand <- cov(study, grade) / (sd(study) * sd(grade))
c(by_hand = r_hand, r_cor = cor(study, grade)) # match by_hand r_cor
0.844582 0.844582 CORRELATIONS
/VARIABLES=study grade.using Statistics
cor(df.study, df.grade)Correlation is covariance stripped of its units, squeezed onto a fixed scale from \(-1\) (perfect negative) through \(0\) (no linear relationship) to \(+1\) (perfect positive). Now the number means the same thing whatever the variables - which is why \(r\), not raw covariance, is what we report.
Two lasting warnings. First, the familiar one: correlation is not causation. Ice-cream sales and drownings correlate (both driven by summer heat), but ice cream does not cause drowning. Second, and less appreciated: \(r\) measures only the linear part of a relationship. A strong curved relationship can produce \(r \approx 0\). This is why we told you to always plot your data - the single number can hide the whole story:
# a perfect U-shape: strong relationship, ~zero correlation
xx <- seq(-3, 3, length.out = 100)
yy <- xx^2
c(correlation = round(cor(xx, yy), 3)) # ~0, yet Y is completely determined by X!correlation
0 plot(xx, yy, pch = 19, col = "darkorange",
main = "r is near zero, but the relationship is perfect (just not linear)")
The U-shape above was one warning. Here is the most famous one, and the reason “always plot your data” belongs on the first page of every statistics course. In 1973, Francis Anscombe built four small datasets that share almost every summary statistic you would normally report - the same mean of \(x\), the same mean of \(y\), the same variances, the same correlation, and the same regression line - yet describe four completely different worlds. R ships them in the built-in anscombe data:
a <- anscombe
sapply(1:4, function(i) round(cor(a[[paste0("x", i)]], a[[paste0("y", i)]]), 3))[1] 0.816 0.816 0.816 0.817Four identical correlations, \(r = 0.816\). Every set also has the same fitted line, \(\hat{y} = 3 + 0.5\,x\). If you reported only the numbers, you would call these datasets interchangeable. Now plot them:
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
for (i in 1:4) {
x <- a[[paste0("x", i)]]; y <- a[[paste0("y", i)]]
plot(x, y, pch = 19, col = "steelblue", xlim = c(3, 20), ylim = c(2, 14),
main = paste("Dataset", i), xlab = "x", ylab = "y")
abline(lm(y ~ x), col = "firebrick", lwd = 2) # the SAME line every time
}
par(op)The same red line runs through all four, but only the first tells the story the line claims:
Identical statistics, four different truths. The picture told you in a moment what the correlation coefficient never could. This is why the graphical display of data is not decoration; it is often the fastest way to catch a relationship that is curved, contaminated, or driven by a single point.
Look again at what the scatter plots did: they showed you where the straight line was wrong. That is exactly the job of a residual. Recall the chain that runs through this book. We began with dispersion, how a variable spreads, and dispersion is the raw material of uncertainty. Once we start making predictions, that uncertainty splits into two pieces: the predicted value \(\hat{y}\), what the model expects, and the residual \(y - \hat{y}\), what the model missed. The residuals are where the model is wrong, and it is the residuals, not the fit, that tell you whether the model is any good and whether your assumption of a straight-line relationship holds.
Fit a line to the curved Dataset 2 and plot the residuals, and the curve the correlation hid comes right back:
m2 <- lm(y2 ~ x2, data = anscombe)
plot(fitted(m2), resid(m2), pch = 19, col = "steelblue",
xlab = "fitted (predicted) value", ylab = "residual (y minus predicted)",
main = "Dataset 2 residuals arch — the misfit the scatter revealed")
abline(h = 0, lty = 2, col = "grey50")
The residuals bow upward at the ends and dip in the middle, the signature of a curved relationship forced through a straight line. The scatter plot and the residual plot told the same story.
Here is why that matters far beyond four small datasets. With two variables you can see a bad linear fit directly in the scatter plot. But real models carry many predictors, and you cannot plot a relationship in five or ten dimensions. Anscombe’s quartet stops being usable the moment you reach multiple regression, because there is no two-dimensional picture left to draw. What survives is the residual. Residual analysis carries the “look at where the line is wrong” habit from the simple scatter plot into models too complex to plot. So treat Anscombe’s quartet as a bridge, not a destination: it teaches, in pictures anyone can read, the very thing you will later have to check with residuals once the pictures run out.
One more idea deserves its own heading, because it runs underneath everything in this book: to covary, two variables must first vary. Covariance is built from deviations, \((x - \bar{x})(y - \bar{y})\), so if \(x\) barely moves there is nothing to multiply and no covariance to find, however strong the true relationship. A blunt but useful heuristic follows: the best predictors are usually the ones with the most variance, because a predictor pinned near a single value cannot explain anything.
The practical consequence is restriction of range: sample only a narrow slice of a variable and you shrink its variance and, with it, the correlation you can observe, even though the real relationship is unchanged. Watch a genuine correlation collapse when we keep only the middle of \(x\):
set.seed(2); n <- 1000
x <- rnorm(n)
y <- 0.6 * x + rnorm(n, sd = 0.8) # a real relationship
mid <- x > -0.5 & x < 0.5 # keep only a narrow band of x
c(full_range = round(cor(x, y), 2),
restricted_range = round(cor(x[mid], y[mid]), 2)) full_range restricted_range
0.62 0.23 Same data-generating process, same underlying relationship, but the restricted-range correlation is a fraction of the full one. This is why a study of only high-scoring students finds weak predictors of success (everyone is already high), and why a lab that tests a narrow slice of the population can miss an effect that is real in the wider world. It cuts the other way too: unreliable measures add error variance that dilutes the true signal, attenuating relationships for the same underlying reason. So a recurring question in careful analysis is simply where is the variance, and how much of it is there? Variance is the fuel every relationship in this book runs on.
Here is the payoff that carries into the rest of the book. The slope of a simple regression is literally covariance divided by variance:
b_from_cov <- cov(study, grade) / var(study)
c(cov_over_var = b_from_cov, lm_slope = coef(lm(grade ~ study))[["study"]])cov_over_var lm_slope
4.332777 4.332777 REGRESSION
/STATISTICS COEFF R ANOVA
/DEPENDENT grade
/METHOD=ENTER study.using GLM
m = lm(@formula(grade ~ study), df)
coeftable(m)They match exactly. Regression, correlation, and everything downstream are built from covariance. When we say the General Linear Model “models covariance structure,” this is what we mean, and it is why this chapter is the real beginning of the modeling half of the book.
cov() and cor(). What is the sign, and why?cov(x, y) / var(x) equals the slope of lm(y ~ x) for your data.We now hold the engine (covariance) but have not yet built the machine. Before we assemble full models, a short detour to make sure our sensors are trustworthy: are our measurements reliable and valid? A model built on a broken ruler predicts nothing, no matter how strong the covariance looks.