This tutorial makes use of the following R package(s): ggplot2, dplyr, lattice (for data only).

This material can be read in conjunction with section 2.4 of Cleveland’s book.

# Introduction

Thus far, we have used the quantile-quantile plots to compare the distributions between two empirical (i.e. observational) datasets. We can also use the q-q plot to compare an empirical observation to a theoretical observation (i.e. one defined mathematically). Examples of popular theoretical observations are the normal distribution (aka the Gaussian distribution), the chi-square distribution, and the exponential distribution just to name a few.

There are many reasons we might want to compare empirical data to theoretical distributions:

• A theoretical distribution is easy to parameterize. For example, if the shape of the distribution of a batch of numbers can be approximated by a normal distribution we can reduce the complexity of our data to just two values: the mean and the standard deviation.

• If data can be approximated by certain theoretical distributions, then many mainstream statistical procedures can be applied to the data.

• In inferential statistics, knowing that a sample was derived from a population whose distribution follows a theoretical distribution allows us to derive certain properties of the population from the sample. For example, if we know that a sample comes from a normally distributed population, we can define confidence intervals for the sample mean using a t-distribution.

• Modeling the distribution of the observed data can provide insight into the underlying process that generated the data.

But very few empirical datasets follow any theoretical distributions exactly. So the questions usually ends up being “how well does theoretical distribution X fit my data?”

The theoretical quantile-quantile plot is a tool to explore how a batch of numbers deviates from a theoretical distribution and to visually assess whether the difference is significant for the purpose of the analysis. In the following examples, we will compare empirical data to the normal distribution using the normal quantile-quantile plot.

# The normal q-q plot

The normal q-q plot is just a special case of the empirical q-q plot we’ve explored so far; the difference being that we assign the normal distribution quantiles to the x-axis.

## Drawing a normal q-q plot from scratch

In the following example, we’ll compare the Alto 1 group to a normal distribution. First, we’ll extract the Alto 1 height values and save them as an atomic vector object using dplyr’s piping operations. However, dplyr’s operations will return a dataframe–even if a single column is select. To force the output to an atomic vector, we’ll pipe the subset to .$height which will extract the height column into a plain vector element. library(dplyr) df <- lattice::singer alto <- df %>% filter(voice.part == "Alto 1") %>% .$height

Next, we need to sort alto in ascending order.

alto <- sort(alto)

Next, we need to find the matching normal distribution quantiles. We first find the quantiles for alto, then use qnorm to find the matching normal distribution values from the quantiles.

i      <- 1:length(alto)
fi     <- (i - 0.5) / length(alto)
x.norm <- qnorm(fi)

Now we can plot the sorted alto values against the normal values.

plot( alto ~ x.norm, type="p", xlab="Normal quantiles", pch=20)

When comparing a batch of numbers to a theoretical distribution on a q-q plot, we are looking for significant deviation from a straight line. To make it easier to judge straightness, we can fit a line to the points. Note that we are not creating a 45° (or x=y) slope; the range of values between both sets of numbers do not match. Here, we are only seeking the straightness of the points.

There are many ways one can fit a line to the data, Cleveland opts to fit a line to the first and third quartile of the q-q plot. The following chunk of code identifies the quantiles for both the alto dataset and the theoretical normal distribution. It then computes the slope and intercept for that line.

# Find 1st and 3rd quartile for the Alto 1 data
y <- quantile(alto, c(0.25, 0.75))

# Find the 1st and 3rd quartile of the normal distribution
x <- qnorm( c(0.25, 0.75))

# Now we can compute the intercept and slope of the line that passes
# through these points
slope <- diff(y) / diff(x)
int   <- y[1] - slope * x[1]

Next, we add the line to the plot.

abline(a=int, b=slope )

## Using R’s built-in functions

R has two built-in functions that facilitate the plot building task when comparing a batch to a normal distribution: qqnorm and qqline.

qqnorm(alto)
qqline(alto)

That’s it. Just two lines of code!

## Using the ggplot2 plotting environment

We can take advantage of the stat_qq() function to plot the points, but the equation for the line must be computed manually (as was done earlier). Those steps will be repeated here.

library(ggplot2)

# Find the slope and intercept of the line that passes through the 1st and 3rd
# quartile of the normal q-q plot

y     <- quantile(alto, c(0.25, 0.75)) # Find the 1st and 3rd quartiles
x     <- qnorm( c(0.25, 0.75))         # Find the matching normal values on the x-axis
slope <- diff(y) / diff(x)             # Compute the line slope
int   <- y[1] - slope * x[1]           # Compute the line intercept

# Generate normal q-q plot

ggplot() + aes(sample=alto) + stat_qq(distribution=qnorm) +
geom_abline(intercept=int, slope=slope) + ylab("Height")           

We can, of course, make use of ggplot’s faceting function to generate trellised plots. For example, the following plot replicates Cleveland’s figure 2.11 (except for the layout which we’ll setup as a single row of plots instead). But first, we will need to compute the slopes for each singer group. We’ll use dplyr’s piping operations to create a new dataframe with singer group name, slope and intercept.

library(dplyr)
intsl <- df %>% group_by(voice.part) %>%
summarize(q25    = quantile(height,0.25),
q75    = quantile(height,0.75),
norm25 = qnorm( 0.25),
norm75 = qnorm( 0.75),
slope  = (q25 - q75) / (norm25 - norm75),
int    = q25 - slope * norm25) %>%
select(voice.part,slope, int) 

The above chunk creates the following table.

voice.part slope int
Bass 2 2.965204 72.000
Bass 1 2.223903 70.500
Tenor 2 1.482602 70.000
Tenor 1 3.706506 68.500
Alto 2 2.223903 65.500
Alto 1 2.594554 64.750
Soprano 2 2.779879 63.875
Soprano 1 1.667927 63.875

It’s important that the voice.part names match those in df letter-for-letter so that when ggplot is called, it will know which facet to assign the slope and intercept values to via geom_abline.

ggplot(df, aes(sample=height)) + stat_qq(distribution=qnorm) +
geom_abline(data=intsl, aes(intercept=int, slope=slope), col="blue") +
facet_wrap(~voice.part,nrow=1) + ylab("Height")