In a recent paper we included data from a survey we conducted. During the publication process, one of the reviewers asked for a more in depth statistical analysis of the data set. He (she?) explicitly expressed a special interest in correlating the variables of the survey in order to spot any interesting correlations. This posed a number of problems:

1. The data set has 35 variables. This translates into a huge correlation matrix.

2. There are many ways to correlate variables. In my field the Pearson product-moment correlation coefficient (Pearson’s r) and the Spearman’s rank correlation coefficient (Spearman’s rho) are the most common ones. Which one should I choose and why?

3. Correlations by themselves are not very useful. You, most likely, want to have at least some information about the statistical significance of the correlation.

The Data Set

For the sake of demonstration we will use the mtcars data set, provided by the stats package. The data set includes 11 variables with 32 observations, hence rendering the task a bit more manageable.

Choosing a Suitable Correlation Algorithm

As most papers in my field use either Pearson’s r or Spearman’s rho I began my research there. As it turns out, there are two types of algorithms: parametric and non-parametric. Parametric algorithms depend on some assumptions about the data set. In case of Pearson’s r you need to know about the distribution, scaling and dependency of the variables of the data set.

Spearman’s rho, on the other hand, is a non-parametric algorithm. It compares two ranked variables of either ordinal, interval or ratio type and describes how well they can be described using a monotonic function. No further assumptions on the data set are necessary which is handy if you do not know too much about the distribution of your data.

Let’s Get Started!

R, a statistical computing programming language, is being used for the statistical analysis with some additional libraries: ggplot2 for plotting, Hmisc to create a correlation matrix, reshape2 to meld the data frame as well as stats to provide the mtcars data set.

library(ggplot2)
library(reshape2)
library(Hmisc)
library(stats)

The first thing is to calculate the correlation of the individual variables to another. As we have our data already in a data frame, merely a call of the rcorr function of the Hmisc library is necessary. This creates a new list with two entries: ”r” the correlation coefficients and ”P” the significance levels.

rcorr Computes a matrix of Pearson’s r or Spearman’s rho rank correlation coefficients for all possible pairs of columns of a matrix. Missing values are deleted in pairs rather than deleting all rows of x having any missing variables. Ranks are computed using efficient algorithms, using midranks for ties.
— R Documentation / rcorr {Hmisc}

d <- mtcars
cormatrix = rcorr(as.matrix(d), type='spearman')

The correlation coefficients can be plotted using a heat map representation. ggplot2 provides the geom_tile geometric object for this purpose. In order to plot the correlation matrix we need to meld the data frame first.

cordata = melt(cormatrix$r)
ggplot(cordata, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile() + xlab("") + ylab("")

Figure 2 shows the resulting heat map, plotted by ggplot2. However, some key data is missing. There is no information about the significance levels nor does the plot include the numeric values of Spearman's rho. In order to efficiently use the little space per tile, I wrote a little function ("abbreviateSTR") which abbreviates the values. The resulting strings are stored in a label column. This column can then be plotted onto the corresponding tile. To further enhance the distinction between significant (P > 0.05) and insignificant correlations, a red, semi transparent, ”X” is being overlaid onto insignificant tiles. The resulting plot is shown in figure 1.

Conclusion

By using the non-parametric Spearman algorithm it is possible to create a preliminary overview of the correlations within a data set. This information can be plotted utilizing a heat map representation in conjunction with significance levels as well as numeric values of Spearman's rho. The resulting plot is a comprehensible representation of the data set, which allows quick identification of significant correlations.

Looking at the heat map, one could get the impression that the gross horsepower has a strong negative correlation (r ≈ -0.89, P < 0.01) with the gasoline consumption. Who would have thought?

Complete R Script

#!/usr/bin/env Rscript

library(ggplot2)
library(reshape2)
library(Hmisc)
library(stats)

abbreviateSTR <- function(value, prefix){  # format string more concisely
  lst = c()
  for (item in value) {
    if (is.nan(item) || is.na(item)) { # if item is NaN return empty string
      lst <- c(lst, '')
      next
    }
    item <- round(item, 2) # round to two digits
    if (item == 0) { # if rounding results in 0 clarify
      item = '<.01'
    }
    item <- as.character(item)
    item <- sub("(^[0])+", "", item)    # remove leading 0: 0.05 -> .05
    item <- sub("(^-[0])+", "-", item)  # remove leading -0: -0.05 -> -.05
    lst <- c(lst, paste(prefix, item, sep = ""))
  }
  return(lst)
}

d <- mtcars

cormatrix = rcorr(as.matrix(d), type='spearman')
cordata = melt(cormatrix$r)
cordata$labelr = abbreviateSTR(melt(cormatrix$r)$value, 'r')
cordata$labelP = abbreviateSTR(melt(cormatrix$P)$value, 'P')
cordata$label = paste(cordata$labelr, "n", 
                      cordata$labelP, sep = "")
cordata$strike = ""
cordata$strike[cormatrix$P > 0.05] = "X"

txtsize <- par('din')[2] / 2
ggplot(cordata, aes(x=Var1, y=Var2, fill=value)) + geom_tile() + 
  theme(axis.text.x = element_text(angle=90, hjust=TRUE)) +
  xlab("") + ylab("") + 
  geom_text(label=cordata$label, size=txtsize) + 
  geom_text(label=cordata$strike, size=txtsize * 4, color="red", alpha=0.4)