Exploratory Analysis of Provenance Data Using R and the Provenance Package

Abstract

The provenance of siliclastic sediment may be traced using a wide variety of chemical, mineralogical and isotopic proxies. These define three distinct data types: (1) compositional data such as chemical concentrations; (2) point-counting data such as heavy mineral compositions; and (3) distributional data such as zircon U-Pb age spectra. Each of these three data types requires separate statistical treatment. Central to any such treatment is the ability to quantify the ‘dissimilarity’ between two samples. For compositional data, this is best done using a logratio distance. Point-counting data may be compared using the chi-square distance, which deals better with missing components (zero values) than the logratio distance does. Finally, distributional data can be compared using the Kolmogorov–Smirnov and related statistics. For small datasets using a single provenance proxy, data interpretation can sometimes be done by visual inspection of ternary diagrams or age spectra. However, this no longer works for larger and more complex datasets. This paper reviews a number of multivariate ordination techniques to aid the interpretation of such studies. Multidimensional Scaling (MDS) is a generally applicable method that displays the salient dissimilarities and differences between multiple samples as a configuration of points in which similar samples plot close together and dissimilar samples plot far apart. For compositional data, classical MDS analysis of logratio data is shown to be equivalent to Principal Component Analysis (PCA). The resulting MDS configurations can be augmented with compositional information as biplots. For point-counting data, classical MDS analysis of chi-square distances is shown to be equivalent to Correspondence Analysis (CA). This technique also produces biplots. Thus, MDS provides a common platform to visualise and interpret all types of provenance data. Generalising the method to three-way dissimilarity tables provides an opportunity to combine several datasets together and thereby facilitate the interpretation of ‘Big Data’. This paper presents a set of tutorials using the statistical programming language R. It illustrates the theoretical underpinnings of compositional data analysis, PCA, MDS and other concepts using toy examples, before applying these methods to real datasets with the provenance package.

1. Introduction

At its most basic level, sedimentary provenance analysis identifies the mineralogical, chemical or isotopic composition of individual grains, or assemblages of multiple grains in siliclastic sediment. These properties can then be used to group samples of similar affinity, and thereby trace the flow of sediment through a sediment routing system, e.g., [1,2,3,4,5]. Different levels of statistical complexity arise when multiple samples are compared to each other, or when multiple provenance proxies are applied to multiple samples.

Using a number of short tutorials, this paper will introduce several simple but effective exploratory data analysis techniques that can help to make geological sense of ‘Big Data’ in a sedimentary provenance context. The term ‘exploratory’ means that these techniques allow the user to explore the data independent of any prior knowledge about the geological setting [6,7,8,9]. It groups a number of graphical methods to visualise the data and reveal patterns of similarity and differences between samples and variables. This paper will not introduce methods such as discriminant analysis that formally assign samples to pre-defined provenance areas or petrotectonic settings [10,11].

These notes do by no means claim to give a comprehensive overview of exploratory data analysis. The selection of methods presented herein is heavily biased towards techniques that are implemented in a software package for (sedimentary) geology that was created by Vermeesch et al. [12].

provenance is available free of charge at the Comprehensive R Archive Network (CRAN, https://cran.r-project.org/package=provenance), on GitHub (http://github.com/pvermees/provenance), or via http://provenance.london-geochron.com. The package is written in the statistical programming language R, which is available for Windows, Mac OS-X and Linux/Unix. The easiest way to install the latest stable version of the package is to first install R from http://r-project.org and then type the following code at the command prompt (i.e., the ‘>’):

Once installed, the package can be loaded by typing:

There are two ways to use provenance. The first of these is through a query-based user interface. To access this interface, type:

The main advantage of the query-based user interface is that it does not require any knowledge of R. Its main disadvantage is the relative lack of flexibility and the difficulty to automate complex and/or repetitive tasks. The second way to use provenance is via the R language itself. This is the quicker and more flexible option, whose only downside is a steeper learning curve compared to the query-based interface. This tutorial will help the reader to climb this learning curve whilst explaining the theoretical underpinnings of the methods that are implemented in the package.

This text assumes that the reader has a basic understanding of the R programming language, although a short tutorial is provided in the Appendix A for readers who lack such prior knowledge. The paper also assumes that the reader has some basic statistical knowledge. More specifically, (s)he is expected to be familiar with the normal distribution, and understand the meaning of the arithmetic mean, standard deviation and confidence intervals. The normal distribution underpins much of ‘conventional’ statistics, but we will see that it rarely applies to provenance data. This, in fact, is the main take-home message of this paper.

There exist three fundamental types of provenance data:

• Chemical data such as major and trace element concentrations are known as compositional data. Section 2 and Section 3 show that the statistical analysis of this class of data is fraught with difficulties. Fortunately, these are easily overcome by means of ‘Aitchison’s logratio transformation’. This transformation is a prerequisite to further statistical treatment, including Principal Component Analysis and compositional biplots of multi-sample datasets (Section 2 and Section 3).
• Categorical data such as bulk petrography and heavy mineral compositions are known as point-counting data. These are closely related to, but are fundamentally different from, compositional data. Compositional data consist of strictly positive real numbers that are subject to a constant-sum constraint and whose analytical precision can generally be ignored. In contrast, point-counting data contain integer values that may be greater than or equal to zero, and whose multinomial uncertainty is significant compared to the underlying compositional dispersion. Section 4 shows that both of these differences can be captured by a combination of logistic normal and multinomial statistics.
• Detrital age spectra form a third class of data that will be referred to as distributional data. Section 5 and Section 7 introduce kernel density estimation, the Kolmogorov–Smirnov statistic, and multidimensional scaling as ways to visualise, compare, and interpret distributional data.

Finally, Section 11 will consider the case where multiple compositional, point-counting and/or distributional datasets are combined. Procrustes analysis and 3-way multidimensional scaling are statistical techniques that aim to extract geologically meaningful trends from such ‘Big Data’ [13].

2. Ratio Data

Summary: 

This tutorial investigates the ratios of two sets of random numbers. It shows that the arithmetic mean and confidence intervals of these synthetic data yield nonsensical results. These problems are solved by a logarithmic transformation. This simple example has important implications because ratio data are common in sedimentary provenance analysis, and are closely related to compositional data, which are introduced in Section 3.

Many statistical operations assume normality. This includes averaging, the construction of confidence intervals, regression, etc. Although Gaussian distributions are common, it would be unwise to assume normality for all datasets. This paper makes the point that, more often than not, the normality assumption is invalid in the context of sedimentary provenance analysis. Ignoring this non-normality can lead to counter-intuitive and plainly wrong results.

To illustrate this point, we will now consider the simple case of ratio data, which are quite common in the Earth Sciences. Take, for example, the ratio of apatite to tourmaline in heavy mineral analysis, which has been used to indicate the duration of transport and storage prior to deposition [14]. In this part of the tutorial, we will investigate the statistics of ratio data using a synthetic example.

• Create two vectors A and B, each containing 100 random numbers between 0 and 1: Intuitively, given that $A/B=1/(B/A)$ and $B/A=1/(A/B)$, we would expect the same to be true for their means $\overline{(A/B)}$ and $\overline{(B/A)}$. However, when we define two new variables for the (inverse) of the (reciprocal) mean ratios: then we find that AB.mean≠inv.BA.mean. So $\overline{(A/B)}\ne 1/\overline{(B/A)}$ and $\overline{(B/A)}\ne 1/\overline{(A/B)}$! This is a counterintuitive and clearly wrong result.
• Calculate the standard deviation of $A/B$ and multiply this by two to obtain a ‘2-sigma’ confidence interval for the data: then we find that $LL<0$, which is nonsensical since A and B are both strictly positive numbers and their ratio is therefore not allowed to take negative values either. Herein lies the root of the problem. The sampling distribution of A/B is positively skewed, whereas the normal distribution is symmetric with tails ranging from $-\infty$ to $+\infty$. Geologists frequently encounter strictly positive numbers. Time, for example, is a strictly positive quantity, expressed by geochronologists as ‘years before present’, where ‘present’ is equivalent to zero.
• The problems caused by applying normal theory to strictly positive data can often be solved by simply taking logarithms [15]. The transformed data are then free to take on any value, including negative values, and this often allows normal theory to be applied with no problems. For example, when we calculate the (geometric) mean after taking the logarithm of the ratio data: then we find that AB.gmean = inv.BA.gmean, which is a far more sensible result.
• Calculating the 2-sigma interval for the log-transformed data: also produces strictly positive values, as expected.

3. Compositional Data

Summary: 

Compositional data such as chemical concentrations suffer from the same problems as the ratio data of Section 2. The tutorial uses a geochemical dataset of Al₂O₃ – (CaO+Na₂O) – K₂O data to demonstrate that the ‘conventional’ arithmetic mean and confidence intervals are inappropriate for data that can be constrained to a constant sum. A logratio transformation solves these problems.

Like the ratios of the previous Section, the chemical compositions of rocks and minerals are also expressed as strictly positive numbers. They, however, do not span the entire range of positive values, but are restricted to a narrow subset of that space, ranging from 0 to 1 (if fractions are used) or from 0 to 100% (using percentage notation). The compositions are further restricted by a constant sum constraint:

$$\sum _{i=1}^n{C}_i=1$$

for an n-component system. Consider, for example, a three-component system $\{x,y,z\}$, where $x+y+z=1$. Such compositions can be plotted on ternary diagrams, which are very popular in geology. Well-known examples are the Q-F-L diagram of sedimentary petrography [16], the A-CN-K diagram in weathering studies [17], and the A-F-M, Q-A-P and Q-P-F diagrams of igneous petrology [18]. The very fact that it is possible to plot a ternary diagram on a two-dimensional sheet of paper already tells us that it really displays only two and not three dimensions worth of information. Treating the ternary data space as a regular Euclidean space with Gaussian statistics leads to incorrect results, as illustrated by the following example.

• Read a compositional dataset containing the Al₂O₃ – (CaO+Na₂O) – K₂O composition of a number of synthetic samples: where row.names=1 indicates that the sample names are contained in the first column; and the header=TRUE and check.names=FALSE arguments indicate that the first column of the input table contains the column headers, one of which contains a special character (‘+’).
• Calculate the arithmetic mean composition and 95% confidence limits for each column of the dataset: and construct the 2-sigma confidence confidence bounds:
• In order to plot the compositional data on a ternary diagram, we will need to first load the provenance package into memory: Now plot the Al₂O₃, (CaO + Na₂O) and K₂O compositions on a ternary diagram alongside the arithmetic mean composition: where ternary(x) creates a ternary data ‘object’ from a variable x, and pch = 20 and pch = 22 produce filled circles and squares, respectively. Notice how the arithmetic mean plots outside the data cloud, and therefore fails to represent the compositional dataset (Figure 1).
• Add a 2-sigma confidence polygon to this figure using the ternary.polygon() function that is provided in the auxiliary helper.R script (see Online Supplement): Note that the polygon partly plots outside the ternary diagram, into physically impossible negative data space. This nonsensical result is diagnostic of the dangers of applying ‘normal’ statistics to compositional data. It is similar to the negative limits for the ratio data in Section 2.

A comprehensive solution to the compositional data conundrum was only found in the 1980s, by Scottish statistician John Aitchison [19]. It is closely related to the solution of the ratio averaging problem discussed in the previous section. The trick is to map the n-dimensional composition to an (n-1)-dimensional Euclidean space by means of a logratio transformation. For example, in the ternary case, we can map the compositional variables x, y and z to two transformed variables v and w:

$$v=ln\left(\frac{x}{z}\right),w=ln\left(\frac{y}{z}\right)$$

(1)

After performing the statistical analysis of interest (e.g., calculating the mean or constructing a 95% confidence region) on the transformed data, the results can then be mapped back to compositional space with the inverse logratio transformation. For the ternary case:

$$x=\frac{e^v}{e^v+e^w+1},y=\frac{e^w}{e^v+e^w+1},z=\frac{1}{e^v+e^w+1}$$

(2)

This transformation is implemented in the provenance package. Let us use this feature to revisit the K-CN-A dataset, and add the geometric mean and 95% confidence region to the ternary diagram for comparison with the arithmetic mean and confidence polygon obtained before.

5.

Compute the logratio mean composition and add it to the existing ternary diagram as a red square: This red square falls right inside the data cloud, an altogether more satisfying result than the arithmetic mean shown in blue (Figure 1).

6.

To add a compositional confidence contour, we must re-read ACNK.csv into memory using the read.compositional() function. This will tell the provenance package to treat the resulting variable as compositional data in subsequent operations: Adding the 95% confidence contour using provenance’s ternary.ellipse() function: creates a 95% confidence ellipse in logratio space, and maps this back to the ternary diagram. This results in a ‘boomerang’-shaped contour that tightly hugs the compositional data whilst staying inside the boundaries of the ternary diagram (Figure 1).

This Section (and Section 8) only touched the bare essentials of compositional data analysis. Further information about this active field of research can be found in Pawlowsky-Glahn et al. [20]. For additional R-recipes for compositional data analysis using the compositions package, the reader is referred to Van den Boogaart and Tolosana-Delgado [21,22].

4. Point-Counting Data

Summary: 

Point-counting data such as heavy mineral counts are underlain by compositional distributions. However, they are not amenable to the logratio transformations introduced in Section 3 because they commonly contain zero values. Averages and confidence intervals for this type of data require hybrid statistical models combining compositional and multinomial aspects.

The mineralogical composition of silicilastic sediment can be determined by tallying the occurrence of various minerals in a representative sample of (200–400, say) grains [23,24]. Such point-counting data are closely related to the compositional data that were discussed in the previous section. However, there are some crucial differences between these two data classes [25].

Point-counting data are associated with significant (counting) uncertainties, which are ignored by classical compositional data analysis. As a consequence, point-counting data often contain zero values, which are incompatible with the log-ratio transformation defined in Equation (1). Although ‘rounding zeros’ also occur in compositional data, where they can be removed by ‘imputation’ methods [26,27], these methods are ill-suited for point-counting datasets in which zeros are the rule rather than the exception.

• Download the auxiliary data file HM.csv from the Online Supplement. This file contains a heavy mineral dataset from the Namib Sand Sea [13]. It consists of 16 rows (one for each sample) and 15 columns (one for each mineral). Read these data into memory and tell provenance to treat it as point-counting data in all future operations: Galbraith [28]’s radial plot is an effective way to visually assess the degree to which the random counting uncertainties account for the observed scatter of binary point-counting data. Applying this to the epidote/garnet-ratio of the heavy mineral data (Figure 2): Each circle on the resulting scatter plot represents a single sample in the HM dataset. Its epidote/garnet-ratio can be obtained by projecting the circle onto the circular scale. Thus, low and high ratios are found at negative and positive angles to the origin, respectively. The horizontal distance of each point from the origin is proportional to the total number of counts in each sample and, hence, to its precision. An (asymmetric) 95% confidence interval for the ep/gt-ratio of each sample can be obtained by projecting both ends of a 2-sigma confidence bar onto the circular scale.

Suppose that the data are underlain by a single true population and random counting uncertainties are the sole source of scatter. Let $\theta$ be the true but unknown proportion of the binary subpopulation that consists of the first mineral (epidote, say). Then $(1-\theta )$ is the fraction of grains that belong to the second mineral (garnet). Further suppose that we have counted a representative sample of N grains from this population. Then the probability that this sample contains n grains of the first mineral and $m=N-n$ grains of the second mineral follows a binomial distribution:

$$p\left(n\right)=\left(\genfrac{}{}{0pt}{}{n+m}{n}\right){\theta }^n{(1-\theta )}^m$$

(3)

If multiple samples in a dataset are indeed underlain by the same fraction $\theta$, then approximately 95% of the samples should fit within a horizontal band of two standard errors drawn on either side of the origin. In this case, $\theta$ can be estimated by pooling all the counts together and computing the proportion of the first mineral as a fraction of the total number of grains counted [25].

However, the ep/gt-ratios in HM scatter significantly beyond the 2-sigma band (Figure 2i). The data are therefore said to be overdispersed with respect to the counting uncertainties. This indicates the presence of true geological dispersion in the compositions that underly the point-counting data. The dispersion can be estimated by a random effects model with two parameters:

$$\beta \equiv ln\left(\frac{\theta }{1-\theta }\right)\approx \mathcal{N}(\mu ,{\sigma }^2)$$

(4)

where $\beta$ is a new variable that follows a normal distribution with mean $\mu$ and standard deviation $\sigma$, both of which have geological significance.

The ‘central ratio’ is given by $exp\left[\widehat{\mu }\right]$ where $\widehat{\mu }$ is the maximum likelihood estimate for $\mu$. This estimates the geometric mean (ep/gt-) ratio of the true underlying composition. The ‘dispersion’ ($\widehat{\sigma }$) estimates the geological scatter [25,29]. In the case of our heavy mineral dataset, the epidote-garnet subcomposition is 75% dispersed. This means that the coefficient of variation (geometric standard deviation divided by geometric mean) of the true epidote/garnet-ratios is approximately 0.75.

2.

The continuous mixtures from the previous section can be generalised from two to three or more dimensions. The following code snippet uses it to construct a 95% confidence contour for the ternary subcomposition of garnet, epidote and zircon (Figure 2ii). Note that this dataset contains four zero values, which would have rendered the logratio approach of Figure 1 unusable.

3.

For datasets comprising more than three variables, the central composition can be simply obtained as follows: This produces a matrix with the proportions of each component; its standard error; the dispersion of the binary subcomposition formed by the component and the amalgamation of all remaining components; and the outcome of a chi-square test for homogeneity.

5. Distributional Data

Summary: 

This Section investigates a 16-sample, 1547-grain dataset of detrital zircon U-Pb ages from Namibia. It uses Kernel Density Estimation and Cumulative Age Distributions to visualise this dataset, and introduces the Kolmogorov–Smirnov statistic as a means of quantifying the dissimilarity between samples.

Compositional data such as the chemical concentrations of Section 3 and Section 8 are characterised by the relative proportions of a number of discrete categories. A second class of provenance proxies is based on the sampling distribution of continuous variables such as zircon U-Pb ages [30,31]. These distributional data do not fit in the statistical framework of the (logistic) normal distribution.

• Download auxiliary data file DZ.csv from the Online Supplement. This file contains a detrital zircon U-Pb dataset from Namibia. It consists of 16 columns—one for each sample—each containing the single grain U-Pb ages of their respective sample. Let us load this file into memory using provenance’s read.distributional() function:DZ now contains an object of class distributional containing the zircon U-Pb ages of 16 Namibian sand samples. To view the names of these samples:
• One way to visualise the U-Pb age distributions is as Kernel Density Estimates. A KDE is defined as:

  $$KD{E}_x\left(t\right)=\frac{1}{n}\sum _{i=1}^n\mathcal{K}\left(t\right|x_i,bw)$$

  (5)

  where $\mathcal{K}$ is the ‘kernel’ and $bw$ is the ‘bandwidth’ [32,33]. The kernel can be any unimodal and symmetric shape (such as a box or a triangle), but is most often taken to be Gaussian (where $x_i$ is the mean and $bw$ the standard deviation). The bandwidth can either be set manually, or selected automatically based on the number of data points and the distance between them. provenance implements the automatic bandwidth selection algorithm of Botev et al. [34] but a plethora of alternatives are available in the statistics literature. To plot all the samples as KDEs: where ncol specifies the number of columns over which the KDEs are divided.
• Alternatively, the Cumulative Age Distribution (CAD) is a second way to show the data [35]. A CAD is a step function that sets out the rank order of the dates against their numerical value:

  $$\mathrm{CAD}\left(t\right)=\sum _{i=1}^{n}1(t<t_i)/n$$

  (6)

  where $1(*)=1$ if * is true and $1(*)=0$ if * is false. The main advantages of CADs over KDEs are that (i) they do not require any smoothing (i.e., there is no ‘bandwidth’ to choose), and (ii) they can superimpose multiple samples on the same plot. Plotting samples N1, N2 and N4 of the Namib dataset: we can see that (1) the CADs of samples N1 and N2 plot close together with steepest sections at 500 Ma and 1000 Ma, reflecting the prominence of those age components; (2) sample N4 is quite different from N1 and N2.
• We can quantify this difference using the Kolmogorov–Smirnov (KS) statistic [36,37,38], which represents the maximum vertical difference between two CADs: This shows that the KS-statistic between N1 and N2 is KS(N1,N2) = 0.18, whereas KS(N1,N4) = 0.44, and KS(N2,N4) = 0.35 (Figure 3). The KS statistic is a non-negative value that takes on values between zero (perfect overlap between two distributions) and one (no overlap between two distributions). It is symmetric because the KS statistic between any sample x and another sample y equals that between y and x. For example, KS(N1,N2) = 0.18 = KS(N2,N1). Finally, the KS-statistic obeys the triangle equality, which means that the dissimilarity between any two samples is always smaller than or equal to the sum of the dissimilarities between those two samples and a third. For example, KS(N1,N2) = 0.18 < KS(N1,N4) + KS(N2,N4) = 0.44 + 0.35 = 0.79. These three characteristics qualify the KS statistics as a metric, which makes it particularly suitable for Multidimensional Scaling (MDS) analysis (see Section 7). The KS statistic is just one of many dissimilarity measures for distributional data. However, not all these alternatives to the KS statistic fulfil the triangle inequality [38].

6. Principal Component Analysis (PCA)

Summary: 

Principal Component Analysis is an exploratory data analysis method that takes a high dimensional dataset as input and produces a lower (typically two-) dimensional ‘projection’ as output. PCA is closely related to Multidimensional Scaling (MDS), compositional PCA, and Correspondence Analysis (CA), which are introduced in Section 7, Section 8 and Section 9. This tutorial introduces PCA using the simplest working example of three two-dimensional points. Nearly identical examples will be used in Section 7, Section 8 and Section 9.

• Consider the following bivariate (a and b) dataset of three (1, 2 and 3) samples:

  $$X=\begin{array}{ccc}\begin{array}{c}\\ 1\\ 2\\ 3\end{array}& \begin{array}{c}a\\ \left[\begin{array}{c}\hfill \mathrm{-1}\\ \hfill 3\\ \hfill 4\end{array}\right.\end{array}& \begin{array}{c}b\\ \left.\begin{array}{c}\hfill 7\\ \hfill 2\\ \hfill 3\end{array}\right]\end{array}\end{array}$$

  (7)

  Generating and plotting X in R: yields a diagram in which two of the three data points plot close together while the third one plots further away.
• Imagine that you live in a one-dimensional world and cannot see the spatial distribution of the three points represented by X. Principal Component Analysis (PCA) is a statistical technique invented by Pearson [39] to represent multi- (e.g., two-) dimensional data in a lower- (e.g., one-) dimensional space whilst preserving the maximum amount of information (i.e., variance). This can be achieved by decomposing X into four matrices (C, S, V and D):

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill X=& 1_{\mathrm{3,1}}C+SVD\hfill \\ \hfill =& \left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]\left[\begin{array}{cc}2& 4\end{array}\right]+\left[\begin{array}{cc}-1.15& 0\\ 0.58& -1\\ 0.58& 1\end{array}\right]\left[\begin{array}{cc}3.67& 0\\ 0& 0.71\end{array}\right]\left[\begin{array}{cc}0.71& -0.71\\ 0.71& 0.71\end{array}\right]\hfill \end{array}\begin{array}{c}\hfill \end{array}\end{array}$$

  (8)

  where C is the centre (arithmetic mean) of the two data columns; S are the normalised scores; the diagonals of V correspond to the standard deviations of the two principal components; and D is a rotation matrix (the principal directions). S, V and D can be recombined to define two more matrices:

  $$\begin{array}{cc}\hfill P& =SV=\left[\begin{array}{cc}-4.24& 0\\ 2.12& -0.71\\ 2.12& 0.71\end{array}\right],\hfill \end{array}$$

  (9)

  $$\begin{array}{cc}\hfill \mathrm{and}L& =VD=\left[\begin{array}{cc}2.6& -2.6\\ 0.5& 0.5\end{array}\right]\hfill \end{array}$$

  (10)

  where P is a matrix of transformed coordinates (the principal components or scores) and L contains the scaled eigenvectors or loadings. Figure 4i shows X as numbers on a scatterplot, C as a yellow square, and $1_{\mathrm{2,1}}C\pm L$ as a cross. Thus, the first principal direction (running from the upper left to the lower right) has been stretched by a factor of $(3.67/0.71)=5.2$ with respect to the second principal component, which runs perpendicular to it. The two principal components are shown separately as Figure 4ii, and their relative contribution to the total variance of the data as Figure 4iv. Figure 4 can be reproduced with the following R-code:
• Although the two-dimensional example is useful for illustrative purposes, the true value of PCA obviously lies in higher dimensional situations. As a second example, let us consider one of R’s built-in datasets. USArrests contains statistics (in arrests per 100,000 residents) for assault, murder, and rape in each of the 50 US states in 1973. Also given is the percentage of the population living in urban areas. Thus, USArrests is a four-column table that cannot readily be visualised on a two-dimensional surface. Applying PCA yields four principal components, the first two of which represent 62% and 25% of the total variance, respectively. Because the four columns of the input data are expressed in different units (arrests per 100,000 or percentage), it is necessary to scale the data to have unit variance before the analysis takes place: The resulting biplot shows that the loading vectors for Murder, Assault and Rape are all pointing in approximately the same direction (dominating the first principal component), perpendicular to UrbanPop (which dominates the second principal component). This tells us that crime and degree of urbanisation are not correlated in the United States.

7. Multidimensional Scaling

Summary: 

Multidimensional Scaling (MDS) is a less restrictive superset of PCA. This tutorial uses a geographical example to demonstrate how MDS re-creates a map of Europe from a table of pairwise distances between European cities. Applying the same algorithm to the synthetic toy-example of Section 6 yields exactly the same output as PCA.

• Multidimensional Scaling (MDS [40,41,42,43]) is a dimension-reducing technique that aims to extract two- (or higher) dimensional ‘maps’ from tables of pairwise distances between objects. This method is most easily illustrated with a geographical example. Consider, for example, the eurodist dataset that is built into R, and which gives the road distances (in km) between 21 cities in Europe (see ?eurodist for further details):
• The MDS configuration can be obtained by R’s built-in cmdscale() function Set up an empty plot with a 1:1 aspect ratio, and then label the MDS configuration with the city names: Note that the map may be turned ‘upside down’. This reflects the rotation invariance of MDS configurations.
• R’s cmdscale() function implements so-called ‘classical’ MDS, which aims to fit the actual distances [42,43]. If these distances are Euclidean, then it can be shown that MDS is equivalent to PCA [44,45,46]. To demonstrate this equivalence, let us apply MDS to the data in Equation (7). First, run the first two lines of code from part 1 in Section 6. Calculating the Euclidean distances between the three samples produces a dissimilarity matrix d. For example, the distance between points 1 and 2 is $\sqrt{{(-1-3)}^2+{(7-2)}^2}=6.4$. This value is stored in d[1,2]. In R: which produces:

  $$d=\begin{array}{cccc}\begin{array}{c}\\ 1\\ 2\\ 3\end{array}& \begin{array}{c}1\\ \left[\begin{array}{c}\hfill 0\\ \hfill 6.4\\ \hfill 6.4\end{array}\right.\end{array}& \begin{array}{c}2\\ 6.4\\ 0\\ 1.4\end{array}& \begin{array}{c}3\\ \left.\begin{array}{c}\hfill 6.4\\ \hfill 1.4\\ \hfill 0\end{array}\right]\end{array}\end{array}$$

  (11)
• Next, calculate the MDS configuration: Finally, plot the MDS configuration as a scatterplot of text labels: which is identical to the PCA configuration of Figure 4iii apart from an arbitrary rotation or reflection.
• An alternative implementation of MDS loosens the Euclidean distance assumption by fitting the relative distances between objects [40,41]. Let us apply this to the dataset of European city distances using the isoMDS function of the ‘Modern Applied Statistics with S’ (MASS) package [47]: To compute and plot the non-metric MDS configuration: where conf3 is a list with two items: stress, which expresses the goodness-of-fit of the MDS configuration; and points, which contains the configuration. The ‘$’ operator is used to access any of these items.
  Non-metric MDS is a less-restrictive superset of classical MDS and, hence, PCA, which opens this methodology up to non-Euclidean dissimilarity measures, such as the KS-distance introduced in Section 5 [48].

8. PCA of Compositional Data

Summary: 

PCA can be applied to compositional data after logratio transformation. This tutorial first applies such compositional PCA to a three sample, three variable dataset that is mathematically equivalent to the three sample two variable toy example of Section 6. Then, it applies the same method to a real dataset of major element compositions from Namibia. This is first done using basic R and then again (and more succinctly) using the provenance package.

Consider the following trivariate (a, b and c) dataset of three (1, 2 and 3) compositions that are constrained to a constant sum (a_i + b_i + c_i = 100% for 1 ≤ i ≤ 3, Figure 5):

$$X=\begin{array}{cccc}\begin{array}{c}\\ 1\\ 2\\ 3\end{array}& \begin{array}{c}\mathrm{a}\\ \left[\begin{array}{c}\hfill 0.034\\ \hfill 69.45\\ \hfill 72.44\end{array}\right.\end{array}& \begin{array}{c}\mathrm{b}\\ 99.88\\ 25.55\\ 26.65\end{array}& \begin{array}{c}\mathrm{c}\\ \left.\begin{array}{c}\hfill 0.091\\ \hfill 5.01\\ \hfill 0.92\end{array}\right]\end{array}\end{array}$$

(12)

It would be wrong to apply conventional PCA to this dataset, because this would ignore the constant sum constraint. As was discussed in Section 6, PCA begins by ‘centering’ the data via the arithmetic mean. Section 3 showed that this yields incorrect results for compositional data. Although the additive logratio transformation (alr) of Equation (1) solves the closure problem, it is not suitable for PCA because it is not isometric. For example, the alr-distance between samples 2 and 3 is 1.74 if b is used as a common denominator, but 2.46 if c is used as a common denominator.

The fact that distances are not unequivocally defined in alr-space spells trouble for PCA. Recall the equivalence of PCA and classical MDS, which was discussed in Section 7. MDS is based on dissimilarity matrices, so if distances are not well defined then neither are the MDS configuration and, hence, the principal components. This issue can be solved by the centred logratio transformation (clr):

$$u_i=\ln \left[\frac{x_i}{g_i}\right], v_i=\ln \left[\frac{y_i}{g_i}\right], \mathrm{and} w_i=\ln \left[\frac{z_i}{g_i}\right]$$

(13)

where g_i is the geometric mean of the ith sample:

$$g_i=\exp \left[\frac{\ln \left[x_i\right]+\ln \left[y_i\right]+\ln \left[z_i\right]}{3}\right]$$

Applying the clr-transformation to the data of Equation (12) yields a new trivariate dataset:

$$Xc=\begin{array}{cccc}\begin{array}{c}\\ 1\\ 2\\ 3\end{array}& \begin{array}{c}ln(a/g)\\ \left[\begin{array}{c}\hfill -3\\ \hfill 1.21\\ \hfill 1.79\end{array}\right.\end{array}& \begin{array}{c}ln(b/g)\\ 5\\ 0.21\\ 0.79\end{array}& \begin{array}{c}ln(c/g)\\ \left.\begin{array}{c}\hfill -2\\ \hfill -1.42\\ \hfill -2.58\end{array}\right]\end{array}\end{array}$$

(14)

where g stands for the geometric mean of each row. Note that each of the rows of X_c adds up to zero. Thanks to the symmetry of the clr-coordinates, the distances between the rows (which are also known as Aitchison distances) are well defined. Subjecting Equation (14) to the same matrix decomposition as Equation (8) yields:

$$X_c=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]\left[\begin{array}{ccc}0& 2& -2\end{array}\right]+\left[\begin{array}{ccc}-1.15& 0& 0\\ 0.58& -1& 0\\ 0.58& 1& 0\end{array}\right]\left[\begin{array}{ccc}3.67& 0& 0\\ 0& 0.71& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}0.71& -0.71& 0\\ 0.71& 0.71& -0.82\\ -0.58& -0.58& -0.58\end{array}\right]$$

(15)

so that

$$P=\left[\begin{array}{ccc}-4.24& 0& 0\\ 2.12& -0.71& 0\\ 2.12& 0.71& 0\end{array}\right]\mathrm{and}L=\left[\begin{array}{ccc}2.59& -2.59& 0\\ 0.29& 0.29& -0.58\\ 0& 0& 0\end{array}\right]$$

(16)

Note that, even though this yields three principal components instead two, the variance of the third component in matrix V is zero. Therefore, all the information is contained in the first two components. Furthermore, note that the first two principal components of the compositional dataset are identical to those of the PCA example shown in Section 6 (Equation (9)). This is, of course, intentional.

• The following script applies compositional PCA to a dataset of major element compositions from Namibia (see Online Supplement) using base R:
• Alternatively, we can also do this more easily in provenance: where the read.compositional function reads the .csv file into an object of class compositional, thus ensuring that logratio statistics are used in all provenance functions (such as PCA) that accept compositional data as input. Also note that the provenance package overloads the plot function to generate a compositional biplot when applied to the output of the PCA function.

9. Correspondence Analysis

Summary: 

Point-counting data can be analysed by MDS using the Chi-square distance. Correspondence Analysis (CA) yields identical results whilst producing biplots akin to those obtained by PCA. This tutorial first uses a simple three sample, three variable toy example that is (almost) identical to those used in Section 6, Section 7 and Section 8, before applying CA to a real dataset of heavy mineral counts from Namibia.

Consider the following three sets of trivariate point-counting data:

$$X=\begin{array}{cccc}\begin{array}{c}1\\ 2\\ 3\end{array}& \begin{array}{c}a\\ \left[\begin{array}{c}\hfill 0\\ \hfill 38\\ \hfill 108\end{array}\right.\end{array}& \begin{array}{c}b\\ 100\\ 13\\ 38\end{array}& \begin{array}{c}c\\ \left.\begin{array}{c}\hfill 0\\ \hfill 1\\ \hfill 0\end{array}\right]\end{array}\end{array}$$

(17)

This dataset intentionally looks similar on a ternary diagram to the compositional dataset of Section 3. The only difference is the presence of zeros, which preclude the use of logratio statistics. This problem can be solved by replacing the zero values with small numbers, but this only works when their number is small [26,27]. Correspondence Analysis (CA) is an alternative approach that does not require such ‘imputation’.

CA is a dimension reduction technique that is similar in many ways to PCA [25,49]. CA, like PCA, is a special case of MDS. Whereas ordinary PCA uses the Euclidean distance, and compositional data can be compared using the Aitchison distance, point-counting data can be compared by means of a chi-square distance:

$$d_{ij}=\sqrt{\sum _{k=1}^K\frac{X_{··}}{X_{·k}}{\left(\frac{X_{ik}}{X_{i·}}-\frac{X_{jk}}{X_{j·}}\right)}^2}$$

(18)

where $X_{·k}={\sum }_{i=1}^m{X}_{ik}$, $X_{i·}={\sum }_{k=1}^K{X}_{ik}$ and $X_{··}={\sum }_{i=1}^m{\sum }_{k=1}^K{X}_{ik}$. Applying this formula to the data of Equation (17) produces the following dissimilarity matrix:

$$\begin{array}{cccc}\begin{array}{c}1\\ 2\\ 3\end{array}& \begin{array}{c}1\\ \left[\begin{array}{c}\hfill 0\\ \hfill 1.5\\ \hfill 1.5\end{array}\right.\end{array}& \begin{array}{c}2\\ 1.5\\ 0\\ 0.33\end{array}& \begin{array}{c}3\\ \left.\begin{array}{c}\hfill 1.5\\ \hfill 0.33\\ \hfill 0\end{array}\right]\end{array}\end{array}$$

(19)

Note that, although these values are different than those in Equation (11), the ratios between them are (approximately) the same. Specifically, $d_{1,2}/d_{\mathrm{1,3}}=1.5/1.5=1$ for Equation (19) and $d_{\mathrm{1,2}}/d_{\mathrm{1,3}}=6.4/6.4=1$ for Equation (11); or $d_{\mathrm{1,2}}/d_{\mathrm{2,3}}=1.5/0.33=4.5$ for Equation 19) and $d_{\mathrm{1,2}}/d_{\mathrm{2,3}}=6.4/1.4=4.5$ for Equation (11). Therefore, when we subject our point-counting data to an MDS analysis using the chi-square distance, the resulting configuration appears nearly identical to the example of Section 7.

The following script applies CA to the heavy mineral composition of Namib desert sand. It loads a table called HM.csv that contains point counts for 16 samples and 15 minerals. To reduce the dominance of the least abundant components, the code extracts the most abundant minerals (epidote, garnet, amphibole and clinopyroxene) from the datasets and amalgamates the ultra-stable minerals (zircon, tourmaline and ru- tile), which have similar petrological significance.

10. MDS Analysis of Distributional Data

Summary: 

This brief tutorial applies MDS to the detrital zircon U-Pb dataset from Namibia, using the Kolmogorov–Smirnov statistic as a dissimilarity measure.

Part 5 in Section 7 introduced non-metric MDS as a less-restrictive superset of classical MDS and, hence, PCA. This opens this methodology up to non-Euclidean dissimilarity measures, such as the KS-distance introduced in part 4 in Section 5 [38,48].

In this case, the overloaded plot function produces not one but two graphics windows. The first of these shows the MDS configuration, whereas the second shows the Shepard plot [40,41]. This is a scatterplot that sets out the Euclidean distances between the samples measured on the MDS configuration against the disparities, which are defined as:

$$\delta [i,j]=f\left(KS\right[i,j\left]\right)$$

(20)

where $KS[i,j]$ is the KS-distance between the i^th and j^th sample and f is a monotonic transformation, which is shown as a step-function. The Shepard plot allows the user to visually assess the goodness-of-fit of the MDS configuration. This can be further quantified using the ‘stress’ parameter:

$$S=\sum _i\sum _j{(d[i,j]-\delta [i,j])}^2/\sum _i\sum _j{\left(d[i,j]\right)}^2$$

(21)

The lower the stress, the better the fit. For moderately sized datasets, stress values should be less than 10% [40]. For larger datasets, a higher dimensional solution may be necessary, using the optional parameter k of provenance’s MDS function [50].

11. ‘Big’ Data

Summary: 

The tutorial jointly analyses 16 Namibian samples using five different provenance proxies, including all three data classes introduced in Section 3, Section 4 and Section 5. It introduces Procrustes Analysis and 3-way MDS as two alternative ways to extract geologically meaningful information from these multivariate ‘big’ dataset.

It is increasingly common for provenance studies to combine compositional, point-counting or distributional datasets together [4,13]. Linking together bulk sediment data, heavy mineral data and single mineral data requires not only a sensible statistical approach, but also a full appraisal of the impact of mineral fertility and heavy mineral concentration in eroded bedrock and derived clast sediment [51,52,53]. Assuming that such an appraisal has been made, this Section introduces some exploratory data analysis tools that can reveal meaningful structure in complex datasets.

• The full Namib Sand Sea study that we have used as a test case for this tutorial comprises five datasets (see Online Supplement):

  (a)

  Major element concentrations (Major.csv, compositional data)

  (b)

  Trace element concentrations (Trace.csv, compositional data)

  (c)

  Bulk petrography (PT.csv, point-counting data)

  (d)

  Heavy mineral compositions (HM.csv, point-counting data)

  (e)

  Detrital zircon U-Pb data (DZ.csv, distributional data)

  All these datasets can be visualised together in a single summary plot: where Major, Trace, QFL and HM are shown as pie charts (the latter two with a different colour map than the former), and DZ as KDEs. Adding DZ instead of KDEs(DZ) would plot the U-Pb age distributions as histograms.
• The entire Namib dataset comprises 16,125 measurements spanning five dimensions worth of compositional, distributional and point-counting information. This complex dataset, which may be rightfully described by the internet-era term of ‘Big Data’, is extremely difficult to interpret by mere visual inspection of the pie charts and KDEs. Applying MDS/PCA to each of the five individual datasets helps but presents the analyst with a multi-plot comparison problem. provenance implements two methods to address this issue [13]. The first of these is called ‘Procrustes Analysis’ [54]. Given a number of MDS configurations, this technique uses a combination of transformations (translation, rotation, scaling and reflection) to extract a ‘consensus view’ for all the data considered together:
• Alternatively, ‘3-way MDS’ is an extension of ‘ordinary’ (2-way) MDS that accepts 3-dimensional dissimilarity matrices as input. provenance includes the most common implementation of this class of algorithms, which is known as ‘INdividual Difference SCALing’ or INDSCAL [55,56]: This code produces two pieces of graphical output (Figure 6). The ‘group configuration’ represents the consensus view of all provenance proxies considered together. This looks very similar to the Procrustes configuration created by the previous code snippet. The second piece of graphical information displays not the samples but the provenance proxies. It shows the weights that each of the proxies attach to the horizontal and vertical axis of the group configuration.
  For example, the heavy mineral compositions of the Namib desert sands can be (approximately) described by stretching the group configuration vertically by a factor of 1.9, whilst shrinking it horizontally by a factor of 0.4. In contrast, the configurations of the major and trace element compositions for the same samples are obtained by shrinking the group configuration vertically by a factor 0.8, and stretching it horizontally by a factor of 1.3. Thus, by combining these weights with the group configuration yields five ‘private spaces’ that aim to fit each of the individual datasets.
  INDSCAL group configurations are not rotation-invariant, in contrast with the 2-way MDS configurations of Section 7. This gives geological meaning to the horizontal and vertical axes of the plot. For example, samples N1 and N10 plot along a vertical line on the group configuration, indicating that they have different heavy mineral compositions, but similar major and trace element compositions. On the other hand, samples N4 and N8 plot along a horizontal line, indicating that they have similar major and trace element compositions but contrasting heavy mineral compositions.
  Closer inspection of the weights reveals that the datasets obtained from fractions of specific densities (HM, PT and DZ) attach stronger weights to the vertical axis, whereas those that are determined on bulk sediment (Major and Trace) dominate the horizontal direction. Provenance proxies that use bulk sediment are more sensitive to winnowing effects than those that are based on density separates. This leads to the interpretation that the horizontal axis separates samples that have been affected by different degrees of hydraulic sorting, whereas the vertical direction separates samples that have different provenance.

12. Summary, Conclusions and Outlook

The statistical toolbox implemented by the provenance package is neither comprehensive nor at the cutting edge of exploratory data analysis. PCA, MDS, CA, and KDEs are tried and tested methods that have been around for many decades. Nothing new is presented here and that is intentional. This paper makes the point that even the most basic statistical parameters like the arithmetic mean and standard deviation cannot be blindly applied to geological data [24,57,58]. Great care must be taken when applying established techniques to sedimentary provenance data such as chemical compositions, point-counts or U-Pb age distributions. Given the difficulty of using even the simplest of methods correctly, geologists may want to think twice before exploring more complicated methods, or inventing entirely new ones.

The set of tutorials presented in this paper did not cover all aspects of statistical provenance analysis. Doing so would fill a book rather than a paper. Some additional topics for such a book could include (1) supervised and unsupervised learning algorithms such as cluster analysis and discriminant analysis, which can group samples into formal groups [10,11,59,60]; (2) the physical and chemical processes that affect the composition of sediment from ‘source to sink’ [5,61,62,63]; and (3) quality checks and corrections that must be made to ensure that the data reveal meaningful provenance trends rather than sampling effects [51,52,64,65,66].

The paper introduced three distinct classes of provenance data. Compositional, point-counting and distributional data each require different statistical treatment. Multi-sample collections of these data can be visualised by Multidimensional Scaling, using different dissimilarity measures (

Table 1. A summary of the three types of provenance data introduced in this paper along with a suitable dissimilarity measure and its corresponding ordination technique.

Data Type	Dissimilarity Measure	Ordination Technique
compositional	Aitchison	Principal Component Analysis
point-counting	Chi-square	Correspondence Analysis
distributional	Kolmogorov–Smirnov	Multidimensional Scaling

). Distributional data can be compared using the Kolmogorov–Smirnov statistic or related dissimilarity measures, and plugged straight into an MDS algorithm for further inspection. Compositional data such as chemical concentrations can be visualised by conventional ‘normal’ statistics after logratio transformation. The Euclidean distance in logratio space is called the Aitchison distance in compositional data space. Classical MDS using this distance is equivalent to Principal Component Analysis. Finally, point-counting data combine aspects of compositional data analysis with multinomial sampling statistics. The Chi-square distance is the natural way to quantify the dissimilarity between multiple point-counting samples. MDS analysis using the Chi-square distance is equivalent to Correspondence Analysis, which is akin to PCA for categorical data.

However, there are some provenance proxies that do not easily fit into these three categories. Varietal studies using the chemical composition of single grains of heavy minerals combine aspects of compositional and distributional data [3,60]. Similarly, paired U-Pb ages and Hf-isotope compositions in zircon [1] do not easily fit inside the distributional data class described above. Using the tools provided by the provenance package, such data can be processed by procustes analysis or 3-way MDS (Section 11). Thus, U-Pb and $ϵ$(Hf)-distributions, say, could be entered into the indscal function as separate entities. However, by doing so, the single-grain link between the two datasets would be lost. Alternative approaches may be pursued to address this issue, and new dissimilarity measures could be developed for this hybrid data type. Novel approaches to matrix decomposition may be a way forward in this direction [8,67,68].