My friend Guy has been monitoring electricity consumption of his old refrigerator, and he wants to know if the refrigerator is still as efficient as it should be, or if it’s time to replace it.
Though the data set is small, it is real data and presents many of the challenges to a data scientist that any real data set will. Here I walk the reader through a complete analysis in R, step-by-step, from data entry to final recommendation.
When analyzing data, I have often found it useful to think of the data as being one of four main types, according to the typology proposed by Stevens.[1] Different types of data have certain characteristics; understanding what type of data you have helps with selecting the analysis to perform perform while preventing basic mistakes.
The types, or “scales of measurement,” are:
The generally-appropriate statistics and mathematical operations for each type are summarized in table 1.
While this is a useful typology for most use, and certainly for initial consideration, there are valid criticisms of Stevens’ typology. For example, percentages and count data have some characteristics of ratio-scale data, but with additional constraints. e.g. the average of the counts 
R[2] recognizes at least fifteen different types of data. Several of these are related to identifying functions and other objects—most users don’t need to worry about most of them. The main types that industrial engineers and scientists will need to use are:
x <- 1 + i2
or
x <- complex(real=1, imaginary=2)
Like type numeric, may be used for all scales of measurement.
The above types of data can be stored in several types, or structures, of variables. The equivalent to a variable in Excel would be rows, columns or tables of data. The main ones that we will use are:
c(1, 2, 1, 3, 4) ## [1] 1 2 1 3 4
is, by default, a numeric vector of type double, but
c(1, 2, 1, 3, 4, "name") ## [1] "1" "2" "1" "3" "4" "name"
will be a character vector, or a vector where all data is stored as type character, and the numbers will be stored as characters rather than numbers. It will not be possible to perform mathematical operations on these numbers-stored-as-characters without first converting them to type numeric.
factor(c("a", "b", "c", "a"), levels=c("a","b","c","d"))
## [1] a b c a
## Levels: a b c d
# letters a - c in 2x4 array array(data=letters[1:3], dim=c(2,4)) ## [,1] [,2] [,3] [,4] ## [1,] "a" "c" "b" "a" ## [2,] "b" "a" "c" "b" # numbers 1 - 3 in 2x4 array array(data=1:3, dim=c(2,4)) ## [,1] [,2] [,3] [,4] ## [1,] 1 3 2 1 ## [2,] 2 1 3 2
matrix(data = rep(1:3, times=2), nrow=2, ncol=3) ## [,1] [,2] [,3] ## [1,] 1 3 2 ## [2,] 2 1 3
list("text", "more", 2, c(1,2,3,2))
## [[1]]
## [1] "text"
##
## [[2]]
## [1] "more"
##
## [[3]]
## [1] 2
##
## [[4]]
## [1] 1 2 3 2
list where all columns must be vectors of the same number of rows (determined with NROW()). However, unlike matrices, different columns can contain different types of data and each row and column must have a name. If not named explicitly, R names rows by their row number and columns according to the data assigned assigned to the column. Data frames are typically used to store the sort of data that industrial engineers and scientists most often work with, and is the closest analog in R to an Excel spreadsheet. Usually data frames are made up of one or more columns of factors and one or more columns of numeric data.
data.frame(rnorm(5), rnorm(5), rnorm(5)) ## rnorm.5. rnorm.5..1 rnorm.5..2 ## 1 0.2939566 1.28985202 -0.01669957 ## 2 0.3672161 -0.01663912 -1.02064116 ## 3 1.0871615 1.13855476 0.78573775 ## 4 -0.8501263 -0.17928722 1.03848796 ## 5 -1.6409403 -0.34025455 -0.62113545
More generally, in R all variables are objects, and R distinguishes between objects by their internal storage type and by their class declaration, which are accessible via the typeof() and class() functions. Functions in R are also objects, and the users can define new objects to control the output from functions like summary() and print(). For more on objects, types and classes, see section 2 of the R Language Definition.
Table 2 summarizes the internal storage and R classes of the main data and variable types.
Many of our statistical tests make assumptions about the distribution of the underlying population. Many of the most common—ImR (XmR) and XbarR control charts, ANOVA, t-tests—assume normal distributions in the underlying population (or normal distributions in the residuals, in the case of ANOVA), and we’re often told that we must carefully check the assumptions.
At the same time, there’s a lot of conflicting advice about how to test for normality. There are the statistical tests for normality, such as Shapiro-Wilk or Anderson-Darling. There’s the “fat pencil” test, where we just eye-ball the distribution and use our best judgement. We could even use control charts, as they’re designed to detect deviations from the expected distribution. We are discouraged from using the “fat pencil” because it will result in a lot of variation from person to person. We’re often told not to rely too heavily on the statistical tests because they are not sensitive with small sample sizes and too sensitive to the tails. In industrial settings, our data is often messy, and the tails are likely to be the least reliable portion of our data.
I’d like to explore what the above objections really look like. I’ll use R to generate some fake data based on the normal distribution and the t distribution, and compare the frequency of p-values obtained from the Shapiro-Wilk test for normality.
First, we need to load our libraries
library(ggplot2) library(reshape2)
To make this easy to run, I’ll create a function to perform a large number of normality tests (Shapiro-Wilk) for sample sizes n = 5, 10 and 1000, all drawn from the same data:
#' @name assign_vector
#' @param data A vector of data to perform the t-test on.
#' @param n An integer indicating the number of t-tests to perform. Default is 1000
#' @return A data frame in "tall" format
assign_vector <- function(data, n = 1000) {
# replicate the call to shapiro.test n times to build up a vector of p-values
p.5 <- replicate(n=n, expr=shapiro.test(sample(my.data, 5, replace=TRUE))$p.value)
p.10 <- replicate(n=n, expr=shapiro.test(sample(my.data, 10, replace=TRUE))$p.value)
p.1000 <- replicate(n=n, expr=shapiro.test(sample(my.data, 1000, replace=TRUE))$p.value)
#' Combine the data into a data frame,
#' one column for each number of samples tested.
p.df <- cbind(p.5, p.10, p.1000)
p.df <- as.data.frame(p.df)
colnames(p.df) <- c("5 samples","10 samples","1000 samples")
#' Put the data in "tall" format, one column for number of samples
#' and one column for the p-value.
p.df.m <- melt(p.df)
#' Make sure the levels are sorted correctly.
p.df.m <- transform(p.df.m, variable = factor(variable, levels = c("5 samples","10 samples","1000 samples")))
return(p.df.m)
}
I want to simulate real-word conditions, where we have an underlying population from which we sample a limited number of times. To start, I’ll generate 100000 values from a normal distribution. To keep runtimes low I’ll have assign_vector() sample from that distribution when performing the test for normality.
n.rand <- 100000 n.test <- 10000 my.data <- rnorm(n.rand) p.df.m <- assign_vector(my.data, n = n.test)
We would expect that normally distributed random data will have an equal probability of any given p-value. i.e. 5% of the time we’ll see p-value ≤ 0.05, 5% of the time we’ll see p-value > 0.05 and ≤ 0.10, and so on through > 0.95 and ≤ 1.00. Let’s graph that and see what we get for each sample size:
ggplot(p.df.m, aes(x = value)) +
geom_histogram(binwidth = 1/10) +
facet_grid(facets=variable ~ ., scales="free_y") +
xlim(0,1) +
ylab("Count of p-values") +
xlab("p-values") +
theme(text = element_text(size = 16))
Histogram of p-values for the normal distribution, for sample sizes 5, 10 and 1000.
This is, indeed, what we expected.
Now let’s compare the normal distribution to a t distribution. The t distribution would pass the “fat pencil” test—it looks normal to the eye:
ggplot(NULL, aes(x=x, colour = distribution)) +
stat_function(fun=dnorm, data = data.frame(x = c(-6,6), distribution = factor(1)), size = 1) +
stat_function(fun=dt, args = list( df = 20), data = data.frame(x = c(-6,6), distribution = factor(2)), linetype = "dashed", size = 1) +
scale_colour_manual(values = c("blue","red"), labels = c("Normal","T-Distribution")) +
theme(text = element_text(size = 12),
legend.position = c(0.85, 0.75)) +
xlim(-4, 4) +
xlab(NULL) +
ylab(NULL)
Starting with random data generated from the t-distribution:
my.data <- rt(n.rand, df = 20)
Histogram of p-values for the t distribution, for sample sizes 5, 10 and 1000.
The tests for normality are not very sensitive for small sample sizes, and are much more sensitive for large sample sizes. Even with a sample size of 1000, the data from a t distribution only fails the test for normality about 50% of the time (add up the frequencies for p-value > 0.05 to see this).
Since the t distribution is narrower in the middle range and has longer tails than the normal distribution, the normality test might be failing because the entire distribution doesn’t look quite normal; we haven’t learned anything specifically about the tails.
To test the tails, we can construct a data set that uses the t distribution for the middle 99% of the data, and the normal distribution for the tails.
my.data <- rt(n.rand, df = 20) my.data.2 <- rnorm(n.rand) # Trim off the tails my.data <- my.data[which(my.data < 3 & my.data > -3)] # Add in tails from the other distribution my.data <- c(my.data, my.data.2[which(my.data.2 < -3 | my.data.2 > 3)])
Histogram of p-values for sample sizes 5, 10 and 1000, from a data set constructed from the t distribution in the range -3 to +3 sigmas, with tails from the normal distribution below -3 and above +3.
Despite 99% of the data being from the t distribution, this is almost identical to our test with data from just the normal distribution. It looks like the tails may be having a larger impact on the normality test than rest of the data
Now let’s flip this around: data that is 99% normally-distributed, but using the t distribution in the extreme tails.
my.data <- rnorm(n.rand) my.data.2 <- rt(n.rand, df = 20) # Trim off the tails my.data <- my.data[which(my.data < 3 & my.data > -3)] # Add in tails from the other distribution my.data <- c(my.data, my.data.2[which(my.data.2 < -3 | my.data.2 > 3)])
Histogram of p-values for sample sizes 5, 10 and 1000, from a data set constructed from the normal distribution in the range -3 to +3 sigmas, with tails from the t-distribution below -3 and above +3.
Here, 99% of the data is from the normal distribution, yet the normality test looks almost the same as the normality test for just the t-distribution. If you check the y-axis scales carefully, you’ll see that the chance of getting p-value ≤ 0.05 is a bit lower here than for the t distribution.
To make the point further, suppose we have highly skewed data:
my.data <- rlnorm(n.rand, 0, 0.4)
This looks like:
For small sample sizes, even this is likely to pass a test for normality:
In other words, if we have enough data to fail a normality test, we always will because our real-world data won’t be clean enough. If we don’t have enough data to reliably fail a normality test, then there’s no point in performing the test, and we have to rely on the fat pencil test or our own understanding of the underlying processes.
Don’t get too hung up on whether your data is normally distributed or not. When evaluating and summarizing data, rely mainly on your brain and use the statistics only to catch really big errors in judgement. When attempting to make predictions about future performance, e.g. calculating Cpk or simulating a process, recognize the opportunities for errors in judgment and explicitly state you assumptions.
The free and open-source R statistics package is a great tool for data analysis. The free add-on package qcc provides a wide array of statistical process control charts and other quality tools, which can be used for monitoring and controlling industrial processes, business processes or data collection processes. It’s a great package and highly customizable, but the one feature I wanted was the ability to manipulate the control charts within the grid graphics system, and that turned out to be not so easy.
I went all-in and completely rewrote qcc’s plot.qcc() function to use Hadley Wickham’s ggplot2 package, which itself is built on top of grid graphics. I have tested the new code against all the examples provided on the qcc help page, and the new ggplot2 version works for all the plots, including X-bar and R, p- and u- and c-charts.
In qcc, an individuals and moving range (XmR or ImR) chart can be created simply:
library(qcc) my.xmr.raw <- c(5045,4350,4350,3975,4290,4430,4485,4285,3980,3925,3645,3760,3300,3685,3463,5200) x <- qcc(my.xmr.raw, type = "xbar.one", title = "Individuals Chart\nfor Wheeler sample data") x <- qcc(matrix(cbind(my.xmr.raw[1:length(my.xmr.raw)-1], my.xmr.raw[2:length(my.xmr.raw)]), ncol = 2), type = "R", title = "Moving Range Chart\nfor Wheeler sample data")
This both generates the plot and creates a qcc object, assigning it to the variable x. You can generate another copy of the plot with plot(x).
To use my new plot function, you will need to have the packages ggplot2, gtable, qcc and grid installed. Download my code from the qcc_ggplot project on Github, load qcc in R and then run source("qcc.plot.R"). The ggplot2-based version of the plotting function will be used whenever a qcc object is plotted.
library(qcc)
source("qcc.plot.R")
my.xmr.raw <- c(5045,4350,4350,3975,4290,4430,4485,4285,3980,3925,3645,3760,3300,3685,3463,5200)
x <- qcc(my.xmr.raw, type = "xbar.one", title = "Individuals Chart\nfor Wheeler sample data")
x <- qcc(matrix(cbind(my.xmr.raw[1:length(my.xmr.raw)-1], my.xmr.raw[2:length(my.xmr.raw)]), ncol = 2), type = "R", title = "Moving Range Chart\nfor Wheeler sample data")
Below, you can compare the individuals and moving range charts generated by qcc and by my new implementation of plot.qcc():
The qcc individuals chart as implemented in the qcc package.
The qcc individuals chart as implemented using ggplot2 and grid graphics.
The qcc moving range chart as implemented in the qcc package.
The qcc moving range chart as implemented using ggplot2 and grid graphics.
In addition to the standard features in qcc plots, I’ve added a few new options.
size or cexgeom_point().font.sizeggplot() and grid’s viewport().title = element_blank()title = NULL a title will be created, or a user-defined text string may be passed to title.TRUE, creates a new graph (grid.newpage()). Otherwise, will write into the existing device and viewport. Intended to simplify the creation of multi-panel or composite charts.digits is provided by the qcc package to control the number of digits printed on the graph, where it either uses the default option set for R or a user-supplied value. I have tried to add some intelligence to calculating a default value under the assumption that we can tell something about the measurement from the data supplied. You can see the results in the sample graphs above.This little project turned out to be somewhat more difficult than I had envisioned, and there are several lessons-learned, particularly in the use of ggplot2.
First, ggplot2 really needs data frames when plotting. Passing discrete values or variables not connected to a data frame will often result in errors or just incorrect results. This is different than either base graphics or grid graphics, and while Hadley Wickham has mentioned this before, I hadn’t fully appreciated it. For instance, this doesn’t work very well:
my.test.data <- data.frame(x = seq(1:10), y = round(runif(10, 100, 300))) my.test.gplot <- ggplot(my.test.data, aes(x = x, y = y)) + geom_point(shape = 20) index.1 <- c(5, 6, 7) my.test.gplot <- my.test.gplot + geom_point(aes(x = x[index.1], y = y[index.1]), col = "red") my.test.gplot
Different variations of this sometimes worked, or sometimes only plotted some of the points that are supposed to be colored red.
However, if I wrap that index.1 into a data frame, it works perfectly:
my.test.data <- data.frame(x = seq(1:10), y = round(runif(10, 100, 300))) my.test.gplot <- ggplot(my.test.data, aes(x = x, y = y)) + geom_point(shape = 20) index.1 <- c(5, 6, 7) my.test.subdata <- my.test.data[index.1,] my.test.gplot <- my.test.gplot + geom_point(data = my.test.subdata, aes(x = x, y = y), col = "red") my.test.gplot
Another nice lesson was that aes() doesn’t always work properly when ggplot2 is called from within a function. In this case, aes_string() usually works. There’s less documentation than I would like on this, but you can search the ggplot2 Google Group or Stack Overflow for more information.
One of the bigger surprises was discovering that aes() searches for data frames in the global environment. When ggplot() is used from within a function, though, any variables created within that function are not accessible in the global environment. The work-around is to tell ggplot which environment to search in, and a simple addition of environment = environment() within the ggplot() call seems to do the trick. This is captured in a stack overflow post and the ggplot2 issue log.
my.test.data <- data.frame(x = seq(1:10), y = round(runif(10, 100, 300))) my.test.gplot <- ggplot(my.test.data, environment = environment(), aes(x = x, y = y)) + geom_point(shape = 20) index.1 <- c(5, 6, 7) my.test.subdata <- my.test.data[index.1,] my.test.gplot <- my.test.gplot + geom_point(data = my.test.subdata, aes(x = x, y = y), col = "blue") my.test.gplot
Finally, it is possible to completely and seamlessly replace a function created in a package and loaded in that package’s namespace. When I set out, I wanted to end up with a complete replacement for qcc’s internal plot.qcc() function, but wasn’t quite sure this would be possible. Luckily, the below code, called after the function declaration, worked. One thing I found was that I needed to name my function the same as the one in the qcc package in order for the replacement to work in all cases. If I used a different name for my function, it would work when I called plot() with a qcc object, but qcc’s base graphics version would be used when calling qcc() with the parameter plot = TRUE.
unlockBinding(sym="plot.qcc", env=getNamespace("qcc"));
assignInNamespace(x="plot.qcc", value=plot.qcc, ns=asNamespace("qcc"), envir=getNamespace("qcc"));
assign("plot.qcc", plot.qcc, envir=getNamespace("qcc"));
lockBinding(sym="plot.qcc", env=getNamespace("qcc"));
For now, the code suits my immediate needs, and I hope that you will find it useful. I have some ideas for additional features that I may implement in the future. There are some parts of the code that can and should be further cleaned up, and I’ll tweak the code as needed. I am certainly interested in any bug reports and in seeing any forks; good ideas are always welcome.
Updated below
The R Function Reference is a mind map that I created as a guide for novice and intermediate users of the R statistics language. When you first open it, I suggest that you collapse all the nodes by clicking on the “Expand/Collapse all nodes” button in the bottom left of the screen to make the map easier to navigate. You can also adjust the zoom level with the slider next to that button.
The top-level nodes of the R Function Reference
The mind map is arranged in eight sections, or main branches, arranged by task. What do you want to do? Each branch covers a general set tasks, such as learning to use R, running R, working with data, statistical analysis or plotting data. The end of each string of nodes is generally a function and example. The Reference provides code fragments, rather than details of the function or complete reproducible code blocks. Once you’ve followed the Reference and have an idea of how to accomplish something, you can look up the details in R’s help system (e.g. “?read.csv” to learn more about using the read.csv() function), or search Google or the online R-Help mailing list archives for answers using the function name.
There are a lot of useful nodes and examples, especially in the “Graphs” section, but the mind map is not complete; some trails end before you get to a useful function reference. I am sorry for that, but it’s a work in progress, and will be slowly updated over time.
Comments and suggestions are welcome.
In comments, several users reported problems opening the mind map. With a little investigation, it appears that the size of the mind map is the problem. To try to fix the problems , I have split the mind map out into several small mind maps, all linked together.
The new main mind map is the R Function Reference, Main. The larger branches on this main map no longer expand to their own content, but contain a link to a “child” mind map. The link looks like a sheet of paper with an arrow pointing to the right, click on it and little cartoon speech bubble will pop up with a link that you have to click on to go to the child mind map. Likewise, the central nodes on the child mind maps contain a link back to the main mind map.
Due to load times and the required extra clicks, this may slightly reduce usability for users who didn’t have a problem with the all-in-one version, but will hopefully make the mind map accessible to a broader audience.
I have to offer praise to the developers of Mind42. Though I couldn’t directly split branches off into their own mind maps or duplicate the mind map, it was very easy to export the mind map as a native Mind42 file and then import it multiple times, editing the copies without any loss of data or links. The ability to link directly between mind maps within Mind42 was also a key enabling feature. Considering that this is a free web app, its capabilities are most impressive. They were also quick to respond when I posted a call for help on the Mind42 forum.
Please let me know how the new, “improved” version works.
The old mind map, containing everything, is still available, but I will not update it.
Recently Chandoo.org posted a question about how to graph data when you have a lot of small values and a few larger values. It’s not the first time that I’ve come across this question, and I’ve seen a lot of answers, many of them really bad. While all solutions involve trade-offs for understanding and interpreting graphs, some solutions are better than others.
Data graphs tell stories by revealing patterns in complex data. Good data graphs let the data tell the story by revealing the patterns, rather than trying to impose patterns on the data.
As William Cleveland discusses in The Elements of Graphing Data and his 1993 paper A Model for Studying Display Methods of Statistical Graphics, there are two basic visual operations that people employ when looking at and interpreting graphs: pattern perception and table look-up. Pattern perception is where we see the geometric patterns in a graph: groupings; relative differences (larger/smaller); or trends (straight/curved or increasing/decreasing). Table look-up is where we explore details of values and names on a graph. These two operations are distinct and complimentary, and it is through these two operations that the data’s story is told.
| month | sales | |
|---|---|---|
| 1 | Feb 09 | 200 |
| 2 | Mar 09 | 300 |
| 3 | Apr 09 | 200 |
| 4 | May 09 | 300 |
| 5 | Jun 09 | 200 |
| 6 | Jul 09 | 300 |
| 7 | Aug 09 | 350 |
| 8 | Sep 09 | 400 |
| 9 | Oct 09 | 450 |
| 10 | Nov 09 | 1200 |
| 11 | Dec 09 | 100000 |
| 12 | Jan 10 | 85000 |
| 13 | Feb 10 | 450 |
So suppose that we have some data like that at right, where we are interested in the patterns of smaller, individual values, but there are also a few extremely large values, or outliers. We describe such data as being skewed. How do we plot this data? First, for such a small data set, a simple table is the best approach. People can see the numbers and interpret them, there aren’t too many numbers to make sense of and the table is very compact. For more complicated data sets, though, a graph is needed. There’s a few basic options:
A bar chart with all data, including outliers, plotted on the same scale.
This is the simplest solution, and if you’re only interested in knowing about the outliers (Dec ’09 and Jan ’10) then it will do. However, it completely hides whatever is happening in the rest of the months. Pattern recognition tells us that two months near the end of the series have the big numbers. Table-lookup tells us the approximate values and that these months are around December ’09 and February ’10, but the way the labels string together and overlap the tick marks, it’s not clear exactly what the labels are, let alone which label applies to which bar (which months are those, precisely? Is that “09 Dec” and “09 Feb?” Do the numbers even go with the text, or are they separate labels?).
For all but the simplest of messages, this rendition defeats both pattern recognition and table look-up. We definitely need a better solution.
Excel gives us an easy solution: break the data into two columns (“small” numbers in one and “large” numbers in the other) and plot them on separate axes. Now we can see all the data, including the patterns in all the months.
Bar chart, with outliers plotted using a secondary axis.
Unfortunately, pattern recognition tells us that the big-sales months are about the same as all the other months. It’s only the table look-up that tells us how big of a difference there is between the two blue columns and the rest of the data. This is why I’ve added data labels to the two columns: to aid table look-up.
Even if we tweaked around with the axes to set the outliers off from the rest of the data, we’d still have the same basic problem: pattern recognition would tell us that there is a much smaller difference than there actually is. By using a secondary axis, we’ve set up a basic conflict between pattern recognition and table look-up. Worse, it’s easy to confuse the axes; which bars go with which axis? Reproduction in black and white or grayscale would make it impossible to correctly connect bars to the correct axis. Some types of color blindness would similarly make it difficult to interpret the graph. Table look-up is easily defeated with secondary axes.
The secondary axis presents so many problems that I always advise against using it. Stephen Few, author of Show Me The Numbers and Information Dashboard Design, calls graphs with secondary axes “dual-scaled graphs.” In his 2008 article Dual-Scaled Axes in Graphs, he concludes that there is always a better way to display data than by using secondary axes. Excel makes it easy to create graphs like this, but it’s always a bad idea.
In scientific applications, skewed data is common, and the usual solution is to plot the logarithm of the values.
Bar chart plotting skewed with logarithmic axis.
With the logarithm, it is easy to plot, and see, all of the data. Trends in small values are not hidden. Pattern perception immediately tells us the overall story of the data. Table look-up is easier than with secondary axes, and immediately tells us the scale of the differences. Plotting the logarithm allows pattern perception and table look-up to compliment each other.
Below, I’ve created the same graph using a dot plot instead of a bar chart. Dot plots have many advantages over bar charts: most obviously, dot plots provide a better arrangement for category labels (e.g. the months); also, dot plots provide a clearer view of the data by plotting the data points rather than filling in the space between the axis and the data point. There are some nice introductions to dot plots, including William Cleveland’s works and a short introduction by Naomi Robbins. The message is clear: any data that you might present with a bar chart (or pie chart) will be better presented using dot plots.
Skewed data plotted on a dot plot using a logarithmic scale.
Another approach, which might be better for audiences unfamiliar with logarithmic scales, is to use a scale break, or broken axis. With some work, we can create a scale break in Excel or OpenOffice.org.
Bar chart with outliers plotted by introducing a subtle scale break on the y-axis.
There are plenty of tutorials for how to accomplish this in Excel. For this example, I created the graph in OpenOffice.org Spreadsheet, using the same graph with the secondary axis, above. I adjusted the two scales, turned off the labels for both y-axes and turned off the tick marks for the secondary y-axis. Then I copied the graph over to the OpenOffice.org Draw application and added y-axis labels and the break marks as drawing objects.
That pretty much highlights the first problem with this approach: it takes a lot of work. The second problem is that those break marks are just too subtle; people will miss them.
The bigger problem is with interpretation. As with the secondary axis, this subtle scale break sets up a basic conflict between the two basic operations of graph interpretation. Pattern recognition tells us that the numbers are comparable; it’s only table look-up that tells us what a large difference there is.
Cleveland’s recommendation, when the logarithm won’t work, is to use a full-panel scale break. In this way, pattern recognition tells that there are two distinct groups of data, and table look-up tells us what they are.
Dot plot with a full scale break to show outliers.
The potential disadvantage of this approach is that pattern perception might be fooled. While the scale break visually groups the “large” values from the “small” ones, the scale also changes, so that the broader panel on the left actually represents a much narrower range of values (about 1100 dollars range) than the narrower panel on the right (about 17000 dollars range). Our audience might have difficulties interpreting this correctly.
Edward Tufte has popularized the idea of small multiples, the emphasis of differences by repeating a graph or image with small changes from one frame to the next. In this case, we could show the full data set, losing fidelity in the smaller values, and then repeat the graph while progressively zooming in on a narrower and narrower slice with each repetition.
The full data, with outliers, is plotted on the left. On the right, a zoomed view showing detail in the smaller values.
This shares many similarities to Cleveland’s full scale break, but provides greater flexibility. With this data, there are two natural ranges: 0 – 100000 and 0 – 1200. If there were more data between 1200 and 85000, we might repeat the graph several times, zooming in more with each repetition to show lower levels of detail.
I think there are two potential pitfalls. As with the full scale break, the audience might fail to appreciate the effect of the changes to scale. Worse, the audience might be fooled into thinking that each graph represented a different set of data, rather than just a different slice of the same data. Some care in preparing such graphs will be needed for successful communication.
When presenting data that is, like the data above, arranged by category, use a dot plot instead of bar charts. When your data is heavily skewed, the best solution is to graph the logarithm of the data. However, if your audience will be unable to correctly interpret the logarithm, try a full scale break or small multiples.
I’m reading a work of fiction about the Knights Templar, based on the same mythos as Dan Brown’s novels. A throw-away line about the history of the Order’s Masters sparked my mathematical curriosity: “For the sixty-six who’d come before, the average tenure was a mere eighteen years.”
Is this reasonable? Did the author calculate the average tenure or just guess? The Knights Templar got their start in 1119. The story takes place around 2000, maybe 2005. So do sixty-six Masters average eighteen-year tenures over some 885 years?
We can do a quick check: had there been eighty-eight Masters, the average tenure would be about ten years. Had there been forty-four Masters, the average tenure would have been twenty years. Sixty-six falls mid-way between the two, so the average should be around fifteen. Eighteen is not completely unreasonable, but it might be too long.
We can be more exact: divide 885 by 66 and we immediately see that the average should be about 13.5 years. But this calculation is really answering the question “how long would the tenures be if all the Masters had the same tenure?” We might expect some short tenures of a year or two and one or two long tenures of perhaps twenty or thirty years.
Such a distribution, with no values less than zero, a few large values and most values clustering somewhere inbetween, might have a very different average due to the lopsided (skewed) distribution. We can approximate such a distribution with the Poisson distribution. Poisson does not have fractional increments, so we’ll only get whole-year tenures, but it should be good enough to determine if 18 is a reasonable average tenure. Also, Poisson is easier to fit than the more precise beta distribution.
So I fired up Minitab, generated a bunch of random Poisson-distributed values with a mean of 18, and then added them up in groups of sixty-six. The average of these sums was 1224 years; much longer than the 885 years required by the story. Eighteen years is too long.
Playing around with different values for the mean, I find that the “right” average for tenure length should be close to 13.5.
To answer the original question, then: eighteen years isn’t completely unreasonable, but it’s definitly wrong. I have to wonder how the author came up with this. If he just pulled “eighteen” and “sixty-six” out of thin air, then I have to say that he guessed pretty well. Unfortunately, it’s clear that he didn’t bother to do even a simple calculation while sitting in front of his computer typing, where a calculator was readily available.
I know: mathematical accuracy isn’t the point of the story. However, the point of fiction isn’t to create a milieu based on logical falsehoods, but to create a fictional milieu that is believable (i.e. our modern world, if the myths surrounding the Templars were real). The author lets us down through such acts of carelessness or laziness.
This also brings me back to my title. You know you’re a geek when you take a throw-away statement in a work of fiction and perform some statistical analysis to fact-check it.
At least it’s a fun journey.
Tom Hopper 
You must be logged in to post a comment.