Hands-on 9.3

Author

Guan Jhen LIN

Published

June 11, 2024

Modified

June 11, 2024

Hands on 9 consists of 5 different exercises, topics covered in this exercise:

9.3 Heatmap for Visualising and Analysing Multivariate Data

When applied to a tabular format, heatmaps are useful for cross-examining multivariate data, through placing variables in the columns and observation (or records) in rowa and colouring the cells within the table. Heatmaps are good for showing variance across multiple variables, revealing any patterns, displaying whether any variables are similar to each other, and for detecting if any correlations exist in-between them.

pacman::p_load(seriation, dendextend, heatmaply, tidyverse)

In this hands-on exercise, the data of World Happines 2018 report will be used. The data set is downloaded from here. The original data set is in Microsoft Excel format. It has been extracted and saved in csv file called WHData-2018.csv.

wh <- read_csv("data/WHData-2018.csv")
row.names(wh) <- wh$Country
wh1 <- dplyr::select(wh, c(3, 7:12))
wh_matrix <- data.matrix(wh)

Notice that wh_matrix is in R matrix format.

9.3.1 Static Heatmap

There are many R packages and functions can be used to drawing static heatmaps, they are:

  • heatmap()of R stats package. It draws a simple heatmap.

  • heatmap.2() of gplots R package. It draws an enhanced heatmap compared to the R base function.

  • pheatmap() of pheatmap R package. pheatmap package also known as Pretty Heatmap. The package provides functions to draws pretty heatmaps and provides more control to change the appearance of heatmaps.

  • ComplexHeatmap package of R/Bioconductor package. The package draws, annotates and arranges complex heatmaps (very useful for genomic data analysis). The full reference guide of the package is available here.

  • superheat package: A Graphical Tool for Exploring Complex Datasets Using Heatmaps. A system for generating extendable and customizable heatmaps for exploring complex datasets, including big data and data with multiple data types. The full reference guide of the package is available here.

In this section, you will learn how to plot static heatmaps by using heatmap() of R Stats package.

heatmap() of Base Stats

wh_heatmap <- heatmap(wh_matrix,
                      Rowv=NA, Colv=NA)

Note:

By default, heatmap() plots a cluster heatmap. The arguments Rowv=NA and Colv=NA are used to switch off the option of plotting the row and column dendrograms. To plot a cluster heatmap, we just have to use the default as shown in the code chunk below.

wh_heatmap <- heatmap(wh_matrix)

Note:

The order of both rows and columns is different compare to the native wh_matrix. This is because heatmap do a reordering using clusterisation: it calculates the distance between each pair of rows and columns and try to order them by similarity. Moreover, the corresponding dendrogram are provided beside the heatmap. Here, red cells denotes small values, and red small ones. This heatmap is not really informative. Indeed, the Happiness Score variable have relatively higher values, what makes that the other variables with small values all look the same. Thus, we need to normalize this matrix. This is done using the scale argument. It can be applied to rows or to columns following your needs.

The code chunk below normalises the matrix column-wise.

wh_heatmap <- heatmap(wh_matrix,
                      scale="column",
                      cexRow = 0.6, 
                      cexCol = 0.8,
                      margins = c(10, 4))

Notice that the values are scaled now. Also note that margins argument is used to ensure that the entire x-axis labels are displayed completely and, cexRow and cexCol arguments are used to define the font size used for y-axis and x-axis labels respectively.

9.3.2 Creating Interactive Heatmap

heatmaply is an R package for building interactive cluster heatmap that can be shared online as a stand-alone HTML file. It is designed and maintained by Tal Galili.

Before we get started, you should review the Introduction to Heatmaply to have an overall understanding of the features and functions of Heatmaply package. You are also required to have the user manualof the package handy with you for reference purposes.

In this section, you will gain hands-on experience on using heatmaply to design an interactive cluster heatmap. We will still use the wh_matrix as the input data.

9.3.3 Working with heatmaply

heatmaply(mtcars)
heatmaply(wh_matrix[, -c(1, 2, 4, 5)])

Note that:

Different from heatmap(), for heatmaply() the default horizontal dendrogram is placed on the left hand side of the heatmap. The text label of each raw, on the other hand, is placed on the right hand side of the heat map. When the x-axis marker labels are too long, they will be rotated by 135 degree from the north.

9.3.4 Data trasformation

When analysing multivariate data set, it is very common that the variables in the data sets includes values that reflect different types of measurement. In general, these variables’ values have their own range. In order to ensure that all the variables have comparable values, data transformation are commonly used before clustering.

Three main data transformation methods are supported by heatmaply(), namely: scale, normalise and percentilse.

9.3.5 Scaling method

When all variables are came from or assumed to come from some normal distribution, then scaling (i.e.: subtract the mean and divide by the standard deviation) would bring them all close to the standard normal distribution. In such a case, each value would reflect the distance from the mean in units of standard deviation. The scale argument in heatmaply() supports column and row scaling. The code chunk below is used to scale variable values columewise.

heatmaply(wh_matrix[, -c(1, 2, 4, 5)],
          scale = "column")

9.3.6 Normalising method

When variables in the data comes from possibly different (and non-normal) distributions, the normalize function can be used to bring data to the 0 to 1 scale by subtracting the minimum and dividing by the maximum of all observations. This preserves the shape of each variable’s distribution while making them easily comparable on the same “scale”. Different from Scaling, the normalise method is performed on the input data set i.e. wh_matrix as shown in the code chunk below.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]))

9.3.7 Percentising method

This is similar to ranking the variables, but instead of keeping the rank values, divide them by the maximal rank. This is done by using the ecdf of the variables on their own values, bringing each value to its empirical percentile. The benefit of the percentize function is that each value has a relatively clear interpretation, it is the percent of observations that got that value or below it. Similar to Normalize method, the Percentize method is also performed on the input data set i.e. wh_matrix as shown in the code chunk below.

heatmaply(percentize(wh_matrix[, -c(1, 2, 4, 5)]))

9.3.8 Clustering algorithm

heatmaply supports a variety of hierarchical clustering algorithm. The main arguments provided are:

  • distfun: function used to compute the distance (dissimilarity) between both rows and columns. Defaults to dist. The options “pearson”, “spearman” and “kendall” can be used to use correlation-based clustering, which uses as.dist(1 - cor(t(x))) as the distance metric (using the specified correlation method).

  • hclustfun: function used to compute the hierarchical clustering when Rowv or Colv are not dendrograms. Defaults to hclust.

  • dist_method default is NULL, which results in “euclidean” to be used. It can accept alternative character strings indicating the method to be passed to distfun. By default distfun is “dist”” hence this can be one of “euclidean”, “maximum”, “manhattan”, “canberra”, “binary” or “minkowski”.

  • hclust_method default is NULL, which results in “complete” method to be used. It can accept alternative character strings indicating the method to be passed to hclustfun. By default hclustfun is hclust hence this can be one of “ward.D”, “ward.D2”, “single”, “complete”, “average” (= UPGMA), “mcquitty” (= WPGMA), “median” (= WPGMC) or “centroid” (= UPGMC).

In general, a clustering model can be calibrated either manually or statistically.

9.3.9 Manual approach

In the code chunk below, the heatmap is plotted by using hierachical clustering algorithm with “Euclidean distance” and “ward.D” method.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          dist_method = "euclidean",
          hclust_method = "ward.D")

9.3.10 Statistical approach

In order to determine the best clustering method and number of cluster the dend_expend() and find_k() functions of dendextend package will be used.

First, the dend_expend() will be used to determine the recommended clustering method to be used.

wh_d <- dist(normalize(wh_matrix[, -c(1, 2, 4, 5)]), method = "euclidean")
dend_expend(wh_d)[[3]]
  dist_methods hclust_methods     optim
1      unknown         ward.D 0.6137851
2      unknown        ward.D2 0.6289186
3      unknown         single 0.4774362
4      unknown       complete 0.6434009
5      unknown        average 0.6701688
6      unknown       mcquitty 0.5020102
7      unknown         median 0.5901833
8      unknown       centroid 0.6338734

The output table shows that “average” method should be used because it gave the high optimum value.

Next, find_k() is used to determine the optimal number of cluster.

wh_clust <- hclust(wh_d, method = "average")
num_k <- find_k(wh_clust)
plot(num_k)

Figure above shows that k=3 would be good.

With reference to the statistical analysis results, we can prepare the code chunk as shown below.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          dist_method = "euclidean",
          hclust_method = "average",
          k_row = 3)

9.3.11 Seriation

One of the problems with hierarchical clustering is that it doesn’t actually place the rows in a definite order, it merely constrains the space of possible orderings. Take three items A, B and C. If you ignore reflections, there are three possible orderings: ABC, ACB, BAC. If clustering them gives you ((A+B)+C) as a tree, you know that C can’t end up between A and B, but it doesn’t tell you which way to flip the A+B cluster. It doesn’t tell you if the ABC ordering will lead to a clearer-looking heatmap than the BAC ordering.

heatmaply uses the seriation package to find an optimal ordering of rows and columns. Optimal means to optimize the Hamiltonian path length that is restricted by the dendrogram structure. This, in other words, means to rotate the branches so that the sum of distances between each adjacent leaf (label) will be minimized. This is related to a restricted version of the travelling salesman problem.

Here we meet our first seriation algorithm: Optimal Leaf Ordering (OLO). This algorithm starts with the output of an agglomerative clustering algorithm and produces a unique ordering, one that flips the various branches of the dendrogram around so as to minimize the sum of dissimilarities between adjacent leaves. Here is the result of applying Optimal Leaf Ordering to the same clustering result as the heatmap above.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          seriate = "OLO")

The default options is “OLO” (Optimal leaf ordering) which optimizes the above criterion (in O(n^4)). Another option is “GW” (Gruvaeus and Wainer) which aims for the same goal but uses a potentially faster heuristic.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          seriate = "GW")

The option “mean” gives the output we would get by default from heatmap functions in other packages such as gplots::heatmap.2.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          seriate = "mean")

The option “none” gives us the dendrograms without any rotation that is based on the data matrix.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          seriate = "none")

9.3.12 Working with colour palettes

The default colour palette uses by heatmaply is viridis. heatmaply users, however, can use other colour palettes in order to improve the aestheticness and visual friendliness of the heatmap.

In the code chunk below, the Blues colour palette of rColorBrewer is used

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          seriate = "none",
          colors = Blues)

9.3.13 The finishing touch

Beside providing a wide collection of arguments for meeting the statistical analysis needs, heatmaply also provides many plotting features to ensure cartographic quality heatmap can be produced.

In the code chunk below the following arguments are used:

k_row is used to produce 5 groups. margins is used to change the top margin to 60 and row margin to 200. fontsizw_row and fontsize_col are used to change the font size for row and column labels to 4. main is used to write the main title of the plot. xlab and ylab are used to write the x-axis and y-axis labels respectively.

heatmaply(normalize(wh_matrix[, -c(1, 2, 4, 5)]),
          Colv=NA,
          seriate = "none",
          colors = Blues,
          k_row = 5,
          margins = c(NA,200,60,NA),
          fontsize_row = 4,
          fontsize_col = 5,
          main="World Happiness Score and Variables by Country, 2018 \nDataTransformation using Normalise Method",
          xlab = "World Happiness Indicators",
          ylab = "World Countries"
          )

9.4 Visual Multivariate Analysis with Parallel Coordinates Plot

Parallel coordinates plot is a data visualisation specially designed for visualising and analysing multivariate, numerical data. It is ideal for comparing multiple variables together and seeing the relationships between them. For example, the variables contribute to Happiness Index. Parallel coordinates was invented by Alfred Inselberg in the 1970s as a way to visualize high-dimensional data. This data visualisation technique is more often found in academic and scientific communities than in business and consumer data visualizations. As pointed out by Stephen Few(2006), “This certainly isn’t a chart that you would present to the board of directors or place on your Web site for the general public. In fact, the strength of parallel coordinates isn’t in their ability to communicate some truth in the data to others, but rather in their ability to bring meaningful multivariate patterns and comparisons to light when used interactively for analysis.” For example, parallel coordinates plot can be used to characterise clusters detected during customer segmentation.

By the end of this hands-on exercise, you will gain hands-on experience on:

  • plotting statistic parallel coordinates plots by using ggparcoord() of GGally package,

  • plotting interactive parallel coordinates plots by using parcoords package, and

  • plotting interactive parallel coordinates plots by using parallelPlot package.

pacman::p_load(GGally, parallelPlot, tidyverse)
wh <- read_csv("data/WHData-2018.csv")

9.4.1 Plotting Static Parallel Coordinates Plot

Code chunk below shows a typical syntax used to plot a basic static parallel coordinates plot by using ggparcoord().

ggparcoord(data = wh, 
           columns = c(7:12))

9.4.2 Plotting a parallel coordinates with boxplot

The basic parallel coordinates failed to reveal any meaning understanding of the World Happiness measures. In this section, you will learn how to makeover the plot by using a collection of arguments provided by ggparcoord().

ggparcoord(data = wh, 
           columns = c(7:12), 
           groupColumn = 2,
           scale = "uniminmax",
           alphaLines = 0.2,
           boxplot = TRUE, 
           title = "Parallel Coordinates Plot of World Happines Variables")

Things to learn from the code chunk above.

  • groupColumn argument is used to group the observations (i.e. parallel lines) by using a single variable (i.e. Region) and colour the parallel coordinates lines by region name.

  • scale argument is used to scale the variables in the parallel coordinate plot by using uniminmax method. The method univariately scale each variable so the minimum of the variable is zero and the maximum is one.

  • alphaLines argument is used to reduce the intensity of the line colour to 0.2. The permissible value range is between 0 to 1.

  • boxplot argument is used to turn on the boxplot by using logical TRUE. The default is FALSE.

  • title argument is used to provide the parallel coordinates plot a title.

9.4.3 Parallel coordinates with facet

Since ggparcoord() is developed by extending ggplot2 package, we can combination use some of the ggplot2 function when plotting a parallel coordinates plot.

library(GGally)

In the code chunk below, facet_wrap() of ggplot2 is used to plot 10 small multiple parallel coordinates plots. Each plot represent one geographical region such as East Asia.

# Check the structure of your dataset
str(wh)
spc_tbl_ [156 × 12] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ Country                     : chr [1:156] "Albania" "Bosnia and Herzegovina" "Bulgaria" "Croatia" ...
 $ Region                      : chr [1:156] "Central and Eastern Europe" "Central and Eastern Europe" "Central and Eastern Europe" "Central and Eastern Europe" ...
 $ Happiness score             : num [1:156] 4.59 5.13 4.93 5.32 6.71 ...
 $ Whisker-high                : num [1:156] 4.7 5.22 5.02 5.4 6.78 ...
 $ Whisker-low                 : num [1:156] 4.48 5.04 4.84 5.24 6.64 ...
 $ Dystopia                    : num [1:156] 1.46 1.88 1.22 1.77 2.49 ...
 $ GDP per capita              : num [1:156] 0.916 0.915 1.054 1.115 1.233 ...
 $ Social support              : num [1:156] 0.817 1.078 1.515 1.161 1.489 ...
 $ Healthy life expectancy     : num [1:156] 0.79 0.758 0.712 0.737 0.854 0.737 0.732 0.578 0.671 0.716 ...
 $ Freedom to make life choices: num [1:156] 0.419 0.28 0.359 0.38 0.543 0.553 0.259 0.448 0.363 0.35 ...
 $ Generosity                  : num [1:156] 0.149 0.216 0.064 0.12 0.064 0.086 0.061 0.274 0.092 0.026 ...
 $ Perceptions of corruption   : num [1:156] 0.032 0 0.009 0.039 0.034 0.174 0.022 0.023 0.066 0.006 ...
 - attr(*, "spec")=
  .. cols(
  ..   Country = col_character(),
  ..   Region = col_character(),
  ..   `Happiness score` = col_double(),
  ..   `Whisker-high` = col_double(),
  ..   `Whisker-low` = col_double(),
  ..   Dystopia = col_double(),
  ..   `GDP per capita` = col_double(),
  ..   `Social support` = col_double(),
  ..   `Healthy life expectancy` = col_double(),
  ..   `Freedom to make life choices` = col_double(),
  ..   Generosity = col_double(),
  ..   `Perceptions of corruption` = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 
# Display the first few rows of your dataset
head(wh)
# A tibble: 6 × 12
  Country         Region `Happiness score` `Whisker-high` `Whisker-low` Dystopia
  <chr>           <chr>              <dbl>          <dbl>         <dbl>    <dbl>
1 Albania         Centr…              4.59           4.70          4.48     1.46
2 Bosnia and Her… Centr…              5.13           5.22          5.04     1.88
3 Bulgaria        Centr…              4.93           5.02          4.84     1.22
4 Croatia         Centr…              5.32           5.40          5.24     1.77
5 Czech Republic  Centr…              6.71           6.78          6.64     2.49
6 Estonia         Centr…              5.74           5.82          5.66     1.46
# ℹ 6 more variables: `GDP per capita` <dbl>, `Social support` <dbl>,
#   `Healthy life expectancy` <dbl>, `Freedom to make life choices` <dbl>,
#   Generosity <dbl>, `Perceptions of corruption` <dbl>
ggparcoord(data = wh, 
           columns = 7:12, 
           groupColumn = 2,
           scale = "uniminmax",
           alphaLines = 0.2,
           boxplot = TRUE, 
           title = "Multiple Parallel Coordinates Plots of World Happiness Variables by Region") +
  facet_wrap(~ Region)

9.4.4 Rotating x-axis text label

To make the x-axis text label easy to read, let us rotate the labels by 30 degrees. We can rotate axis text labels using theme() function in ggplot2 as shown in the code chunk below

ggparcoord(data = wh, 
           columns = c(7:12), 
           groupColumn = 2,
           scale = "uniminmax",
           alphaLines = 0.2,
           boxplot = TRUE, 
           title = "Multiple Parallel Coordinates Plots of World Happines Variables by Region") +
  facet_wrap(~ Region) + 
  theme(axis.text.x = element_text(angle = 30))

Thing to learn from the code chunk above:

To rotate x-axis text labels, we use axis.text.x as argument to theme() function. And we specify element_text(angle = 30) to rotate the x-axis text by an angle 30 degree.

9.4.5 Adjusting the rotated x-axis text label

Rotating x-axis text labels to 30 degrees makes the label overlap with the plot and we can avoid this by adjusting the text location using hjust argument to theme’s text element with element_text(). We use axis.text.x as we want to change the look of x-axis text.

ggparcoord(data = wh, 
           columns = c(7:12), 
           groupColumn = 2,
           scale = "uniminmax",
           alphaLines = 0.2,
           boxplot = TRUE, 
           title = "Multiple Parallel Coordinates Plots of World Happines Variables by Region") +
  facet_wrap(~ Region) + 
  theme(axis.text.x = element_text(angle = 30, hjust=1))

9.4.6 Plotting Interactive Parallel Coordinates Plot: parallelPlot methods

parallelPlot is an R package specially designed to plot a parallel coordinates plot by using ‘htmlwidgets’ package and d3.js. In this section, you will learn how to use functions provided in parallelPlot package to build interactive parallel coordinates plot.

The code chunk below plot an interactive parallel coordinates plot by using parallelPlot().

wh <- wh %>%
  select("Happiness score", c(7:12))
parallelPlot(wh,
             width = 320,
             height = 250)

9.4.7 Rotate axis label

In the code chunk below, rotateTitle argument is used to avoid overlapping axis labels.

parallelPlot(wh,
             rotateTitle = TRUE)

One of the useful interactive feature of parallelPlot is we can click on a variable of interest, for example Happiness score, the monotonous blue colour (default) will change a blues with different intensity colour scheme will be used.

9.4.8 Changing the colour scheme

We can change the default blue colour scheme by using continousCS argument as shown in the code chunl below.

parallelPlot(wh,
             continuousCS = "YlOrRd",
             rotateTitle = TRUE)

9.4.9 Parallel coordinates plot with histogram

In the code chunk below, histoVisibility argument is used to plot histogram along the axis of each variables.

histoVisibility <- rep(TRUE, ncol(wh))
parallelPlot(wh,
             rotateTitle = TRUE,
             histoVisibility = histoVisibility)