What Useful Information Can the Principal Analysis Provide?

1. Introduction

Large datasets are increasingly widespread in many disciplines. In club to translate such datasets, methods are required to drastically reduce their dimensionality in an interpretable mode, such that most of the information in the data is preserved. Many techniques have been developed for this purpose, but principal component analysis (PCA) is one of the oldest and most widely used. Its idea is simple—reduce the dimensionality of a dataset, while preserving every bit much 'variability' (i.east. statistical information) as possible.

Although it is used, and has sometimes been reinvented, in many different disciplines it is, at heart, a statistical technique and hence much of its development has been by statisticians.

This means that 'preserving as much variability as possible' translates into finding new variables that are linear functions of those in the original dataset, that successively maximize variance and that are uncorrelated with each other. Finding such new variables, the primary components (PCs), reduces to solving an eigenvalue/eigenvector problem. The earliest literature on PCA dates from Pearson [1] and Hotelling [2], merely information technology was not until electronic computers became widely bachelor decades after that information technology was computationally feasible to use it on datasets that were non trivially small. Since then its utilise has burgeoned and a large number of variants have been adult in many unlike disciplines. Substantial books have been written on the subject area [3,4] and there are even whole books on variants of PCA for special types of information [5,6]. In §2, the formal definition of PCA will be given, in a standard context, together with a derivation showing that it tin exist obtained every bit the solution to an eigenproblem or, alternatively, from the singular value decomposition (SVD) of the (centred) information matrix. PCA tin be based on either the covariance matrix or the correlation matrix. The choice between these analyses will be discussed. In either case, the new variables (the PCs) depend on the dataset, rather than being pre-divers basis functions, and so are adaptive in the broad sense. The main uses of PCA are descriptive, rather than inferential; an case will illustrate this.

Although for inferential purposes a multivariate normal (Gaussian) distribution of the dataset is usually assumed, PCA as a descriptive tool needs no distributional assumptions and, as such, is very much an adaptive exploratory method which can exist used on numerical information of diverse types. Indeed, many adaptations of the basic methodology for dissimilar data types and structures have been developed, two of which will exist described in §3a,d. Some techniques give simplified versions of PCs, in order to aid estimation. Two of these are briefly described in §3b, which also includes an case of PCA, together with a simplified version, in atmospheric science, illustrating the adaptive potential of PCA in a specific context. Section 3c discusses one of the extensions of PCA that has been most agile in contempo years, namely robust PCA (RPCA). The explosion in very large datasets in areas such as image assay or the analysis of Web data has brought virtually of import methodological advances in data analysis which often find their roots in PCA. Each of §3a–d gives references to recent work. Some concluding remarks, emphasizing the breadth of application of PCA and its numerous adaptations, are made in §4.

ii. The bones method

(a) Principal component assay as an exploratory tool for data analysis

The standard context for PCA equally an exploratory information analysis tool involves a dataset with observations on p numerical variables, for each of n entities or individuals. These data values define p n-dimensional vectors x ₁,…,10 _p or, equivalently, an north×p information matrix 10, whose jth column is the vector 10 _j of observations on the jth variable. We seek a linear combination of the columns of matrix X with maximum variance. Such linear combinations are given by , where a is a vector of constants a _ane,a ₂,…,a _p . The variance of whatsoever such linear combination is given past var(X a)=a′S a, where S is the sample covariance matrix associated with the dataset and ′ denotes transpose. Hence, identifying the linear combination with maximum variance is equivalent to obtaining a p-dimensional vector a which maximizes the quadratic form a′Southward a. For this trouble to take a well-defined solution, an additional restriction must be imposed and the most common restriction involves working with unit of measurement-norm vectors, i.due east. requiring a′a=1. The problem is equivalent to maximizing a′Southward a−λ(a′a−i), where λ is a Lagrange multiplier. Differentiating with respect to the vector a, and equating to the goose egg vector, produces the equation

2.1

Thus, a must be a (unit of measurement-norm) eigenvector, and λ the respective eigenvalue, of the covariance matrix South. In particular, we are interested in the largest eigenvalue, λ _i (and corresponding eigenvector a ₁), since the eigenvalues are the variances of the linear combinations defined by the respective eigenvector a: var(10 a)=a′S a=λ a′a=λ. Equation (2.1) remains valid if the eigenvectors are multiplied by −one, and and then the signs of all loadings (and scores) are arbitrary and just their relative magnitudes and sign patterns are meaningful.

Any p×p real symmetric matrix, such as a covariance matrix S, has exactly p real eigenvalues, λ _thousand (k=one,…,p), and their respective eigenvectors tin can exist divers to grade an orthonormal ready of vectors, i.due east. a′ _m a _k′=1 if k=k′ and cipher otherwise. A Lagrange multipliers arroyo, with the added restrictions of orthogonality of different coefficient vectors, can as well be used to testify that the full set of eigenvectors of S are the solutions to the problem of obtaining upwards to p new linear combinations , which successively maximize variance, subject to uncorrelatedness with previous linear combinations [4]. Uncorrelatedness results from the fact that the covariance between ii such linear combinations, X a _k and Ten a _m′, is given by a′ _thousand′ S a _k =λ _k a′ _k′ a _k =0 if k′≠k.

It is these linear combinations X a _k that are called the primary components of the dataset, although some authors confusingly likewise use the term 'main components' when referring to the eigenvectors a _{one thousand} . In standard PCA terminology, the elements of the eigenvectors a _m are commonly called the PC loadings, whereas the elements of the linear combinations X a _k are chosen the PC scores, as they are the values that each individual would score on a given PC.

It is common, in the standard approach, to define PCs as the linear combinations of the centred variables x* _j , with generic element , where denotes the mean value of the observations on variable j. This convention does not change the solution (other than centring), since the covariance matrix of a prepare of centred or uncentred variables is the same, but it has the advantage of providing a straight connectedness to an culling, more geometric arroyo to PCA.

Cogent past X* the n×p matrix whose columns are the centred variables 10* _j , nosotros have

2.ii

Equation (2.ii) links upwards the eigendecomposition of the covariance matrix S with the singular value decomposition of the column-centred data matrix X*. Any arbitrary matrix Y of dimension northward×p and rank r (necessarily, ) can be written (due east.k. [4]) as

2.3

where U,A are north×r and p×r matrices with orthonormal columns (U′U=I _r =A′A, with I _r the r×r identity matrix) and Fifty is an r×r diagonal matrix. The columns of A are called the right singular vectors of Y and are the eigenvectors of the p×p matrix Y′Y associated with its non-nothing eigenvalues. The columns of U are called the left singular vectors of Y and are the eigenvectors of the north×due north matrix Y Y′ that represent to its non-zilch eigenvalues. The diagonal elements of matrix L are called the singular values of Y and are the non-negative square roots of the (common) not-nothing eigenvalues of both matrix Y′Y and matrix Y Y′. We assume that the diagonal elements of L are in decreasing order, and this uniquely defines the lodge of the columns of U and A (except for the case of equal atypical values [4]). Hence, taking Y=X*, the correct singular vectors of the cavalcade-centred data matrix X* are the vectors a _yard of PC loadings. Due to the orthogonality of the columns of A, the columns of the matrix product X*A=ULA′A=UL are the PCs of X*. The variances of these PCs are given by the squares of the atypical values of Ten*, divided by n−one. Equivalently, and given (ii.2) and the above properties,

2.four

where 50 ² is the diagonal matrix with the squared singular values (i.due east. the eigenvalues of (n−1)S). Equation (2.4) gives the spectral decomposition, or eigendecomposition, of matrix (n−1)S. Hence, PCA is equivalent to an SVD of the column-centred information matrix X*.

The properties of an SVD imply interesting geometric interpretations of a PCA. Given whatsoever rank r matrix Y of size n×p, the matrix Y _q of the same size, but of rank q<r, whose elements minimize the sum of squared differences with corresponding elements of Y is given [vii] by

2.5

where L _q is the q×q diagonal matrix with the first (largest) q diagonal elements of L and U _q , A _q are the due north×q and p×q matrices obtained by retaining the q corresponding columns in U and A.

In our context, the north rows of a rank r column-centred information matrix 10* define a scatterplot of n points in an r-dimensional subspace of , with the origin as the centre of gravity of the scatterplot. The higher up issue implies that the 'best' n-point approximation to this scatterplot, in a q-dimensional subspace, is given by the rows of X* _q , divers as in equation (two.v), where 'best' means that the sum of squared distances between corresponding points in each scatterplot is minimized, every bit in the original arroyo by Pearson [1]. The system of q axes in this representation is given past the first q PCs and defines a primary subspace. Hence, PCA is at eye a dimensionality-reduction method, whereby a gear up of p original variables can be replaced by an optimal set of q derived variables, the PCs. When q=two or q=3, a graphical approximation of the due north-betoken scatterplot is possible and is oft used for an initial visual representation of the full dataset. It is important to note that this result is incremental (hence adaptive) in its dimensions, in the sense that the best subspace of dimension q+1 is obtained past adding a farther cavalcade of coordinates to those that defined the best q-dimensional solution.

The quality of any q-dimensional approximation can be measured past the variability associated with the set up of retained PCs. In fact, the sum of variances of the p original variables is the trace (sum of diagonal elements) of the covariance matrix S. Using simple matrix theory results it is straightforward to bear witness that this value is also the sum of the variances of all p PCs. Hence, the standard mensurate of quality of a given PC is the proportion of total variance that information technology accounts for,

2.six

where tr(S) denotes the trace of S. The incremental nature of PCs besides means that we can speak of a proportion of total variance explained by a set of PCs (usually, but not necessarily, the first q PCs), which is often expressed as a percent of full variance deemed for: .

It is common do to use some predefined percentage of total variance explained to determine how many PCs should be retained (70% of total variability is a mutual, if subjective, cut-off point), although the requirements of graphical representation frequently lead to the use of just the offset two or three PCs. Even in such situations, the percentage of total variance accounted for is a fundamental tool to assess the quality of these low-dimensional graphical representations of the dataset. The emphasis in PCA is almost ever on the first few PCs, just there are circumstances in which the concluding few may be of interest, such as in outlier detection [4] or some applications of image analysis (see §3c).

PCs can also be introduced as the optimal solutions to numerous other problems. Optimality criteria for PCA are discussed in detail in numerous sources (see [4,eight,ix], amidst others). McCabe [10] uses some of these criteria to select optimal subsets of the original variables, which he calls master variables. This is a dissimilar, computationally more complex, problem [11].

(b) Example: fossil teeth data

PCA has been applied and found useful in very many disciplines. The ii examples explored hither and in §3b are very different in nature. The offset examines a dataset consisting of ix measurements on 88 fossil teeth from the early on mammalian insectivore Kuehneotherium, while the 2nd, in §3b, is from atmospheric science.

Kuehneotherium is 1 of the primeval mammals and remains have been plant during quarrying of limestone in S Wales, UK [12]. The bones and teeth were done into fissures in the rock, nigh 200 one thousand thousand years ago, and all the lower tooth teeth used in this analysis are from a unmarried fissure. Still, it looked possible that there were teeth from more than than 1 species of Kuehneotherium in the sample.

Of the nine variables, iii mensurate aspects of the length of a tooth, while the other vi are measurements related to height and width. A PCA was performed using the prcomp command of the R statistical software [13]. The first two PCs business relationship for 78.viii% and 16.seven%, respectively, of the total variation in the dataset, so the two-dimensional scatter-plot of the 88 teeth given past figure 1 is a very good approximation to the original scatter-plot in nine-dimensional space. It is, past definition, the best variance-preserving 2-dimensional plot of the data, representing over 95% of total variation. All of the loadings in the first PC have the same sign, and so it is a weighted average of all variables, representing 'overall size'. In figure 1, big teeth are on the left and small teeth on the right. The 2d PC has negative loadings for the three length variables and positive loadings for the other six variables, representing an attribute of the 'shape' of teeth. Fossils near the summit of figure ane take smaller lengths, relative to their heights and widths, than those towards the lesser. The relatively meaty cluster of points in the bottom half of figure 1 is idea to correspond to a species of Kuehneotherium, while the broader group at the top cannot be assigned to Kuehneotherium, but to some related, only as however unidentified, creature.

Figure i. The two-dimensional master subspace for the fossil teeth data. The coordinates in either or both PCs may switch signs when dissimilar software is used.

• Download figure
• Open up in new tab
• Download PowerPoint

(c) Some central bug

(i) Covariance and correlation matrix main component analysis

So far, PCs have been presented as linear combinations of the (centred) original variables. However, the properties of PCA have some undesirable features when these variables have unlike units of measurement. While in that location is nothing inherently incorrect, from a strictly mathematical point of view, with linear combinations of variables with dissimilar units of measurement (their apply is widespread in, for example, linear regression), the fact that PCA is defined past a criterion (variance) that depends on units of measurement implies that PCs based on the covariance matrix S will change if the units of measurement on one or more of the variables change (unless all p variables undergo a common change of scale, in which case the new covariance matrix is merely a scalar multiple of the old one, hence with the aforementioned eigenvectors and the same proportion of total variance explained past each PC). To overcome this undesirable feature, it is common do to begin by standardizing the variables. Each data value 10 _ij is both centred and divided by the standard divergence southward _j of the n observations of variable j,

2.vii

Thus, the initial data matrix X is replaced with the standardized data matrix Z, whose jth column is vector z _j with the n standardized observations of variable j (2.7). Standardization is useful because most changes of scale are linear transformations of the data, which share the same ready of standardized data values.

Since the covariance matrix of a standardized dataset is merely the correlation matrix R of the original dataset, a PCA on the standardized data is too known as a correlation matrix PCA. The eigenvectors a _k of the correlation matrix R define the uncorrelated maximum-variance linear combinations of the standardized variables z ₁,…,z _p . Such correlation matrix PCs are not the aforementioned every bit, nor are they direct related to, the covariance matrix PCs divers previously. Besides, the percentage variance deemed for by each PC will differ and, quite oftentimes, more correlation matrix PCs than covariance matrix PCs are needed to account for the aforementioned percentage of total variance. The trace of a correlation matrix R is merely the number p of variables used in the analysis, hence the proportion of total variance deemed for past whatsoever correlation matrix PC is just the variance of that PC divided past p. The SVD approach is as well valid in this context. Since (due north−1)R=Z′Z, an SVD of the standardized data matrix Z amounts to a correlation matrix PCA of the dataset, forth the lines described after equation (ii.2).

Correlation matrix PCs are invariant to linear changes in units of measurement and are therefore the appropriate option for datasets where different changes of scale are believable for each variable. Some statistical software assumes by default that a PCA means a correlation matrix PCA and, in some cases, the normalization used for the vectors of loadings a ₁₀₀₀ of correlation matrix PCs is not the standard a′ _k a _chiliad =1. In a correlation matrix PCA, the coefficient of correlation between the jth variable and the kth PC is given past (run across [four])

ii.8

Thus, if the normalization is used instead of a′ _thousand a=ane, the coefficients of the new loading vectors are the correlations between each original variable and the kth PC.

In the fossil teeth information of §2b, all nine measurements are in the same units, so a covariance matrix PCA makes sense. A correlation matrix PCA produces similar results, since the variances of the original variable do not differ very much. The get-go two correlation matrix PCs account for 93.7% of total variance. For other datasets, differences tin be more substantial.

(2) Biplots

One of the most informative graphical representations of a multivariate dataset is via a biplot [xiv], which is fundamentally connected to the SVD of a relevant information matrix, and therefore to PCA. A rank q approximation Ten* _q of the full column-centred data matrix X*, defined past (two.5), is written as X* _q =GH′, where G=U _q and H=A _q L _q (although other options are possible, meet [iv]). The n rows m _i of matrix G define graphical markers for each private, which are usually represented by points. The p rows h _j of matrix H define markers for each variable and are commonly represented by vectors. The properties of the biplot are best discussed bold that q=p, although the biplot is defined on a low-rank approximation (usually q=2), enabling a graphical representation of the markers. When q=p the biplot has the following properties:

• — The cosine of the angle betwixt whatever two vectors representing variables is the coefficient of correlation between those variables; this is a straight result of the fact that the matrix of inner products between those markers is HH′=AL ² A′=(north−1)S (2.four), so that inner products betwixt vectors are proportional to covariances (variances for a common vector).

• — Similarly, the cosine of the bending between any vector representing a variable and the centrality representing a given PC is the coefficient of correlation between those 2 variables.

• — The inner production betwixt the markers for individual i and variable j gives the (centred) value of individual i on variable j. This is a direct issue of the fact that GH′=X*. The applied implication of this result is that orthogonally projecting the signal representing individual i onto the vector representing variable j recovers the (centred) value .

• — The Euclidean altitude betwixt the markers for individuals i and i′ is proportional to the Mahalanobis distance between them (see [iv] for more details).

Every bit stated above, these results are only exact if all q=p dimensions are used. For q<p, the results are merely approximate and the overall quality of such approximations tin be measured past the percentage of variance explained past the q largest variance PCs, which were used to build the marker matrices Thou and H.

Effigy ii gives the biplot for the correlation matrix PCA of the fossil teeth data of §2b. The variable markers are displayed as arrows and the molar markers as numbers. The grouping of three nigh horizontal and very tightly knit variable markers for two width variables and one height variable, WIDTH, HTMDT and TRIWIDTH, suggests a grouping of highly correlated variables, which are likewise strongly correlated with the start PC (represented past the horizontal axis). The very loftier proportion of variability explained by the ii-dimensional principal subspace provides solid grounds for these conclusions. In fact, the smallest of the 3 truthful coefficients of correlation between these three variables is 0.944 (HTMDT and TRIWIDTH), and the smallest magnitude correlation between PC1 and any of these variables is 0.960 (TRIWIDTH). The sign departure in PC2 loadings between the iii length variables (towards the lesser left of the plot) and the other variables is clearly visible. Projecting the marking for private 58 onto the positive directions of all variable markers suggests that fossil tooth 58 (on the left of the biplot) is a big molar. Inspection of the data matrix confirms that it is the largest private on half-dozen of the nine variables, and close to largest on the remaining three. Likewise, individuals 85–88 (on the right) are minor-sized teeth. Individuals whose markers are close to the origin have values close to the mean for all variables.

Figure 2. Biplot for the fossil teeth data (correlation matrix PCA), obtained using R'due south biplot command. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

(3) Centrings

As was seen in §2, PCA amounts to an SVD of a column-centred data matrix. In some applications [xv], centring the columns of the data matrix may be considered inappropriate. In such situations, it may be preferred to avert whatever pre-processing of the information and to field of study the uncentred data matrix to an SVD or, equivalently, to carry out the eigendecomposition of the matrix of non-centred second moments, T, whose eigenvectors ascertain linear combinations of the uncentred variables. This is oftentimes referred to as an uncentred PCA and there has been an unfortunate tendency in some fields to equate the name SVD merely with this uncentred version of PCA.

Uncentred PCs are linear combinations of the uncentred variables which successively maximize not-central second moments, subject to having their crossed non-central second moments equal to zero. Except when the vector of cavalcade means (i.e. the center of gravity of the original n-point scatterplot in p-dimensional space) is nigh zero (in which example centred and uncentred moments are like), information technology is not immediately intuitive that in that location should exist similarities between both variants of PCA. Cadima & Jolliffe [15] have explored the relations between the standard (column-centred) PCA and uncentred PCA and establish them to be closer than might be expected, in detail when the size of vector is large. Information technology is oft the instance that there are great similarities between many eigenvectors and (absolute) eigenvalues of the covariance matrix Due south and the corresponding matrix of non-centred second moments, T.

In some applications, row centrings, or both row- and column-centring (known every bit double-centring) of the data matrix, have been considered advisable. The SVDs of such matrices give rise to row-centred and doubly centred PCA, respectively.

(iv) When north<p

Datasets where there are fewer observed entities than variables (n<p) are becoming increasingly frequent, cheers to the growing ease of observing variables, together with the loftier costs of repeating observations in some contexts (such as microarrays [16]). For example, [17] has an example in genomics in which n=59 and p=21 225.

In general, the rank of an north×p data matrix is . If the data matrix has been column-centred, it is . When due north<p, information technology is the number of observed individuals, rather than the number of variables, that usually determines the matrix rank. The rank of the column-centred data matrix X* (or its standardized counterpart Z) must equal the rank of the covariance (or correlation) matrix. The practical implication of this is that in that location are only r not-zero eigenvalues; hence r PCs explicate all the variability in the dataset. Cypher prevents the employ of PCA in such contexts, although some software, equally is the case with R'southward princomp (but not the prcomp) control, may balk at such datasets. PCs can exist adamant as usual, past either an SVD of the (centred) information matrix or the eigenvectors/values of the covariance (or correlation) matrix.

Contempo enquiry (e.yard. [18,19]) has examined how well underlying 'population' PCs are estimated by the sample PCs in the example where n≪p, and it is shown that in some circumstances there is little resemblance between sample and population PCs. However, the results are typically based on a model for the information which has a very small number of structured PCs, and very many noise dimensions, and which has some links with recent work in RPCA (meet §3c).

3. Adaptations of chief component analysis

The bones idea of PCA, leading to depression-dimensional representations of large datasets in an adaptive and insightful way, is simple. However, the subsections in §2 have shown a number of subtleties that add together some complexity. Going further, there are many ways to adjust PCA to reach modified goals or to analyse data of different types. Because PCA is used in a big number of areas, research into modifications and adaptations is spread over literatures from many disciplines. Four such adaptations, chosen fairly arbitrarily from the many that be, namely functional PCA, modifications of PCA to simplify interpretations, RPCA and symbolic data PCA, are described in the following subsections. Other adaptations are briefly mentioned in §4.

(a) Functional principal component analysis

In some applications, such as chemic spectroscopy, observations are functional in nature, irresolute with some continuous variable which, for simplicity, we assume is time. The dataset is then a collection of n functions x _i (t).

How to incorporate such functional features in the analysis is the goal of functional information analysis [twenty]. Early on work on functional PCA (e.g. [21]) performed a standard PCA on an n×p data matrix obtained past sampling n curves x _i (t) at each of p points in time (t _j , with j=1,…,p), then that the element in row i, column j, of the data matrix is x _i (t _j ). The resulting p-dimensional vectors of loadings from a PCA of this data matrix are then viewed as sampled primary functions, which tin can exist smoothed to recover functional form and can exist interpeted as principal sources of variability in the observed curves [twenty]. The above approach does not make explicit use of the functional nature of the northward observations x _i (t). To do so requires adapting concepts. In the standard setting, we consider linear combinations of p vectors, which produce new vectors. Each element of the new vectors is the issue of an inner product of row i of the data matrix, (10 _i1,x _i2,…,ten _ip ), with a p-dimensional vector of weights, a=(a ₁,…,a _p ): . If rows of the information matrix go functions, a functional inner product must be used instead, between a 'loadings function', a(t), and the ith functional observation, ten _i (t). The standard functional inner product is an integral of the form , on some appropriate meaty interval. Likewise, the analogue of the p×p covariance matrix S is a bivariate function S(south,t) which, for any 2 given time instants s and t, returns the respective covariance, defined as

3.one

where is the mean function and is the ith centred role.

The analogue of the eigen-equation (2.ane) involves an integral transform, which reflects the functional nature of S(s,t) and of inner products

3.2

The eigenfunctions a(t) which are the analytic solutions of equation (3.2) cannot, in general, be determined. Ramsay & Silverman [20] hash out estimate solutions based on numerical integration. An alternative approach, which they explore in greater particular, involves the supposition that the curves ten _i (t) can exist written as linear combinations of a set of G ground functions ϕ ₁(t),…,ϕ _G (t), then that, for any information part i,

3.iii

These basis functions can be called to reflect characteristics that are considered relevant in describing the observed functions. Thus, Fourier series functions may exist called to describe periodic traits and splines for more than general trends (B-splines are recommended). Other footing functions that take been used and can be considered are wavelets, exponential, power or polynomial bases. In theory, other bases, adapted to specific properties of a given set of observed functions, may exist considered, although the computational problems that arise from any such choice must be kept in heed. The advantage of the basis function approach lies in the simplification of the expressions given previously. Cogent the n-dimensional vector of functions x _i (t) as x(t), the G-dimensional vector of footing functions as ϕ(t) and the northward×G matrix of coefficients c _ij as C, the n information functions in equation (3.3) tin can be written as a single equation ten(t)=C ϕ(t). The eigenfunction a(t) can also be written in terms of the basis functions, with a(t)=ϕ(t)′b for some Thou-dimensional vector of coefficients b=(b ₁,…,b _G ). Assuming furthermore that ten(t) and ϕ(t) are centred, the covariance function at time (s,t) becomes

and eigen-equation (3.ii) becomes, after some algebraic manipulation (run across [4,xx] for details),

where Westward is the G×Thou matrix of inner products between the basis functions. Since this equation must agree for all values of south, it reduces to

three.4

If the basis functions are orthonormal, W is the Yard×G identity matrix and we stop up with a standard eigenvalue trouble which provides the solutions a(t)=ϕ(t)′b to equation (three.two).

Ramsay & Silverman [xx] further explore methods in which data functions x _i (t) are viewed every bit solutions to differential equations, an approach which they call main differential analysis, in order to highlight its close connections with PCA.

Research on functional PCA has continued apace since the publication of Ramsay and Silverman's comprehensive text. Often this research is parallel to, or extends, similar ideas for data of non-functional class. For instance, deciding how many PCs to retain is an important topic. A large number of suggestions have been made for doing so [4] and many pick criteria are based on intuitive or descriptive ideas, such equally the obvious 'proportion of total variance'. Other approaches are based on models for PCs. The problem of 'how many functional PCs?' is addressed in [22] using a model-based approach and criteria based on data theory.

As with other statistical techniques, it is possible that a few outlying observations may have a disproportionate event on the results of a PCA. Numerous suggestions have been fabricated for making PCA more than robust to the presence of outliers for the usual information structure (see [4] and also §3c). One suggestion, using and then-chosen South-estimators, is extended to functional PCA in [23].

Sometimes, as well as correlations between the p variables, there is a dependence structure between the n observations. A 'dynamic' version of functional PCA is proposed in [24], which is relevant when there are correlations between the observed curves, also every bit the obvious correlation within the curves. Information technology is based on an idea commencement suggested by Brillinger [25] for vector time serial and uses frequency domain assay.

(b) Simplified principal components

PCA gives the best possible representation of a p-dimensional dataset in q dimensions (q<p) in the sense of maximizing variance in q dimensions. A disadvantage is, yet, that the new variables that it defines are usually linear functions of all p original variables. Although information technology was possible to interpret the first two PCs in the fossil teeth instance, information technology is frequently the case for larger p that many variables take non-trivial coefficients in the first few components, making the components difficult to interpret. A number of adaptations of PCA have been suggested that endeavor to brand interpretation of the q dimensions simpler, while minimizing the loss of variance due to not using the PCs themselves. There is a trade-off between interpretability and variance. Two such classes of adaptations are briefly described here.

Rotation. The thought of rotating PCs is borrowed from factor assay [26] (a different method, which is sometimes dislocated with PCA—run across [4] for a fuller discussion). Suppose, as before, that A _q is the p×q matrix, whose columns are the loadings of the kickoff q PCs. Then XA _q is the due north×q matrix whose columns are the scores on the first q PCs for the n observations. Now let T be an orthogonal (q×q) matrix. Multiplication of A _q past T performs an orthogonal rotation of the axes within the space spanned past the commencement q PCs, so that B _q =A _q T is a p×q matrix whose columns are loadings of q rotated PCs. The matrix XB _q is an n×q matrix containing the corresponding rotated PC scores. Whatsoever orthogonal matrix T could exist used to rotate the components, but if it is desirable to brand the rotated components easy to interpret, then T is chosen to optimize some simplicity criterion. A number of such criteria have been suggested, including some that include non-orthogonal (oblique) rotation [26]. The most popular is perhaps the varimax criterion in which an orthogonal matrix T is chosen to maximize , where b _jk is the (j,grand)thursday element of B _q .

Rotation can considerably simplify interpretation and, when viewed with respect to the q-dimensional space that is rotated, no variance is lost, equally the sum of variances of the q rotated components is the aforementioned as for the unrotated components. What is lost is the successive maximization of the unrotated PCs, and then that the total variance of the q components is more than evenly distributed between components after rotation.

Drawbacks of rotation include the need to choose from the plethora of possible rotation criteria, though this choice ofttimes makes less difference than the selection of how many components to rotate. The rotated components can await quite different if q is increased by i, whereas the successively divers nature of unrotated PCs means that this does not happen.

Adding a constraint. Another arroyo to simplification of PCs is to impose a constraint on the loadings of the new variables. Again, there are a number of variants of this arroyo, one of which adapts the LASSO (to the lowest degree absolute shrinkage and selection operator) arroyo from linear regression [27]. In this approach, called SCoTLASS (simplified component technique–LASSO), components are found which successively solve the same optimization trouble every bit PCA, but with the additional constraint , where τ is a tuning parameter. For , the constraint has no effect and PCs are obtained, but every bit τ decreases more and more loadings are driven to zero, thus simplifying interpretion. These simplified components necessarily account for less variance than the respective number of PCs, and ordinarily several values of τ are tried to determine a skilful trade-off between added simplicity and loss of variance.

A departure betwixt the rotation and constraint approaches is that the latter has the reward for interpretation of driving some loadings in the linear functions exactly to zero, whereas rotation normally does not. Adaptations of PCA in which many coefficients are exactly zero are generally known as thin versions of PCA, and there has been a substantial corporeality of research on such PCs in recent years. A skillful review of such work tin be found in Hastie et al. [28] (see as well §3c).

A related technique to SCoTLASS adds a penalisation function to the variance benchmark maximized, so that the optimization problem becomes to successively find a _{one thousand} , k=1,ii,…,p, that maximize , bailiwick to a′ _k a _k =one, where ψ is a tuning parameter [29]. One of the nowadays authors has recently reviewed a paper in which it is demonstrated that these apparently equivalent constraint and penalization approaches actually have quite distinct properties.

The original SCoTLASS optimization trouble is not-convex and is likewise not solvable by simple iterative algorithms, although it is possible to re-express SCoTLASS as an equivalent, though still non-convex, optimization problem for which simple algorithms tin can be used [30]. Another arroyo, due to d'Aspremont et al. [31], reformulates SCoTLASS in a more circuitous manner, but then drops one of the constraints in this new formulation in order to make the trouble convex.

Achieving sparsity is important for large p and specially when n≪p. A number of authors take investigated versions of sparse PCA for this situation using models for the data in which the vast majority of the variables are completely unstructured racket [18,nineteen,32]. These papers and others suggest and investigate the backdrop of algorithms for estimating sparse PCs when information are generated from their models. Lee et al. [17] utilize a different blazon of model, this fourth dimension a random effects model for PC loadings, to derive an alternative penalty part to that used past SCoTLASS, giving some other sparse PCA method. Additionally incorporating shrinkage of eigenvalues leads to yet another method, deemed super-sparse PCA in [17]. Comparisons are given between their methods, SCoTLASS and the elastic net [28] for imitation data and a large genomic example.

(i) Example: sea-level pressure data

I subject in which PCA has been widely used is atmospheric science. It was first suggested in that field past Obukhov [33] and Lorenz [34] and, uniquely to that discipline, it is usually known as empirical orthogonal function (EOF) assay. The book by Preisendorfer & Mobley [35] discusses many aspects of PCA in the context of meteorology and oceanography.

The format of the information in atmospheric science is unlike from that of most other disciplines. This example is taken from [36]. The data consist of measurements of winter (December, January and February) monthly mean sea-level pressure level (SLP) over the Northern Hemisphere north of 20° North. The dataset is available on a 2.five°×two.5° regular grid and spans the catamenia from January 1948 to December 2000. Some preprocessing is done to adjust for the annual bike and the different areas covered by grid squares at different latitudes. In many atmospheric science examples, the variables are measurements at grid points, and the loadings, known as EOFs, are displayed as smooth spatial patterns, as in figure three for the commencement two correlation-based EOFs for the SLP data [36]. There are 1008 variables (grid-points) in this dataset, and the get-go two PCs account for 21% and thirteen% of the variation in these 1008 variables. Figure iii gives a pattern which is commonly known every bit the Arctic Oscillation (AO). Information technology is a measure of the n–s pressure gradient in the Atlantic Ocean and, to a bottom extent, in the Pacific Bounding main and is a major source of variation in atmospheric condition patterns. The second EOF is dominated by variation in the Pacific Body of water. The PCs for examples of this blazon are time serial and then the first PC, for instance, will display which years have high values of the AO and which have depression values.

Figure 3. (a,b) The first two correlation-based EOFs for the SLP information business relationship for 21% and 13% of total variation. (Adapted from [36].)

• Download figure
• Open in new tab
• Download PowerPoint

Figure 4 shows simplified EOFs based on SCoTLASS [36]. The master deviation from the EOFs in figure 3 is for the kickoff EOF, which is now completely dominated by the north–south pressure level gradient in the Atlantic (the North Atlantic Oscillation) with exactly nada loadings for many grid-points. The simplification is paid for past a reduction in percentage of variation explained for the corresponding simplified PC (17% compared with 21%). The second simplified PC is very similar to the original second EOF, also explaining 13% of variation.

Figure 4. (a,b) LASSO-based simplified EOFs for the SLP data. Grey areas are filigree-points with exactly aught loadings. (Adapted from [36].)

• Download figure
• Open in new tab
• Download PowerPoint

(c) Robust principal component assay

By its very nature, PCA is sensitive to the presence of outliers and therefore likewise to the presence of gross errors in the datasets. This has led to attempts to define robust variants of PCA and the expression RPCA has been used for dissimilar approaches to this trouble. Early work past Huber [37,38] discussed robust alternatives to covariance or correlation matrices and means in which they tin exist used to define robust PCs. This work was extended in [39,xl]; see likewise [41].

The demand for methods to deal with very large datasets in areas such as epitome processing, machine learning, bioinformatics or Web data analysis has generated a recent renewed involvement in robust variants of PCA and has led to one of the most vigorous lines of inquiry in PCA-related methods. A word of this event can be plant in [42]. Wright et al. [43] defined RPCA equally a decomposition of an n×p data matrix X into a sum of two n×p components: a low-rank component L and a sparse component S. More precisely, a convex optimization problem was defined equally identifying the matrix components of X=L+South that minimize a linear combination of two dissimilar norms of the components:

3.5

where , the sum of the singular values of 50, is the nuclear norm of Fifty, and is the ℓ₁-norm of matrix S. The motivation for such a decomposition is the fact that, in many applications, low-rank matrices are associated with a full general pattern (east.thousand. the 'correct' data in a corrupted dataset, a face in facial recognition, or a groundwork image in video surveillance information), whereas a thin matrix is associated with disturbances (eastward.g. corrupted information values, furnishings of light or shading in facial recognition, a moving object or person in the foreground of information surveillance images). Sparse components are besides called 'noise', in what tin be confusing terminology since in some applications information technology is precisely the 'racket' component that is of involvement. Trouble 3.5 has obvious points of contact with some of the give-and-take in §3b. Candès et al. [44] return to this problem, as well chosen master component pursuit, and requite theoretical results proving that, under not very stringent conditions, it is possible to exactly recover the depression-rank and sparse components with high probability and that the choice of works well in a general setting, avoiding the demand to cull a tuning parameter. Results are extended to the case of information matrices with missing values. Algorithms for the identification of the components are also discussed in [44], an important consequence given the computational complexity involved. Further variations consider more complex structures for the 'noise' component. Some such proposals are reviewed in [45], where the results of alternative algorithms in the presence of different types of 'dissonance' are compared in the context of epitome-processing and facial recognition problems. Their results show that classical PCA performs adequately well, when compared with these new methods, in terms of both time and the quality of low-rank solutions that are produced. A fairly recent review of work in this expanse can be found in [46].

(d) Symbolic information master component analysis

There is a recent body of work with so-called symbolic data, which is a general designation for more complex data structures, such as intervals or histograms [47,48].

Interval data arise when one wishes to retain a measure of underlying variability in the observations. This may occur if nosotros wish to reflect the lack of precision of a measuring instrument or, more fundamentally, because the data are summary observations for which associated variability is considered inherent to the measurement. This is ofttimes the case when each observation corresponds to a group, rather than an private, equally would be the case with measurements on species, for which a range of values is considered part of the group value. If all p observed variables are of this blazon, each observation is represented past a hyper-rectangle, rather than a bespeak, in p-dimensional space. Extensions of PCA for such data [47,49] seek PCs that are also of interval type, and which therefore also reflect ranges of values.

Another mutual type of symbolic data is given by histograms, which tin be considered a generalization of interval-valued data where for each observation there are several intervals (the histogram bins) and associated frequencies. A contempo review [50] covers several proposed definitions of PCA-type analyses for histogram data. Virtually of them require the definition of concepts such equally distances between histograms (the Wasserstein distance being a common choice) or the sum and hateful of histograms.

4. Conclusion

Although PCA in its standard form is a widely used and adaptive descriptive data assay tool, it as well has many adaptations of its ain that make information technology useful to a wide variety of situations and data types in numerous disciplines. Adaptations of PCA accept been proposed, among others, for binary data, ordinal data, compositional data, discrete data, symbolic information or data with special structure, such every bit time series [4] or datasets with common covariance matrices [6,40]. PCA or PCA-related approaches have also played an important direct role in other statistical methods, such every bit linear regression (with primary component regression [iv]) and even simultaneous clustering of both individuals and variables [51]. Methods such as correspondance assay, approved correlation analysis or linear discriminant assay may be only loosely connected to PCA, but, insofar as they are based on factorial decompositions of certain matrices, they share a mutual approach with PCA. The literature on PCA is vast and spans many disciplines. Space constraints mean that it has been explored very superficially here. New adaptations and methodological results, as well as applications, are still appearing.

Data accessibility

The fossil teeth data are bachelor from I.T.J. The atmospheric scientific discipline data were taken from the publicly accessible NCEP/NCAR reanalysis database (encounter [36] for details).

Authors' contributions

Both authors were every bit involved in drafting the manuscript.

Competing interests

We have no competing interests.

Funding

Research past J.C. is partially supported by the Portuguese Science Foundation FCT - PEst-OE/MAT/UI0006/2014.

Acknowledgements

We thank Pamela Gill and Abdel Hannachi for helpful discussions regarding their data and results.

Footnotes

One contribution of 13 to a theme consequence 'Adaptive data assay: theory and applications'.

© 2022 The Author(south)

Published by the Regal Society. All rights reserved.

References

• 1
  Pearson Yard

  . 1901 On lines and planes of closest fit to systems of points in space. Phil. Mag. 2 , 559–572. (doi:x.1080/14786440109462720) Crossref, Google Scholar

• 2
  Hotelling H

  . 1933 Analysis of a complex of statistical variables into principal components. J. Educ. Psychol. 24 , 417–441, 498–520. (doi:10.1037/h0071325) Crossref, Google Scholar

• 3
  Jackson JE

  . 1991 A user's guide to principal components . New York, NY: Wiley. Crossref, Google Scholar

• iv
  Jolliffe It

  . 2002 Principal component analysis , 2nd edn. New York, NY: Springer-Verlag. Google Scholar

• v
  Diamantaras KI, Kung SY

  . 1996 Principal component neural networks: theory and applications . New York, NY: Wiley. Google Scholar

• six
  Flury B

  . 1988 Mutual principal components and related models . New York, NY: Wiley. Google Scholar

• vii
  Horn R, Johnson C

  . 1985 Matrix analysis . Cambridge, UK: Cambridge University Press. Crossref, Google Scholar

• 8
  Hudlet R, Johnson RA

  . 1982 An extension of some optimal properties of principal components. Ann. Inst. Statist. Math. 34 , 105–110. (doi:ten.1007/BF02481011) Crossref, ISI, Google Scholar

• ix
  Okamoto 1000

  . 1969Optimality of main components. In Multivariate analysis Two (ed. PR Krishnaiah), pp. 673–685. New York, NY: Academic Printing. Google Scholar

• x
  McCabe GP

  . 1984 Main variables. Technometrics 26 , 137–144. (doi:10.1080/00401706.1984.10487939) Crossref, ISI, Google Scholar

• eleven
  Cadima J, Cerdeira JO, Minhoto Grand

  . 2004 Computational aspects of algorithms for variable selection in the context of primary components. Comp. Stat. Information Anal. 47 , 225–236. (doi:10.1016/j.csda.2003.11.001) Crossref, ISI, Google Scholar

• 12
  Gill PG, Purnell MA, Crumpton N, Brown KR, Gostling NJ, Stampanoni Grand, Rayfield EJ

  . 2014 Dietary specializations and diversity in feeding ecology of the earliest stem mammals. Nature 512 , 303–305. (doi:10.1038/nature13622) Crossref, PubMed, ISI, Google Scholar

• 13

  R Development Core Team. 2015 R: A language and environment for statistical computing. Vienna, Republic of austria: R Foundation for Statistical Computing. See http://www.R-project.org. Google Scholar

• fourteen
  Gabriel KR

  . 1971 The biplot graphical display of matrices with application to principal component analysis. Biometrika 58 , 453–467. (doi:10.1093/biomet/58.3.453) Crossref, ISI, Google Scholar

• xv
  Cadima J, Jolliffe IT

  . 2009 On relationships betwixt uncentred and column-centred principal component analysis. Pak. J. Stat. 25 , 473–503. Google Scholar

• 16
  Ringner M

  . 2008 What is principal component assay? Nat. Biotechnol. 26 , 303–304. (doi:10.1038/nbt0308-303) Crossref, PubMed, ISI, Google Scholar

• 17
  Lee D, Lee W, Lee Y, Pawitan Y

  . 2010 Super-sparse master component analyses for loftier-throughput genomic information. BMC Bioinform. 11 , 296. (doi:10.1186/1471-2105-eleven-296) Crossref, PubMed, ISI, Google Scholar

• 18
  Birnbaum A, Johnstone IM, Nadler B, Paul D

  . 2013 Minimax bounds for sparse PCA with noisy loftier-dimensional data. Ann. Stat. 41 , 1055–1084. (doi:x.1214/12-AOS1014) Crossref, PubMed, ISI, Google Scholar

• 19
  Johnstone IM, Lu AY

  . 2009 On consistency and sparsity for principal components analysis in high dimensions. J. Am. Stat. Assoc. 104 , 682–693. (doi:x.1198/jasa.2009.0121) Crossref, PubMed, ISI, Google Scholar

• twenty
  Ramsay JO, Silverman BW

  . 2006 Functional data analysis , second edn. Springer Series in Statistics. New York, NY: Springer. Google Scholar

• 21
  Rao CR

  . 1958 Some statistical methods for comparison of growth curves. Biometrics xiv , 1–17. (doi:x.2307/2527726) Crossref, ISI, Google Scholar

• 22
  Li Y, Wang Northward, Carroll RJ

  . 2013 Selecting the number of main components in functional information. J. Am. Stat. Assoc. 108 , 1284–1294. (doi:ten.1080/01621459.2013.788980) Crossref, ISI, Google Scholar

• 23
  Boente G, Silibian-Barrera M

  . 2015 S-estimators for functional primary components. J. Am. Stat. Assoc. 110 , 1100–1111. doi:x.1080/01621459.2014.946991) Crossref, ISI, Google Scholar

• 24
  Hörmann S, Kidziński L, Hallin Yard

  . 2015 Dynamic functional master components. J. R. Stat. Soc. B 77 , 319–348. (doi:10.1111/rssb.12076) Crossref, Google Scholar

• 25
  Brillinger DR

  . 1981 Time serial: data analysis and theory , Expanded edn. San Francisco, CA: Holden-24-hour interval. Google Scholar

• 26
  Cattell RB

  . 1978 The scientific utilise of factor analysis in behavioral and life sciences . New York, NY: Plenum Printing. Crossref, Google Scholar

• 27
  Jolliffe It, Trendafilov N, Uddin M

  . 2003 A modified chief component technique based on the LASSO. J. Comput. Graph. Stat. 12 , 531–547. (doi:x.1198/1061860032148) Crossref, ISI, Google Scholar

• 28
  Hastie T, Tibshirani R, Wainwright K

  . 2015 Statistical learning with sparsity: the LASSO and generalizations . Boca Raton, FL: CRC Printing. Crossref, Google Scholar

• 29
  Zou H, Hastie T, Tibshirani R

  . 2006 Sparse principal components. J. Comput. Graph. Stat. 15 , 262–264. (doi:10.1198/jcgs.2006.s7) Crossref, ISI, Google Scholar

• 30
  Witten D, Tibshirani R, Hastie T

  . 2009 A penalized matrix decomposition, with applications to sparse chief components and canonical correlation analysis. Biostatistics ten , 515–534. (doi:x.1093/biostatistics/kxp008) Crossref, PubMed, ISI, Google Scholar

• 31
  d'Aspremont A, El Ghaoui L, Jordan MI, Lanckriet GRG

  . 2007 A straight formulation for sparse PCA using semidefinite programming. SIAM Rev. 49 , 434–448. (doi:ten.1137/050645506) Crossref, ISI, Google Scholar

• 32
  Lei J, Vu VQ

  . 2015 Sparsistency and agnostic inference in sparse PCA. Ann. Stat. 43 , 299–322. (doi:10.1214/14-AOS1273) Crossref, ISI, Google Scholar

• 33
  Obukhov AM

  . 1947 Statistically homogeneous fields on a sphere. Usp. Mat. Navk. 2 , 196–198. Google Scholar

• 34
  Lorenz EN

  . 1956Empirical orthogonal functions and statistical weather prediction. Technical report, Statistical Forecast Project Report 1, Dept. of Shooting star. MIT: 49. Google Scholar

• 35
  Preisendorfer RW, Mobley CD

  . 1988 Principal component assay in meteorology and oceanography . Amsterdam, The Netherlands: Elsevier. Google Scholar

• 36
  Hannachi A, Jolliffe IT, Stephenson DB, Trendafilov North

  . 2006 In search of simple structures in climate: simplifying EOFs. Int. J. Climatol. 26 , vii–28. (doi:ten.1002/joc.1243) Crossref, ISI, Google Scholar

• 37
  Huber PJ

  . 1977 Robust statistical procedures . Philadelphia, PA: Social club for Industrial and Applied Mathematics. Google Scholar

• 38
  Huber PJ

  . 1981 Robust statistics . New York, NY: Wiley. Crossref, Google Scholar

• 39
  Ruymagaart FH

  . 1981 A robust chief component analysis. J. Multivariate Anal. 11 , 485–497. (doi:10.1016/0047-259X(81)90091-9) Crossref, ISI, Google Scholar

• 40
  Hallin Thousand, Paindaveine D, Verdebout T

  . 2014 Efficient R-interpretation of master and mutual principal components. J. Am. Stat. Assoc. 109 , 1071–1083. (doi:x.1080/01621459.2014.880057) Crossref, ISI, Google Scholar

• 41
  Huber PJ, Ronchetti EM

  . 2009 Robust statistics , 2nd edn. Wiley Series in Probability and Statistics. New York, NY: Wiley. Crossref, Google Scholar

• 42
  De la Torre F, Black MJ

  . 2003 A framework for robust subspace learning. Int. J. Comput. Vis. 54 , 117–142. (doi:10.1023/A:1023709501986) Crossref, ISI, Google Scholar

• 43
  Wright J, Peng Y, Ma Y, Ganesh A, Rao S

  . 2009Robust main component analysis: exact recovery of corrupted low-rank matrices past convex optimization. In Proc. of Neural Data Processing Systems 2009 (NIPS 2009), Vancouver, BC, Canada, 7–x December 2009. See http://papers.nips.cc/paper/3704-robust-principal-component- analysis-exact-recovery-of-corrupted-low-rank-matrices-via-convex-optimization.pdf . Google Scholar

• 44
  Candès EJ, Li X, Ma Y, Wright J

  . 2011 Robust chief component assay? J. ACM 58 , 11:1–eleven:37. Crossref, ISI, Google Scholar

• 45
  Zhao Q, Meng D, Xu Z, Zuo W, Zhang 50

  . 2014Robust master component analysis with complex dissonance. In Proc. of the 31st Int. Conf. on Machine Learning, Beijing, Red china, 21–26 June 2014. Encounter http://jmlr.org/proceedings/papers/v32/zhao14.pdf. Google Scholar

• 46
  Bouwmans T, Zahzah E

  . 2014Robust PCA via principal component pursuit: a review for a comparative evaluation in video surveillance. Comput. Vis. Prototype Underst. 122 , 22–34. (doi:10.1016/j.cviu.2013.11.009) Google Scholar

• 47
  Bock H-H, Diday Due east

  . 2000 Analysis of symbolic data . Berlin, Germany: Springer. Crossref, Google Scholar

• 48
  Brito P

  . 2014 Symbolic information analysis: another look at the interaction of data mining and statistics. WIREs Data Mining Knowl. Discov. iv , 281–295. (doi:10.1002/widm.1133) Crossref, ISI, Google Scholar

• 49
  Ichino One thousand, Yaguchi H

  . 1994 Generalized Minkowski metrics for mixed feature type information analysis. IEEE Trans. Syst. Man Cybern. 24 , 698–708. (doi:10.1109/21.286391) Crossref, Google Scholar

• 50
  Makosso-Kallyth S

  . In press. Principal axes analysis of symbolic histogram variables. Stat. Anal. Data Mining. (doi:10.1002/sam.11270) Google Scholar

• 51
  Vichi Chiliad, Saporta G

  . 2009 Clustering and disjoint principal component assay. Comp. Stat. Data Anal. 53 , 3194–3208. (doi:10.1016/j.csda.2008.05.028) Crossref, ISI, Google Scholar

riosmervagands.blogspot.com

Source: https://royalsocietypublishing.org/doi/10.1098/rsta.2015.0202

Share this post