“Graphs are essential to good statistical analysis.” That’s the introductory comment from Graphs in Statistical Analysis by FJ Anscombe. Published in 1973 in The American Statistician, the article is still being referenced 50+ years later. The easy explanation for this longevity is that Anscombe provided a compelling case, in the form of four sets of paired data, that have come to be known as the “Anscombe Quartet.” Two x variables (X,X4) and four y variables (Y1,Y2,Y3,Y4) are combined to form four paired sets (X with each of Y1, Y2 and Y3; X4 with Y4), with the same correlation coefficients and regression equations (slope and intercept). Some of the summary statistics:

• \(N = 11\) for each variable
• \(\overline{x}\) is 9.0
  
• Variance of x is 11.0
  
• \(\overline{y}\) is 7.5
  
• Variance of y is 4.125
  
• Correlation is 0.816
  
• Regression line is \(y = 0.5x + 3.0\)
  

The plot twist, so to speak, is that the scatterplots don’t look anything alike. While the summary statistics lead us in one direction, the scatterplots turn us around to take a deeper look. The goal in this report is to get a better sense of how Anscombe developed the datasets.

When studying these figures, it’s helpful to keep a few ideas and formulas handy. For example, any simple linear regression equation goes through the value \(\overline{x},\overline{y}\). Here, each regression line goes through (9, 7.5).

Also, the slope and correlation coefficient are related. This can be expressed in a number of ways, one of which is: \[slope = r\cdot\sqrt{Var(y)/Var(x)}\] Rearranged, we have \[r = slope\cdot\sqrt{Var(x)/Var(y)}\]

For these datasets, \(r = 0.5\cdot\sqrt{11/4.125}\)

\(11/4.125 = 8/3\) so \(r = 0.5\cdot\sqrt{8/3}\)

The squareroot of 8 is 2 times the squareroot of 2 so \(r = 0.5\cdot2\cdot\sqrt{2/3}\) = \(\sqrt{2/3}\)

In shorter form, the nominal r-squared value is \(2/3\) and the nominal r value is 0.8164966.

Given the X values (the integers 4, 5, 6, …14), how would one obtain a set of Y values consistent with the results from Figure 1? That is, we want a set of Y values consistent with \(y = 0.5x + 3.0\), \(r = 0.816\), \(mean = 7.5\), \(variance = 4.125\), and so on.

Note that Figure 1 can be thought of in terms of two groups of data: 1) points that fall on or very close to the regression line, 2) points that fall distinctly above or below the regression line. My sense is that Anscombe made a Y dataset that was linearly related to X, with correlation 1.0, and then simply “moved” data points away from the line until he had the dataset he needed.

In the Semistructured tab, the plot shows data pairs that I formed after thinking in terms of “on the line” or “away from the line.” I made a set of data with a mean of 7.5 and a variance of 4.125 (uniform data), initally with a correlation of 1.0 with the X values. I then selectively swapped Y values (keeping the mean and variance intact) to approach the desired correlation and slope/intercept properties. I’m showing some R code in the next tab to give a sense of what this looks like in practice.

Freewheeling it, I formed a set of data that has the properties of the existing Anscombe datasets but with a distinct appearance.

If you want a dataset with a given mean and standard deviation it’s a simple task. For these purposes, it’s reasonable to make a uniform dataset that, initally, will correlate (r = 1.0) with the X dataset. The starting point doesn’t matter, but let’s begin with the integers 1,2,3, …, 10,11. Call these values Z. The mean of these values is 6 and the standard deviation is 3.32 (squareroot of 11, the variance). We want a dataset with mean 7.5 and a standard deviation of 2.03 (squareroot of 4.125, the variance of the Y datasets). We then do the following:

\({y}_{i} = ({z}_{i} - 6)\cdot\sqrt{Var(4.125)/Var(11)} +7.5\)

This results in the following values: 4.44, 5.05, 5.66, 6.28, 6.89, 7.5, 8.11, 8.72, 9.34, 9.95, 10.56.

We can pair each of these values with an X value (after we sort the X’s from low to high). After some swapping of y values, I modified a few of the values slightly to result in a dataset that keeps the properties of the original Y1 dataset. Sorted: 4.44, 5.05, 5.66, 6.27, 6.88, 7.51, 8.11, 8.75, 9.32, 9.95, 10.56.

While I have some R code to just to illustrate the basics, my actual approach involved a spreadsheet and covariance calculations. For those interested in doing something similar: Those covariance calculations are a little easier to follow when the X dataset is all integers and has an integer mean.

This is repreresentative of the work I did in a spreadsheet to come up with my own dataset for Y1. The spreadsheet included columns to compute the covariance, as well as the summary statistics that we need to match the original Y statistics.

Given the X values (the integers 4, 5, 6, …14), how would one obtain a set of Y values consistent with the results from Figure 2? That is, we want a set of Y values consistent with \(y = 0.5x + 3.0\), \(r = 0.816\), \(mean = 7.5\), \(variance = 4.125\), and so on, but we want values that are a function of x, such as \(y = ax^2+bx+c\).

The Basis My best guess was to start by using the regression equation to yield two data points. To form a dataset consistent with a quadratic equation and the desired regression equation, it seemed reasonable to use X4 values that weren’t the lowest or highest in the dataset. The figure in The Basis tab shows the starting values, (6,6) and (12,9). To build a quadratic equation with two starting points, I picked a constant value (the ‘c’ in \(y = ax^2+bx+c\)) and then adjusted it again and again until the overall equation made sense.

The initial results using this method yielded a dataset that has all the qualities we’re looking for, except for the mean (7.373). This is trivial since we need only shift the generated values by the difference (desired minus generated, or 7.5 - 7.373). The initial quadratic parameters: a = -0.1267, b = 2.78, c* = -6.12; adjusted to fit the mean of all the Y datasets: a = -0.1267, b = 2.78, c = -5.993

Wait, There’s More It might be inevitable, but it is interesting that there’s a distinct pattern to the Y2 dataset. For example, the lowest X value (4) is paired with 3.10, and the highest X value (14) is paired with 8.10. We see the same digits after the decimal (.10) and if we connect a line segment between the two points the slope is 0.5. This pattern is repeated as we move inward from the min and max X values. Note, too, that the slope of the quadratic is 0.5 at the mean/median of X, 9.00.

It is also interesting that ultimately it did not seem to matter which data points I used as my starting points. That is, I tried using x=4 (paired with 0.54+3=5) and x=14 (paired with 0.514+3=10) and got the same dataset as the first try. Other starting points also results in the Y4 dataset.

Just One More Thing We generated the Y2 dataset by decrementing the quadratic’s constant term (c, the y intercept) and solving a series of simultaneous equations. There is another variable that is readily available to us that fits the Anscombe criteria. Instead of determining an equation that gives us a concave function, we could solve for an equation that gives us a convex function.

I’m not showing the code, but I made a couple of minor changes such that the constant term counts up from zero instead of down. With the original Y2 variable, we get a ‘c’ term of around -6. With an alternate solution, we get a term of around 12. The initial quadratic parameters: a = 0.1267, b = -2.78, c* = 12.12; adjusted to fit the mean of all the Y datasets: a = 0.1267, b = -2.78, c = 11.993.

Here’s the dataset: Y2 (convex) = 6.86, 5.86, 10.26, 6.24, 7.74, 11.90, 5.87, 6.90, 8.87, 5.74, 6.26.

The mean is 7.5, the variance is 4.125, the slope is 0.5, the intercept is 3.0, and the correlation coefficient is 0.816. This Y2 variable has all the same properties as shown in the Wait, There’s More figure.

This is a simple approach to iterate the constant term, c, to yield two equations in two unknowns. \(y = 0.5x + 3.0\) at x=6 gives y=6; at x=12 the result is y=9. Set \(ax^2+bx+c\) = 6 for x = 6, and \(ax^2+bx+c\) = 9 for x = 12. Or, \(a6^2+b6+c\) = \(a36+b6+c\) = 6 and \(a12^2+b12+c\) = \(a144+b12+c\) = 9. Note the final results are the same as the original Y2 dataset.

In one manner, Figure 3 is the least interesting of the figures. FJA gave us the equation to fit 10 of the Y values (and even if he hadn’t, it would be easy to derive). The Y3 data values add to 82.50 (average 7.5), the only Y dataset that is exact this way. Still, we want a set of Y values consistent with \(y = 0.5x + 3.0\), \(r = 0.816\), \(mean = 7.5\), \(variance = 4.125\), and so on, but we want 10 values that are a function of x that is decidedly not \(y = 0.5x + 3.0\).

My basic approach to generate the Y3 data consists of a routine to ‘guess’ a starting slope and intercept (for the 10 linear values), and then iterate until all the necessary statistics are achieved (or within reach). Since 10 of the Y values fall on a line, the 11th value need simply be the difference between 82.50 and the sum of the 10 values. This approach is shown in The Code section (next tab).

While this approach proved fruitful, it can be augmented by restricting the slope and intercept values. I’m not including it in the R code, but see The Math tab for more information.

Note that the code will yield the slope (0.346) and intercept(4.0) that Anscombe states in the article. Also note that the Y3 variable initially generated with this method is very close to but not quite the actual dataset. The one value that doesn’t match is the lowest Y value, 5.38. Anscombe’s result is 5.39. Call this initial set of values Y3a. The sum of the Y3a values is 82.49, the variance is 4.127, and the correlation with X3 is 0.8166. The sum of the Y3a values was 82.50 (it is forced in the code) but it went down after rounding all the values. The simplest way to adjust the Y3a values is to add 0.01 to 5.38, which gives us a sum of 82.50, a variance of 4.123, and a correlation of 0.8163. It also gives us the original Y3 values.

A Similar Figure shows a different dependent variable from the one Anscombe gave us. In this case, instead of searching for a linear relationship among the 10 “well behaved” data pairs with slope less than 0.5, I modified the R code to look for a relationship with a slope greater than 0.5 (I’m not showing the code, but it’s a simple modification). The resulting slope (0.654) is greater than the overall slope (0.5) by the same amount the overall slope is greater than the original value (0.346); similarly, the intercept (2.0) is less than the overall intercept (3.0) by the same amount the overall intercept is less than the original value (4.0).

The alternate Y3 data: 8.54, 7.23, 6.26, 7.89, 9.19, 11.16, 5.92, 4.62, 9.85, 6.57, 5.27. The summary statistics are the same for this Y3 variable as all the other Y variables.

This is a simple approach to find slope (m) and intercept (b) terms for the 10 x/y pairs that form a linear relationship.

x <- c(4,5,6,7,8,9,10,11,12,13,14)  #simpler with sorted data
b <- 3.5  #pick an intercept starting value
rdelta <- 0  #correlation
vdelta <- 0  #variance
best <- 0.0025  #for correlation
slopebest <- 0.0 #track best slope fit
intbest <- 0.0 #track best intercept
vb <- 0.010  #for variance

for (j in 1:500)  {
  m <- 0.24  #slope starting point
  b <- b + 1/100
  for (k in 1:1000) {
    m <- m + 1/2500
    y <- x*m + b
    partial <- sum(y) - y[10]
    y[10] <- 82.5 - partial
    r <- cor(x, y, method = "pearson")
    rdelta <- abs(r-0.8165)
    vdelta <- abs(var(y) - 4.125)
    if ((rdelta < best) && (vdelta < vb)) {
      best <- rdelta
      slopebest <- m
      intbest <- b
    }
  } 
}
y <- x*slopebest + intbest
partial <- sum(y) - y[10]
y[10] <- 82.5 - partial

print(paste("Slope =",round(slopebest,digits=3), "and intercept =",round(intbest,digits=3)))

[1] "Slope = 0.346 and intercept = 4"

y <- round(y,2)
y

 [1]  5.38  5.73  6.08  6.42  6.77  7.11  7.46  7.81  8.15 12.74  8.84

ybar <- round(mean(y),digits=2)
yvar <- round(var(y),digits=3)
r <- cor(x, y, method = "pearson")
slint <- lm(y~x)
cf <- coef(slint)
Intercept <- round(cf[[1]],digits=2)
Slope <- round(cf[[2]],digits=2)
Y3summ <- data.frame(ybar,yvar,round(r,digits=4),round(Slope,digits=3),round(Intercept,digits=3))
Y3summ <- setNames(Y3summ,c("avg", "var", "corr", "slope", "int"))
Y3summ

  avg   var   corr slope int
1 7.5 4.127 0.8166   0.5   3

Given a linear relationship between 10 of the 11 data points, we can restrict the possibilities in the code to those that adhere to the linear relationship.

Figure 4 is based on X data distinct from the other three figures. X (figures 1, 2 and 3) is the set of integers [4,14], while X4 consists of 10 values equal to 8 and one value equal to 19. Why these values, and what do they suggest about the Y4 values?

It’s a good bet that Anscombe was doing some of the math in his head or with paper/pencil, but what follows is a methodical approach. With 10 values equal to what we’ll call \({x}_1\) and one value equal to \({x}_2\), we know that the sum of these values needs to be 99 (to yield an average of 9.0). The sum of squared differences (each value minus the mean, quantity squared) needs to be 110 (to give a sample variance of 11).

The math on the next tab shows how to obtain \({x}_1\) and \({x}_2\) values that result in the X4 dataset. It is interesting though that this rather lengthy approach gives us another combination of \({x}_1\) and \({x}_2\) consistent with the desired mean and variance: Ten values equal to 10 and one value equal to -1 result in a mean of 9.0 and a variance of 11.

To pair this new dataset with Y values that are consistent with the other Y values (mean, variance, correlation, regression results) we don’t need to work hard. At x = -1, 0.5(-1) + 3 = 2.5. The remaining ten y values – to give an overall mean of 7.5 – need to sum to 80. Since the original Y4 values at x = 8 sum to 70, we need only to add 1.0 to each of them to give us a new dataset consistent with the original statistics. The results are shown in the right-side tab Wait, There’s More.

This is a bit anticlimactic now that we have a new pairing to consider, but given the original X4 values, how difficult would it be for us to generate Y4 values that meet all the necessary criteria? If we start with an x value of 19 and a y value of 12.5 as a given, consider a set of 10 y values that fit the basic narrative of \(\bar{y}=7.5\). That is, without initial concern for the overall variance or other statistics, what are 10 y values that average 7.0 (so that the overall average of the 11 Y4 values is 7.5)? Here is my initial set: 5.2, 5.6, 6.0, 6.4, 6.8, 7.2, 7.6, 8.0, 8.4, 8.8. These are equally-spaced values that average 7.0. They are similar in range to the Y3 variable (ignoring Y3’s value at 12.74, which itself is similar to Y4’s value at 12.5). If we calculate the other key statistics, we see the following:

• \(\overline{y}\) is 7.5
  
• Variance of y is 4.07
  
• Correlation is 0.822
  
• Regression line is \(y = 0.5x + 3.0\)
  

The variance is less than the target value of 4.125 by about 1.34%, and the correlation is greater than the target value (\(sqrt{(2/3)}\)) by about 0.7%. We’re not surprised by the regression line results, which are identical to the desired values (average of x = 8, average of y at 8 = 7, …). With just a little bit of tweaking, we could build a Y4 variable from this starting point that fulfills all the requirements.

While Anscombe’s emphasis was on the power of graphs, it’s interesting to note a few things about the data that don’t require a scatterplot.

If you sort the data (low to high) by X, it’s hard not to notice that X is the set of integers [4,14]. It may be that Anscombe simply intended his “observation order” to make the variables look more like an observed dataset than a contrived one, but from a calculation standpoint I would guess that he did most of his work in sorted order. If you sort the Figure 4 data by Y4, the data point (19,12.5) stands out. It’s also apparent that Y3 and Y4 have something of a similar feel to each other (I’m not showing the calculations, but if you sort Y3 and Y4 each from low to high the correlation coefficient is 0.995).

Also note that there are other summary statistics that tell us there are differences across the variables.

Kurtosis for X (Figures 1,2 and 3) is -1.2. This is a strong clue that the dataset is uniformly distributed in nature. The kurtosis for X4 is 11.0 and skewness is 3.32. For a set of data with 11 values, these are essentially the largest magnitude values possible, and strong clues that there’s an outlier in the dataset. Given X4’s standard deviation of 3.2, we could also note that the F-pseudosigma value (based on the interquartile range) of 0.0 suggests a dataset with interesting properties.

The Y variables have a range of skewness and kurtosis values. We did’t emphasize it in this report, but only Y3 has a mean value of exactly 7.5. The other Y variables sum to 82.51, so their means are ever so slightly greater than 7.5. For those interested, you will find that Y1, Y2, and Y4 each can be modified by subtracting 0.01 from a single data value while maintaining all the desired summary statistics, including regression and correlation.

“Graphs are essential to good statistical analysis.” It’s a safe bet that Anscombe made a strong case for his opening claim, as one needs only to search for “Anscombe’s Quartet” online to see a myriad of results. There may or may not be higher honors, but the first page of Edward Tufte’s legendary The Visual Display of Quantitiative Information shows the Anscombe data values; the second page shows the scatterplots. The Anscombe dataset is part of the R package datasets, while the first chapter of the NIST/SEMATECH e-Handbook of Statistical Methods uses the quartet to illustrate the benefits of exploratory data analysis and graphics. The quartet shows that a small set of data can still have an impact in a ‘big data’ world.