Principal Component Analysis and biplots. A Back-to-Basics Comparison of Implementations
Ettore Settanni




 
    
        
       
            
       
       
       
      











HighlightingVerbatimcommandchars=
{}



Shaded









































 




























Principal Component Analysis and biplots. A Back-to-Basics
Comparison of Implementations
    
    Ettore Settanni
   

    Institute for Manufacturing 

    University of Cambridge 

  Cambridge, UK 

   

  
    April 29, 2024
=======================================================================================================================





Principal Component Analysis and biplots are so well-established and
readily implemented that it is just too tempting to give for granted
their internal workings. In this note I get back to basics in comparing
how PCA and biplots are implemented in base-R and contributed R
packages, leveraging an implementation-agnostic understanding of the
computational structure of each technique. I do so with a view to
illustrating discrepancies that users might find elusive, as these arise
from seemingly innocuous computational choices made under the hood. The
proposed evaluation grid elevates aspects that are usually disregarded,
including relationships that should hold if the computational rationale
underpinning each technique is followed correctly. Strikingly, what is
expected from these equivalences rarely follows without caveats from the
output of specific implementations alone.




introduction


§ INTRODUCTION


A food scientist, a bioinformaticist, and a survey researcher walk into
a bar… Jokes aside, what could they possibly have in common? The
answer is probably PCA—Principal Component Analysis.

Many will find the concept familiar, perhaps even trite. Overviews of
this technique leave little room for doubt about its popularity and
relevance across disciplines—see e.g., Knox
(ref-KnoxStevenW.2018Ml:t2018, Ch. 11); Bro and
Smilde (ref-chemometric20142014); Abdi and
Williams (ref-AbdiWilliams2010). Often, PCA is
complemented by a biplot—a closely related, yet distinct technique for
jointly visualizing observations and variables that typically make up a
data matrix (ref-DuToit1986GEDAdu Toit, Steyn, and
Stumpf 1986 Ch. 6; ref-Gower2011Gower, Lubbe, and
LeRoux 2011).

There is no shortage of options at the analysts' fingertip when
performing PCA: Tab.<ref> provides a non-exhaustive list.
Base- alone provides two built-in functions besides several
contributed packages' own implementations. For a hands-on summary of
some of these functions see e.g., Mayor
(ref-MayorEric2015Lpaw2015, Ch. 6); Kumar and Paul
(ref-Kumar2016Mtmw2016, Ch. 4). Several references
on multivariate statistics cover both PCA and biplots in the context of
base  (e.g.,
ref-Venables.2002MaswVenables and Ripley 2002, Ch.
11; ref-EverittBrian2011AitaEveritt and Hothorn
2011, Ch. 3) or dedicated  packages (e.g.,
ref-Gower2011Gower, Lubbe, and LeRoux 2011, Ch. 3;
ref-Pages2014MFAbPagès 2014, Ch. 1).



With PCA and biplots being so readily implemented it is tempting to give
for granted their internal workings and go about it mechanically through
a canned routine of choice. The exception that proves the rule is the
package
https://CRAN.R-project.org/package=LearnPCALearnPCA
(ref-LearnPCAHanson and Harvey 2022). Its rich set
of vignettes is unique in its intent to address the self-directed
learner with a view to unpicking methodological aspects of PCA (but not
biplots) that  users rarely engage with, and whose importance
may be underplayed. It also provides comparative insights into
alternative ways of implementing PCA in base .

Against this backdrop, what seems to be missing is—to the best of my
knowledge—a comparison of how PCA and biplots are implemented in base
 and contributed  packages, and how that compares
with an implementation-agnostic understanding of the computational
structure of each technique. In this note I attempt such comparison with
a view to illustrating discrepancies that users might find elusive, as
these arise from seemingly innocuous computational choices made under
the hood. By getting back to basics in PCA and biplots I elevate aspects
that are usually disregarded.

The remainder is structured as follows. The next section outlines an
implementation-agnostic understanding of the computational building
blocks in PCA and biplots with the aid of an illustrative example. Using
these insights as an evaluation grid, selected implementations are then
compared, pinpointing possible points of departure from what is
reasonably expected. A closing section summarises key practical
implications of these findings.

buildblocks


§ IMPLEMENTATION-AGNOSTIC BUILDING BLOCKS


The seasoned practitioner might scoff at the idea of yet another
overview on PCA and swiftly move past. At the risk of disappointing
advanced readers, this section makes a point of reviewing the
computational building blocks of PCA and biplots—especially how they
come about—as an agnostic stance for comparing the output of specific
implementations.

A useful place to start is the visual intuition behind how PCA
works—i.e., that a set of data-points can be represented in a
lower-dimensional space while preserving relevant information about
them. The analogy with image compression often comes to mind
(ref-PooleDavid2015Poole 2014, 607;
ref-Peng_PCAPeng 2020, Ch. 3). In the simplest
case, the data-points are 2-dimensional and the task at hand is to
determine their orthogonal projections onto an appropriately defined
line, thus obtaining a 1-dimensional representation of those points.
While simplistic, the pedagogic merits of a 2-dimensional visual example
are emphasised by the popularity of on-line resources such as Starmer
(ref-StatQuestPCA2018) and amoeba
(ref-stackexchangePCA2015). This is the case I
will consider throughout this section, by means of an illustrative
numerical example. Similar examples can be found elsewhere e.g., Hanson
and Harvey (ref-LearnPCA2022); Abdi and Williams
(ref-AbdiWilliams2010). The example is restricted
to numerical features, or variables, as PCA does not immediately apply
to categorical or binary features—see Kassambara and Mundt
(ref-PCAalternatives2020)'s taxonomy in the
context of package
https://CRAN.R-project.org/package=factoextrafactoextra.

a-motivating-example


 §.§ A motivating example


Consider the example in Tab.<ref>, and
Fig.<ref>. The first two columns of Tab.<ref>
represent the raw data, which are arranged in a so-called data matrix
𝐗 =  x_ij_n × m consisting of
i=1,...,n observations and j=1,...,m features, or variables. In
the example n=6 and m=2. For each feature j the last two
rows of Tab.<ref> give its mean
x̅_j=1/n∑_ix_ij, and sample variance
s_x_j^2=1/n-1∑_i (x_ij - x̅_j)^2. In
 this is just:  and
, respectively.





The columns of Tab.<ref> denoted as 𝐘
correspond to the centred data matrix, which is also of size
n × m, with generic element y_ij = x_ij - x̅_j.
The operation of centring 𝐗 is often expressed in matrix
notation as
𝐘=𝐗-1/n11^T𝐗,
where 1 is a unit vector of appropriate dimension
(ref-Venables.2002MaswVenables and Ripley 2002,
302; ref-Gower2011Gower, Lubbe, and LeRoux 2011).
In  this is accomplished with:
. The
mean and sample variance for the centred data are also shown in the last
two rows of Tab.<ref>. For each feature j the mean is
characteristically zero:
y̅_j=1/n∑_iy_ij=1/n∑_ix_ij-x̅_j=0.
Yet the sample variance is the same as the raw data's:
s_y_j^2=1/n-1∑_i (y_ij - y̅_j)^2=1/n-1∑_i (x_ij - x̅_j - 0)^2=s_x_j.

Contrary to what many think, so long as the quantities represented in
𝐗 are expressed in the same units, or if they scale well
e.g., through a logarithmic transformation, it is not compulsory to
manipulate the data further (see e.g.,
ref-Venables.2002MaswVenables and Ripley 2002).

Having centred the data, each row of 𝐘 can be interpreted
as the coordinates of a point in the Cartesian plane corresponding to an
observation (row header). The n=6 observations can then be plotted
as shown in Fig.<ref>, left-hand side. The same figure also
shows two lines passing through the origin—what is commonly referred
to as principal components—and the orthogonal projections of
the centred data-points' onto each. For a given point, its projection is
found from intersecting a component, say the red line, and its
perpendicular passing through that point (see e.g.,
ref-intermedAlgebraMarecek 2017, Ch. 3). The
position of an orthogonal projection onto a principal component, being
1-dimensional, is just its distance from the origin—or PCA
score. The sign of a score, however, depends on which quadrant in
which the projection falls with respect to the origin. For the two
principal components in our example, the respective scores are
denoted in Tab.<ref> by vectors 𝐳_∙1
(red line) and 𝐳_∙2 (blue line).

What we see at work in Fig.<ref> is a key tenet of PCA,
albeit streamlined: a set of 2-dimensional data-points is transformed
into a set of 1-dimensional points, which are appropriately arranged
along a line, or principal component. In the context of PCA a set of
data-points is “appropriately” represented in a lower-dimensional
space when most information about its variance is retained. In our
example, the more dispersed the projected data-points along a
1-dimensional line, the better.

But how are the lines and projections in Fig.<ref> arrived
at, so that one can claim they are optimal? This aspect will be
discussed next, and will inform an implementation-agnostic evaluation
grid.

PCAsection


 §.§ Computational aspects of PCA


To understand how the data-points' projections are positioned along a
1-dimensional principal component, one can imagine rotating the
lines shown in Fig.<ref> (left) until the variance of the
projections' distances from the origin—i.e., the variance of the
scores in 𝐳_∙1 and
𝐳_∙1—is maximal. For the first principal
component in our example, or PC1 (blue line), the last row in
Tab.<ref> shows that such variance is 21.3. This equals the
sum of the variances along the projections' coordinates, see columns 5
and 6 of Tab.<ref>. Typically, the variance of a principal
component's scores is expressed as a fraction of the total variance
along the the columns of 𝐘. For PC1 in our example, the
maximal projections' variance achieved accounts for about 96% of the
projected points'.

A well-established approach to obtain the optimal scores is to
seek appropriate linear combinations of the original data-points'
coordinates (e.g., ref-JolliffeI.T.2004PcaJolliffe
2004, Ch. 1 & 3; ref-DuToit1986GEDAdu Toit,
Steyn, and Stumpf 1986, Ch. 9). In our example, the coordinates are
represented by the columns of the centred data matrix 𝐘,
and a linear combination of interest is the matrix-vector product:


    𝐳 =
      𝐘𝐚 = [ a_1 y_11 + a_2 y_12 = z_1;                         ⋮; a_1 y_i1 + a_2 y_i2 = z_i;                         ⋮; a_1 y_n1 + a_2 y_n2 = z_n ]


where 𝐚=[ a_1, a_2 ]^T is the
vector of unknown linear combination weights, commonly referred to as
loadings; and 𝐳 is a vector of unknown principal
component scores that position the projections of the data-points
in 𝐘 along an optimally defined 1-dimensional line or
principal component, as discussed.

The process of finding 𝐳 revolves around the sample
variance of its elements. The data being centred, 𝐳 has
mean z̅=1/n∑_iz_i=0, from which it follows that:


    Var(𝐳) =      1/n-1∑_i (z_i - z̅)^2 
        = 1/n-1(𝐘𝐚 - z̅)^T ·(𝐘𝐚 - z̅) 
        = 1/n-1𝐚^T𝐘^T 𝐘𝐚
        = 𝐚^T𝐒𝐚


where 𝐒=1/n-1𝐘^T𝐘 is the sample
covariance matrix of 𝐘. In  this is simply:
.

Based on the above, one must choose 𝐚 to maximize
Var( 𝐳) while respecting orthogonality
constraints (ref-JolliffeI.T.2004PcaJolliffe 2004,
4; ref-Venables.2002MaswVenables and Ripley 2002,
303; ref-EverittBrian2011AitaEveritt and Hothorn
2011, Ch. 3). Since. from Equation <ref>,
Var( 𝐳) = 𝐚^T𝐒𝐚,
this is typically framed as a mathematical program:


    max_𝐚     𝐚^T𝐒𝐚
    s.t.     𝐚^T𝐚=1


As it turns out, solving the program in Equation <ref> for
𝐚 is equivalent to finding the principal eigenvector of
𝐒. To see how that is so, I examine the first-order
conditions for a constrained maximum using Lagrange multipliers
(ref-JolliffeI.T.2004PcaJolliffe 2004, 4):


    d𝐚^T𝐒𝐚/d𝐚 - d/d𝐚λ( 𝐚^T𝐚 -1) =       0 
    
           2𝐚^T𝐒 - 2λ𝐚^T  =       0 
    ( 𝐒 - λ𝐈)𝐚 =       0


where the vector derivative
d𝐚^T𝐒𝐚/d𝐚=2𝐚^T𝐒
is a special result of the general rule
d𝐚^T𝐒𝐚/d𝐚 = 𝐚^T(𝐒^T + 𝐒)
when 𝐒 is symmetric
(ref-BinmoreDaviesBinmore and Davies 2001, 162),
which is the case for covariance matrices. The form of Equation
<ref> is typical of an eigenvalue problem in linear algebra
(e.g. ref-PooleDavid2015Poole 2014, Ch. 4), and
promptly solved in  using the function ,
which in our example yields the eigenvalues of 𝐒


[]
eigen(S)values

denoted as the vector
λ=[ λ_1 λ_2 ]^T;
and the corresponding eigenvectors, juxtaposed

[]
eigen(S)vectors

denoted by the matrix
𝐕 = [ 𝐚_∙1 𝐚_∙2 ].
The above suggests there are two linear combinations—i.e., as many as
there are features—that could satisfy Equation <ref>.
Jolliffe (ref-JolliffeI.T.2004Pca2004, 5)
demonstrates that the same reasoning leading to Equation
<ref> applies recursively beyond the first principal
component.

If one considers only the linear combination associated with the largest
eigenvalue λ = λ_1 =21.28 and the corresponding
eigenvector
𝐚 = 𝐚_∙1=[ -0.94, -0.34 ]^T
of 𝐒, Equation <ref> yields the scores
vector 𝐳=𝐳_∙1 shown in Tab.
<ref>. Choosing
𝐚 = 𝐚_∙2 =[ 0.34, -0.94 ]^T
yields 𝐳=𝐳_∙2, instead.

These vectors may be juxtaposed to form a scores matrix
𝐙 = [ 𝐳_∙ ,1  𝐳_∙,2 ] = 𝐘𝐕.
Each column of the scores matrix 𝐙 is a
lower-dimensional representation from which one can reconstruct the
2-dimensional data-points in 𝐘
(ref-Gower2011Gower, Lubbe, and LeRoux 2011). For
example, columns 5 and 6 in Tab.<ref> are reconstructed as
𝐘_PC1 = 𝐳_∙1𝐚_∙1^T.
In  this is:
.

Combining Equations <ref> and <ref> brings about an
often overlooked equivalence between the variance along a given
scores vector 𝐳 and the corresponding eigenvalue
λ of 𝐒:


    Var[𝐳]       =  𝐚^T𝐒𝐚
          =  𝐚^Tλ𝐚
          = λ


For example, assuming 𝐳=𝐳_∙1 one verifies
in  that:

[]
Z - Y %*%eigen(S)vectorsall.equal(apply(Z,2,var)[1], eigen(S)values[1])

Although it is simplistic to assume two features, the first two
principal components have well-established uses for visualising the
outputs of PCA, even when there are more than two features.
Traditionally, a 2-dimensional visualisation of PCA is attained by
combining the first two principal components in a plot as shown in
Fig.<ref> (right), where the vectors
𝐳_∙,1 and 𝐳_∙,2 provide the
“transformed” system of coordinates and the original n=6
observations are plotted accordingly. In the presence of more than two
features one would typically retain the first two components with the
caveat that the quality of such an approximation is given by the
relative weight of the two largest eigenvalues of the covariance matrix
𝐒. From Equation <ref>, this is equivalent to the
proportion of variance associated with the first two principal
components relative to the variance of the projected data-points.

SVDsection


 §.§ Alternative approaches: Singular Value
Decomposition


The computational building blocks described so far emphasise the logical
path from seeking a variance-maximising linear combination in Equation
<ref> to solving the eigenvalue problem in Equation
<ref>. In practice, the preferred way to obtain
loadings and scores is by doing a Singular Value
Decomposition (SVD) of the rectangular matrix 𝐘 (e.g.,
ref-JolliffeI.T.2004PcaJolliffe 2004, Ch. 3;
ref-LearnPCAHanson and Harvey 2022):


    𝐘=𝐔𝐃𝐕^T


which is equivalent to the command  in . With
reference to our simple example, Equation <ref> generates
the following standard outputs
(ref-PooleDavid2015Poole 2014, Ch. 7):


  * 
  A matrix
  𝐔=[ 𝐮_∙,1 𝐮_∙,3 ]
  formed with the “left” eigenvectors of 𝐘:

[]
 svd(Y)u
  * 
  A matrix
  𝐕=[ 𝐯_∙1 𝐯_∙2 ]
  formed with the “right” eigenvectors of 𝐘[]
svd(Y)v

  which is analogous to the juxtaposed loadings vectors. Yet
  using 𝐕 to indicate both
  [ 𝐚_∙1 𝐚_∙2 ]
  from Equation <ref> and
  [ 𝐯_∙1 𝐯_∙2 ]
  is an abuse of notation, since the elements' signs may differ.

  * 
  A vector
  ℓ=[ l_1 l_2 ]^T of
  “singular values” of 𝐘, which is typically diagonalised
  to form the matrix
  𝐃 = [ l_1   0;   0 l_2 ] in
  Equation <ref>:

[]
svd(Y)d

It is worth emphasising that the singular values in
ℓ are distinct from, but related to the eigenvalues
in λ (e.g.,
ref-LearnPCAHanson and Harvey 2022;
ref-JolliffeI.T.2004PcaJolliffe 2004, 38). In our
example, such relationship is:


    ℓ      = (n-1)^1/2[ λ_1^1/2 λ_2^1/2 ]
          = (n-1)^1/2[ √(Var[𝐳_∙,1]) √(Var[𝐳_∙,1]) ]
          = (n-1)^1/2[ σ_1 σ_2 ]


where σ_j = s.d.[𝐳_∙ j]; and
λ_j = Var[ 𝐳_∙ j], based on Equation
<ref>. One can verify in  that:

[]
all.equal(sqrt((nrow(Y)-1)*eigen(S)values), sqrt(nrow(Y)-1)*apply(Z,2,sd))

To obtain the scores matrix 𝐙 from an SVD, one
right-multiplies both sides of Equation <ref> by
𝐕:


    𝐘𝐕    = 𝐔𝐃𝐕^T 𝐕
        = 𝐔𝐃
        = 𝐙_SVD


In practice, 𝐙_SVD and
𝐙 = [ 𝐳_∙ ,1  𝐳_∙,2 ],
which is shown in Tab.<ref>, are often used interchangeably.
Yet they are only equivalent in absolute values. In our example, one can
verify that:

[]
 all.equal(abs(unname(U %*% D)), abs(unname(Y %*%eigen(S)vector)))sectionBiplots


 §.§ Principal components biplots


The typical PCA plot in Fig.<ref> (right) may be enriched to
become a biplot. The difference between these concepts might seem
cosmetic: both rely on a “transformed” system of coordinates based on
PCA scores, but differ in scope as to what they seek to
visualise. The former only visualises observations i.e., the rows of the
data-matrix; whereas the latter aims to jointly represent observations
and features i.e., both the rows and the columns of the data-matrix
(ref-JolliffeI.T.2004PcaJolliffe 2004, Ch. 5).

On a deeper level, the task of attaining an overlaid representation of a
plot of features and a plot of observations in a biplot is
underpinned by a distinct computational strategy. Whilst it was
straightforward to derive the plot in Fig.<ref> from
previous computations, a biplot is inherently associated with an
SVD of the centred data-matrix (ref-Gower2011Gower,
Lubbe, and LeRoux 2011;
ref-Venables.2002MaswVenables and Ripley 2002, Ch.
11; ref-DuToit1986GEDAdu Toit, Steyn, and Stumpf
1986, Ch. 6; ref-ggbiplot_vignetteVu and Friendly
2024).

In the context of principal component biplots, the SVD in Equation
<ref> serves as the starting point, but is modified so that
the diagonal matrix of singular values is “split” based on a parameter
0≤α≤1. Such modification is illustrated below in the
context of our simplified example, with just two features:


    𝐘    = 𝐔𝐃𝐕^T 
        = 𝐔𝐋^α𝐋^1-α𝐕^T 
        = [ 𝐮_∙,1 𝐮_∙,2 ][ l_1   0;   0 l_2 ]^α[ l_1   0;   0 l_2 ]^1-α[ 𝐯_∙,1; 𝐯_∙,2 ]^α
        = [ u_11l_1^α u_12l_2^α; u_21l_1^α u_22l_2^α;         ⋮         ⋮; u_n1l_1^α u_n2l_2^α;           ][ v_11l_1^(1-α) v_21l_1^(1-α); v_12l_2^(1-α) v_22l_2^(1-α);               ]
        = 𝐀𝐁^T


In the presence of more than two features, the right-hand side of
Equation <ref> yields a rank-2 approximation of the centred
data-matrix i.e.,


    𝐘≈𝐘_(2) = 𝐀𝐁^T


When α=0 Equation <ref> is referred to as a
principal component biplot, which is of interest here. In our
example with just two features, the matrices 𝐀 and
𝐁 take on the following meaning:


  * 
  The coordinates of the observations along the principal components are
  given by 𝐀=𝐔. Based on Equation
  <ref>, and using interchangeably 𝐙 and
  𝐙_SVD for simplicity, it can be shown that this is
  equivalent to scaling the PCA plot coordinates given by the
  scores matrix:

    𝐀    = 𝐔
        = 𝐙𝐃^-1
        =[ 𝐳_∙,1l_1^-1 𝐳_∙,2l_2^-1 ]

  * 
  Simultaneously, the coordinates of the features are given by the
  matrix product 𝐁^T=𝐃𝐕^T between the
  right eigenvectors and the singular values of an SVD. Yet this result
  can be re-expressed in terms of the loadings and
  scores's St.Dev. Based on Equation <ref>:

    𝐁^T     = 𝐃𝐕^T 
        = [ 𝐯_∙,1 l_1; 𝐯_∙,2 l_2 ]
        = (n-1)^1/2[ 𝐯_∙,1σ_1; 𝐯_∙,2σ_2 ]


In our simple example the matrix 𝐀 of observations
coordinates is equivalent to the matrix 𝐔 of left
eigenvectors previously obtained from Equation <ref>, as
expected:

[]
all.equal(Z_svd %*%diag(1/svd(Y)d), svd(Y)u)

whereas the coordinates for the features given by matrix
𝐁^T are:

[]
D %*%t(V)

For our example, a joint representation of observations and features
based on 𝐀 and 𝐁,which are characteristic of a
principal components biplot, is shown in Fig.<ref>.



What is clear from Equations <ref> and <ref>
is that there is more to Fig.<ref> than just overlaying a
crude plot of loadings over a plot of PCA scores—which
happens frequently in practice. If the observations' and features'
coordinates are jointly computed as shown, one can verify certain
identities between the columns of matrix 𝐁 and the standard
deviation of the corresponding features, or their correlation. Such
properties are often overlooked, but useful for comparative purposes and
will be illustrated next.

biplot-properties-linked-to-features


 §.§ Biplot properties linked to
features


The literature often states that, for a given feature j the length
of the vector formed by its coordinates on a biplot should be equivalent
to the feature's standard deviation (St.Dev.) σ_j (e.g.,
ref-Venables.2002MaswVenables and Ripley 2002,
312; ref-DuToit1986GEDAdu Toit, Steyn, and Stumpf
1986, 108). This statement is hardly self-evident and might generate
some confusion. To demonstrate it, I follow Jolliffe
(ref-JolliffeI.T.2004Pca2004, 77) and
left-multiply both sides of Equation <ref> by
𝐘^T:


    𝐘^T𝐘    = 𝐘^T𝐀𝐁^T 
        = 𝐁(𝐔𝐋^0 )^T(𝐔𝐋^0 )𝐁^T 
        = 𝐁𝐔^T𝐔𝐁^T 
    
      (n-1)𝐒    = 𝐁𝐁^T


where 𝐀=𝐔𝐋^α from Equation
<ref>, and α=0 is characteristic of a principal
components biplot. Also, the columns of 𝐔 are orthonormal
hence 𝐔^T𝐔=𝐈. As before, 𝐒
is the covariance matrix of 𝐘.

Next I focus on a generic diagonal element s_jj of the covariance
matrix 𝐒 on the left-hand side of Equation
<ref> to show how the St.Dev. σ^feat_j
of the j-th feature relates to the length of vector
𝐛_j ∙, which corresponds to the j-th column of
matrix 𝐁^T in Equation <ref>:


    (n-1)s_jj      = 𝐛_j ∙^T 𝐛_j ∙
    σ^feat_j       = (n-1)^-1/2‖𝐛_j ∙‖_2 
          = √(1/n-1∑_k( v_jkl_k)^2)
          = √(∑_k (v_jk√(λ_k))^2)


where l_j = ( n-1)^1/2√(λ_j) based on
Equation <ref>. For our example, one can verify in 
that:

[]
all.equal(sqrt(diag(S)), apply(t(B),2,function(x){norm(x,"2")})/sqrt(nrow(Y)-1))

Confusions may arise due to the form of the right-hand side in Equation
<ref>. For instance, it would be imprecise to conclude—as
the phrasing in some references might suggest—that the vector whose
length equals the St.Dev. of a given feature j is the corresponding
column of the matrix 𝐁^T in Equation <ref>. Also,
it would be misleading to expect that Equation <ref> holds for
a rank-2 approximation when there are more than two features, see
Equation <ref>.

A second property associated with principal components biplots is that
the cosine similarity between a pair of vectors of features coordinates
is equivalent to the correlation coefficient between these features,
assuming these have been centred
(ref-JolliffeI.T.2004PcaJolliffe 2004, 77;
ref-DuToit1986GEDAdu Toit, Steyn, and Stumpf 1986,
108). In our example, the presence of two features simplifies the task
of illustrating this point.

The two features vectors shown in Fig.<ref>
(purple-coloured) correspond to the columns of 𝐁 in
Equation <ref>. I denote these as 𝐛_∙1 and
𝐛_∙2. Squaring the 2-norm of their difference
‖𝐛_∙1 - 𝐛_∙2‖_2^2
leads to the following “textbook” equivalence—for details see
Binmore and Davies (ref-BinmoreDavies2001, 18):


    cosθ      = ⟨𝐛_∙1, 𝐛_∙2⟩/‖𝐛_∙1‖_2 ·‖𝐛_∙2‖_2


which is the definition of cosine similarity. In Equation <ref>
the angle in radians between two features vectors is θ; and
⟨·,·⟩ denotes the inner product between them.
What is rarely shown is that Equation <ref> is equivalent to the
correlation between the two features in our simple example. This
equivalence can be proved starting from the correlation coefficient
between two columns of matrix 𝐘:


    corr_𝐲_∙1,𝐲_∙2 = ⟨𝐲_∙ 1, 𝐲_∙ 2⟩/‖𝐲_∙1‖_2 ·‖𝐲_∙2‖_2


Recalling Equations <ref> and <ref>, I denote
the i-th row of matrix 𝐀 as 𝐮_i∙
and rephrase the inner product in Equation <ref> as:


    ⟨𝐲_∙ 1, 𝐲_∙ 2⟩      =   ⟨[ 𝐮_1∙·𝐛_∙1         … 𝐮_i∙·𝐛_∙1         … 𝐚_n∙·𝐛_∙1 ],
           [ 𝐮_1∙·𝐛_∙2         … 𝐮_i∙·𝐛_∙2         … 𝐮_n∙·𝐛_∙2 ]⟩
          =(𝐔𝐛_∙1)^T 𝐔𝐛_∙2


The equivalence in Equation <ref> also applies to the
denominator of Equation <ref> since, for a given feature j,
‖𝐲_∙ j‖_2 = √(⟨𝐲_∙ j, 𝐲_∙ j⟩).
Hence, substituting Equation <ref> in Equation <ref>
yields:


    ⟨𝐲_∙ 1, 𝐲_∙ 2⟩/‖𝐲_∙1‖_2 ·‖𝐲_∙2‖_2      = 𝐛_∙1^T𝐔^T 𝐔𝐛_∙2/√(𝐛_∙ 1^T𝐔^T 𝐔𝐛_∙ 1)√(𝐛_∙ 2^T𝐔^T 𝐔𝐛_∙ 2)
          = 𝐛_∙1^T𝐛_∙2/√(𝐛_∙ 1^T𝐛_∙ 1)√(𝐛_∙ 2^T𝐛_∙ 2)
          = ⟨𝐛_∙ 1, 𝐛_∙ 2⟩/‖𝐛_∙1‖_2 ·‖𝐛_∙2‖_2
    corr_𝐲_∙1,𝐲_∙2      =  cosθ


where 𝐔^T 𝐔=𝐈 is a property of the
“left” eigenvectors matrix from Equation <ref>. In our
example, one can verify that the cosine similarity between
𝐛_∙1 and 𝐛_∙2 is
cosθ=0.84, which indeed equivalent to the correlation
coefficient between 𝐲_∙1 and
𝐲_∙2:

[]
  cos_theta_ab - (B[1,] %*% B[2,]) / (norm(B[1,], "2")*norm(B[2,], "2"))all.equal(as.numeric(cos_theta_ab), as.numeric(cor(Y[,1], Y[,2])))discussion-of-comparative-insights


§ DISCUSSION OF COMPARATIVE
INSIGHTS


This section develops the computational building blocks summarised
earlier, into an implementation-agnostic evaluation grid, and discusses
how specific implementations compare within such grid. The
implementations examined here are those listed in Tab.<ref>,
with key comparative insights summarised upfront in
Tab.<ref> and Tab.<ref> below. Within these
tables, the proposed evaluating grid informs the chosen column headers.
The following subsections detail how these insights came about.



To provide a benchmark, Tab.<ref> extends the illustrative
example to include more than two features. Fig.<ref> shows
the corresponding biplot for a rank-2 approximation assuming that the
first two components are retained. The tables show that specific
functions may share the same algebraic rationale, yet differ in how the
relevant computations are carried out for each building block.

The grid includes verification equivalences corresponding to Equations
<ref>, <ref>, and <ref> i.e.,
relationships that should hold if the computational rationale
underpinning each technique—illustrated in previous sections—is
followed correctly. Strikingly, these equivalences rarely follow without
caveats from the output of specific implementations alone.

group-1-svd


 §.§ Group 1: SVD


The first group of  implementations considered here has a
common trait: the function  is at work under the hood,
providing what is needed to obtain a scores matrix, as per
Equation <ref>, as well as the coordinates for both
observations and features in a principal component biplot, as per
Equation <ref>. This is specifically the case for
Base-, functions like  and ,
although some caveats apply in the latter case.

For example, the command
 applies
 under the hood to the centred data matrix 𝐘.
The returned output  corresponds to the PCA scores in
columns 7-8 in Tab.<ref>, whereas , returns
the matrix of right eigenvectors or loadings.

One aspect that stands out from Tab.<ref> is that, unlike
other implementations in this group,  satisfies
Equation <ref>, which verifies the relationship between the
eigenvalues of the covariance matrix 𝐒, the St.Dev. of the
scores, and the singular values from an SVD:

[]
all.equal(test_pca1sdev*sqrt(nrow(Y)-1), svd(Y)d)

As a corollary, the function also satisfies Equation <ref>:

[]
all.equal(as.numeric(test_pca1sdev2), as.numeric(apply(test_pca1x,2,var)))

The script  reveals how this is achieved:
the output  returns the singular values from the SVD
divided by √(n-1), consistently with Equation <ref>:

[]
all.equal(svd(Y)d/sqrt(nrow(Y)-1), test_pca1sdev)

In base-, a principal component biplot can be obtained from
 via the  function. Although users are
presented with a graphical output,  does not return the
jointly computed coordinates for features and observations. The
function's script  reveals good alignment
with Equations <ref> and <ref>, except for an
imprecision in how the singular values are arrived at—which explains
the caveats for this function in Tab.<ref>.

Specifically,  computes the singular values from
the output  of a  object. How it does so
depends on the argument : the default setting is
, with  denoting that a principal component
biplot is desired. Contrary to what one may expect, when
 the function sets the singular values equal to
, disagreeing with Equation <ref>, which in turn
affects Equations <ref>, <ref> and
<ref>. Surprisingly, when  the
singular values are computed in a fashion similar to Equation
<ref>, although it is erroneously assumed that 
returns the singular values divided by √(n) instead of
√(n-1), the latter being the case for .

Due to these assumptions, when the argument  is
(contestably)  the singular values returned by
 are close to, but different from
ℓ from Equation <ref>:

[]
lam - test_pca1sdev*sqrt(nrow(test_pca1x))all.equal(lam, svd(Y)d)

With the singular values thus obtained,  computes
the matrix of observations coordinates 𝐀 from the scores
matrix consistently with Equation <ref>. Yet the
result differs from what Equation <ref> suggests, and
hence from what is shown in Tab.<ref>and
Fig.<ref>:

[]
all.equal(unname(sweep(test_pca1x,2,lam,"/")), svd(Y)u)

The caveat in Tab.<ref> acknowledges that the equivalence
𝐀 = 𝐔 could be attained with minor changes in how
singular values are computed by , namely:

[]
sv - test_pca1sdev*sqrt(nrow(Y)-1)all.equal(unname(sweep(test_pca1x,2,sv,"/")), unname(svd(Y)u))

Similar caveats apply to the features' coordinates. In principle, the
underlying computations in  are consistent with
Equation <ref>, which yields matrix 𝐁, except that
these results are also affected by how singular values are computed. In
our example, the features' coordinates produced by  are
different from matrix 𝐁 in Equation <ref>—which
is displayed in Fig.<ref>:

[]
B -t(D %*%t(V))all.equal(unname(sweep(test_pca1rotation,2,lam,"*")), B)

Just like the base- functions examined so far, several curated
packages rely on conventional SVD, too. In particular,
https://bioconductor.org/packages/3.18/PcaMethodsPcaMethods
(ref-pcaMethods_paperStacklies et al. 2007) and
https://bioconductor.org/packages/3.18/PCAtoolsPCAtools
(ref-PCAtoolsBlighe and Lun 2023) extend
traditional PCA methods to handle computational challenges that are
specific to bioinformatics data. While both packages offer a homonym
function , they differ in the extent to which they rely on
base- function. For the implementation
, reliance on  is explicit
upon examination of the function's script, and promptly ascertained. For
example:

[]
test_pca2 - pcaMethods::pca(Y, method ="svd")all.equal(test_pca2@scores, test_pca1x[,1:2] )

So long as the selected method in  is
 the equivalence with  can be verified in a
similar fashion for  and .

With regards to , the analogy with
 is less straightforward to identify as it relies on
another package,
https://bioconductor.org/packages/3.18/BiocSingularBiocSingular
, to carry out an SVD. Yet the results associated with a so-called
“exact” SVD are, once again, equivalent to . For
example:

[]
test_pca3 - PCAtools::pca(X, transposed = T)all.equal(as.matrix(test_pca3rotated), test_pca1x)

The equivalence with  can be verified in a similar
fashion for  and .

Unlike their base- counterpart, the biplot functions in
packages
https://bioconductor.org/packages/3.18/pcaMethodspcaMethods
and
https://bioconductor.org/packages/3.18/PCAtoolsPCAtools
do not strictly speaking generate a principal components biplot as per
Equation <ref>. The function 
displays two separate plots of scores and loadings whose
coordinates are derived, respectively, from he the outputs 
and  of the function  without further
processing. Despite its name, the function ,
too, overlays a scores and a loadings plot.

The packages
https://CRAN.R-project.org/package=ggbiplotggbiplot
(ref-ggbiplot_vignetteVu and Friendly 2024) is
dedicated to biplots. Like  in base-, also
 assumes that an SVD-based PCA be carried out
separately e.g., by . The documentation acknowledges
that PCA implementations may differ in how they go about an SVD, and
that  seeks to provide a unifying interface. Yet
 suffers from some imprecisions already encountered
for  e.g., it assumes that the output
 of a  object is equivalent to the
singular values divided by √(n) instead of √(n-1). For
the singular values thus defined, both  and
 return the same observations' coordinates:

[]
test_pca4 - ggbiplot::get_SVD(test_pca1)all.equal(test_pca4U[,1:2], sweep(test_pca1x,2,lam,"/")[,1:2])

If Equation <ref> was followed, instead, the output
 from  would represent the biplot
coordinates for the observations, as well as the matrix of left
eigenvalues from an SVD. Instead, such output is further processed by
the function  as follows: if the desired output is a
principal component biplot i.e., if the argument  is
true,  multiplies  by  and,
continuing with assumpion the input is a  object, by
√(n-1) as well. Ultimately, this process reverts back to the PCA
scores, multiplied by the ratio √(n-1)/ √(n), hence the
caveat in Tab.<ref>:

[]
df.u -as.data.frame(sweep(test_pca4U[, 1:2], 2, test_pca4D[1:2], FUN ="*"))df.u - df.u *sqrt(nrow(test_pca1x) -1)all.equal(test_pca1x[,1:2] *sqrt(nrow(Y)-1) /sqrt(nrow(Y)), as.matrix(df.u))

Moving on to the features' coordinates,  operates
similarly to Equation <ref>, except that the scores' St.Dev.
 from  are used instead of the singular
values 𝐃. Hence the results not only disagree with Equation
<ref>, but also with :

[]
df.v -as.data.frame(sweep(test_pca1rotation[, 1:2], 2, test_pca4D[1:2], FUN ="*"))all.equal(sweep(test_pca1rotation[,1:2],2,lam[1:2],"*"), as.matrix(df.v))group-2-eigenvalue-problem


 §.§ Group 2: eigenvalue
problem


The second group of implementations considered here follows more closely
the variance optimisation rationale underlying PCA, and the
eigendecomposition problem in Equation <ref> concerning the
covariance matrix 𝐒. Whilst less prevalent than SVD,
 implementations adopting this approach are commonly adopted
in the literature—e.g., Mayor
(ref-MayorEric2015Lpaw2015); Hanson and Harvey
(ref-LearnPCA2022). Yet there may be discrepancies
between our agnostic evaluation grid for both base- and
curated packages.

A natural starting point is the function  in
base-. Nearly a homonym of , this
implementation relies upon the eigenvalues and eigenvectors of the
covariance matrix 𝐒, rather than the outputs of a
data-matrix SVD. Yet  disagrees with Equation
<ref> in that it uses the number of observations n, not
n-1, as the divisor for the covariance matrix. What is more, the
signs of the scores returned by  differ from those
computed directly using  due to further processing within
the function:

[]
test_pca5 -princomp(X, fix_sign = F)test_pca5_eig -eigen(1/nrow(Y)*t(Y) %*% Y, symmetric =TRUE)vectorsall.equal(abs(unname(Y %*% test_pca5_eig)), abs(unname(test_pca5scores)))

Hence, the scores returned by  differ in sign from
𝐙 shown in Tab.<ref> despite both being
obtained from .

When PCA is implemented using  instead of
, the eigenvalues do not equal the variance of the
scores as in Equation <ref>:

[]
 all.equal(as.numeric(test_pca5sdev2), as.numeric(apply(test_pca5scores,2,var)))

An object form  can serve as an input to the
functions  and  examined
earlier. Hence, the same caveats discussed before apply, with the
difference that now both functions are correct in assuming that the the
singular values from an SVD are equivalent to multiplying the output
 of  times √(n):

[]
all.equal(as.numeric(svd(Y)d), as.numeric(test_pca5sdev*sqrt(nrow(Y))))

An option which is very close to  in base-
is the function  from the package
https://CRAN.R-project.org/package=ade4ade4(ref-ade4bookThioulouse
et al. 2018). Both functions define the covariance matrix of the
centred data-points in similar terms, thus generating comparable outputs
for PCA scores

[]
test_pca6 - ade4::dudi.pca(X, center =TRUE, scale =FALSE, scannf =FALSE, nf =ncol(Y))all.equal(unname(as.matrix(test_pca6li)), unname(as.matrix(test_pca5scores)))

Similar equivalences hold for the eigenvectors (loadings) and
eigenvalues of covariance matrix. The implementation in
, too, disagrees with Equation <ref>:

[]
 all.equal(as.numeric(test_pca6eig), as.numeric(apply(test_pca6li,2,var)))

Package https://CRAN.R-project.org/package=amapamap
(ref-amapLucas 2022) provides a different take on
the covariance matrix with its function , which
computes the eigenvalues and eigenvectors of 𝐘^T𝐘,
unlike the previous options:

[]
test_pca7 - amap:::acp(X, reduce =FALSE)test_pca7_eig -eigen(t(Y) %*% Y, symmetric =FALSE)all.equal(unname(test_pca7_eigvectors), unname(test_pca7loadings))

Confusingly, the function's output  returns the square roots
of the eigenvalues:

[]
all.equal(sqrt(test_pca7_eigvalues), test_pca7eig)

In absolute values, the scores and loadings generated by 
are comparable to those obtained in the previous cases, as well as with
those in Tab.<ref>. Oddly, the output  of the
function  is computed directly as the St.Dev. of such
scores:

[]
all.equal(test_pca7sdev, apply(test_pca7scores, 2, sd)) 

Yet the equivalence in Equation <ref> is not satisfied:

[]
all.equal(unname(test_pca7_eigvalues), apply(unname(test_pca7scores), 2, var))

Unlike the implementations described so far, the function
 from the package
https://CRAN.R-project.org/package=psychpsych
(ref-psychRevelle, W. 2024) runs an
eigendecomposition of the covariance matrix 𝐒, consistently
with our evaluation grid. Except that the meaning of “loadings” is
peculiarly defined as the eigenvectors rescaled by the square root of
the eigenvalues. This is reminiscent of—but not identical to—the
feature's coordinate in a biplot obtained in Equation <ref>.

If  is applied to the centred data matrix, and if
the implicit scaling just mentioned is corrected, then the scores are
equivalent, in absolute values, with those shown in
Tab.<ref>:

[]
test_pca8 - psych::principal(Y, nfactors =0, cor="cov", rotate ="none")test_pca8_adjscor -sweep(unname(test_pca8scores), 2, sqrt(test_pca8values), "*")all.equal(abs(test_pca8_adjscor), abs(unname(Z)))

Similar considerations apply to the outputs  and
. The latter returns the covariance matrix's eigenvalues,
and satisfies the equivalence in Equation <ref> provided that
the scores are once again corrected from their implicit scaling:

[]
all.equal(unname(test_pca8values), apply(test_pca8_adjscor, 2, var))

All three packages considered in this section offer visualisation
facilities for generating a combined plot of the PCA scores and
loadings. Examples include  for
 objects; , for 
objects; and  for  objects. Yet
none of these functions seems to jointly compute the features' and
observations' coordinates according to Equation <ref>,
hence one could question whether their visual output is a principal
component biplot.

What is more, objects generated by  can also be an
input to the function  discussed in the previous
section. A word of caution, however, is necessary as this biplot
implementation may not agree with Equation <ref> due to how
 chooses the equivalent of the left and right
eigenvector matrices needed to generate a biplot—in addition to the
caveats already mentioned.

group-3-generalised-svd


 §.§ Group 3: generalised SVD


The final group of implementations shares an intent to move beyond some
of the restrictions concerning what type of data can be handled by PCA,
allowing for the inclusion, for example, of “mixed” data. Another
commonality is that the underlying computational device is the so-called
“generalised” SVD (ref-AbdiWilliamsAbdi and
Williams 2010, Appendix B). None of these characteristics seem to have
equivalent implementations in base-.

The most prominent package in this group is probably
https://CRAN.R-project.org/package=FactoMineRFactoMineR
(ref-factominerPaperLê, Josse, and Husson 2008),
whose functions  is of particular interest here (e.g.,
ref-Pages2014MFAbPagès 2014, Ch. 1). This
implementation is underpinned by a generalised SVD function
, which also operates under the hood of packages
such as
https://CRAN.R-project.org/package=PCAmixdataPCAmixdata
(ref-chavent2022multivariateChavent et al. 2022).

As in previous cases, the PCA scores are equivalent, in absolute values,
to those in Tab.<ref>. Since the results of the generalised
SVD are available as part of the outputs of , users can
directly verify that the PCA scores are obtained consistently with
Equation <ref>:

[]
 test_pca9 - FactoMineR::PCA(X, scale.unit =FALSE)[]
D -diag(test_pca9svdvs)all.equal(test_pca9indcoord, test_pca9svdU %*% D)

The singular values in this generalised SVD are related to the
eigenvalues of the covariance matrix whose elements are divided by
√(n) instead of √(n-1), which we encountered previously:

[]
all.equal(test_pca9eig, test_pca9svdvs2)[]
all.equal(test_pca9eig, unname(eigen(1/nrow(Y)*t(Y) %*% Y)values))

However, this implementation too disagrees with the equivalence between
the scores' variance and covariance matrix's eigenvalues in Equation
<ref>:

[]
all.equal(apply(test_pca9indcoord,2,var) , test_pca9eig)

Moving on to the biplot implementations in this group, it is worth
noting that  generates coordinates for both features and
observations and it automatically produces a plot for each when called.
The observations features are simply the scores, whereas the features'
coordinates correspond to matrix 𝐁 in Equation
<ref>:

[]
all.equal(test_pca9varcoord, test_pca9svdV %*% D)

Based on the above it seems reasonable to conclude that visualisation
devices like  and
 do not generate, strictly speaking, a
principal component biplot. The package
https://CRAN.R-project.org/package=factoextrafactoextra
(ref-PCAalternativesKassambara and Mundt 2020)
provides additional visualisation facilities building on the output
generated by , such as the function
. Yet I could not ascertain an immediate
connection with Equation <ref>.

It is worth noting that  and 
hint to the relationship between how features correlate and the angles
formed by the vectors that represent their coordinates on a biplot (see
e.g., ref-Pages2014MFAbPagès 2014, Ch. 1;
ref-STHDA_PCAKassambara 2017). Despite the
terminology used, these metrics bear little resemblance with the
properties concerning the the cosine of the angle between features'
vectors on a principal components biplot, or their lengths, from
Equations <ref> and <ref>.

concluding-remarks


§ CONCLUDING REMARKS


PCA and biplots are well-established techniques, and when it comes to
resources illustrating how they are done in practice,  users
are spoiled for choice. Whilst most resources hint to the underpinning
matrix decomposition approaches, it is not customary to dwell on how key
computational building blocks come about. Users are often left with a
rather procedural understanding of PCA and biplots tied to specific
implementations. Also, the overwhelming prevalence of SVD-based
approaches can obscure the rationale with which one arrives at such
concepts as scores and loadings.

In a context where implementing PCA and biplots has become somewhat a
conditioned reflex, one can hardly resist the temptation to give for
granted the internal workings. Yet in this note I have taken a
back-to-basics approach to comparing PCA and biplots implementations in
base  and a selection of contributed  packages. This
approach highlighted useful equivalences that should hold if the
computational rationale underpinning each technique is followed
correctly.

Findings suggest that discrepancies from an implementation-agnostic
understanding of PCA and biplots do arise, in both base  and
contributed  packages, from seemingly innocuous computational
choices made under the hood. Surprisingly, the identified verification
equivalences rarely follow without caveats from the output of specific
implementations alone. What is more, biplots are often just a misnomer
due to imprecisions in how the underlying computational aspects are
dealt with.

Base- implementations such as  and
 appear to be most convincingly aligned with what one
expects from the fundamental algebra of PCA and biplots, although some
caveats apply in the case of  that might call for minor
amendments. Despite being based on an SVD, functions like
 adhere more closely to our evaluation grid than
 or . Functions with
richer capabilities and better visualisation devices, which are
understandably more popular amongst user, are harder to reconcile with
an implementation-agnostic understanding of PCA and biplots. This is the
case for packages like
https://CRAN.R-project.org/package=FactoMineRFactoMineR
and https://CRAN.R-project.org/package=ggbiplotggbiplot.

This work has no pretence of comprehensiveness, and only some of the
scripts could be reviewed in an attempt to pinpoint possible
discrepancies from the proposed evaluation grid. Other resources provide
more hands-on comparative outlooks that might speak to practitioners
better than this note does; yet they rarely include biplots or offer an
implementation-agnostic perspective.

Despite its limitations, the hoped for outcome of this note is to raise
awareness that getting back-to-basic in PCA and biplot—two techniques
practically given for granted—helps to elevate aspects that are
usually disregarded for comparative purposes, and to address
discrepancies that users continue to find elusive despite the extensive
resources available for their implementation.

references


§ REFERENCES
tocsectionReferencesrefs10 preref-AbdiWilliams
Abdi, H., and L. J. Williams. 2010. “Principal Component Analysis.”WIREs Computational Statistics 2 (4): 433–59.
https://doi.org/<https://doi.org/10.1002/wics.101>.

 preref-stackexchangePCA
amoeba. 2015. “Answer to: Making Sense of Principal Component
Analysis, Eigenvectors and Eigenvalues.” Cross Validated. 2015.
<https://stats.stackexchange.com/q/140579>.

 preref-BinmoreDavies
Binmore, K., and J. Davies. 2001. Calculus. Cambridge: Cambridge
University Press.

 preref-PCAtools
Blighe, K., and A. Lun. 2023. PCAtools: PCAtools: Everything
Principal Components Analysis.
<https://doi.org/10.18129/B9.bioc.PCAtools>.

 preref-chemometric2014
Bro, R., and A. K. Smilde. 2014. “Principal Component Analysis.”Anal. Methods 6: 2812–31.
<https://doi.org/10.1039/C3AY41907J>.

 preref-chavent2022multivariate
Chavent, M., V. Kuentz-Simonet, A. Labenne, and J. Saracco. 2022.
“Multivariate Analysis of Mixed Data: The r Package PCAmixdata.”<https://arxiv.org/abs/1411.4911>.

 preref-DuToit1986GEDA
du Toit, S. H. C., A. G. W. Steyn, and R. H. Stumpf. 1986.
Graphical Exploratory Data Analysis. Springer.

 preref-EverittBrian2011Aita
Everitt, B., and T. Hothorn. 2011. An Introduction to Applied
Multivariate Analysis with r. Springer.

 preref-Gower2011
Gower, J. C., S. Lubbe, and N. LeRoux. 2011. Understanding
Biplots. Chichester: Wiley.

 preref-LearnPCA
Hanson, Bryan A., and David T. Harvey. 2022. LearnPCA: Functions,
Data Sets and Vignettes to Aid in Learning Principal Components Analysis
(PCA). <https://CRAN.R-project.org/package=LearnPCA>.

 preref-JolliffeI.T.2004Pca
Jolliffe, I. T. 2004. Principal Component Analysis. 2nd ed.
Springer.

 preref-STHDA_PCA
Kassambara, A. 2017. “PCA - Principal Component Analysis
Essentials.” STHDA. 2017.
<http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/#r-packages>.

 preref-PCAalternatives
Kassambara, A., and F. Mundt. 2020. Factoextra : Extract and
Visualize the Results of Multivariate Data Analyses.
<https://cran.r-project.org/web/packages/factoextra/readme/README.html>.

 preref-KnoxStevenW.2018Ml:t
Knox, S. W. 2018. Machine Learning : Topics and Techniques. 1st
ed. Wiley.

 preref-Kumar2016Mtmw
Kumar, A., and A. Paul. 2016. Mastering Text Mining with r. Pakt
Publishing.

 preref-factominerPaper
Lê, S., J. Josse, and F. Husson. 2008. “FactoMineR: A Package for
Multivariate Analysis.”Journal of Statistical Software 25 (1):
1–18. <https://doi.org/10.18637/jss.v025.i01>.

 preref-amap
Lucas, A. 2022. Amap: Another Multidimensional Analysis Package.
<https://CRAN.R-project.org/package=amap>.

 preref-intermedAlgebra
Marecek, L. 2017. Intermediate Algebra. OpenStax.
<https://openstax.org/books/intermediate-algebra/pages/preface>.

 preref-MayorEric2015Lpaw
Mayor, E. 2015. Learning Predictive Analytics with r. Pakt
Publishing.

 preref-Pages2014MFAb
Pagès, J. 2014. Multiple Factor Analysis by Example Using r. 1st
ed. Chapman; Hall/CRC.

 preref-Peng_PCA
Peng, R. D. 2020. “Exploratory Data Analysis with r - SVD and PCA.”
2020.
<https://bookdown.org/rdpeng/exdata/dimension-reduction.html#svd-and-pca>.

 preref-PooleDavid2015
Poole, D. 2014. Linear Algebra : A Modern Introduction. 4th ed.
Brooks Cole.

 preref-psych
Revelle, W. 2024. Psych: Procedures for Psychological,
Psychometric, and Personality Research. Evanston, Illinois:
Northwestern University. <https://CRAN.R-project.org/package=psych>.

 preref-pcaMethods_paper
Stacklies, W., H. Redestig, M. Scholz, D. Walther, and J. Selbig. 2007.
“pcaMethods – a Bioconductor Package Providing PCA Methods for
Incomplete Data.”Bioinformatics 23: 1164–67.
<https://doi.org/10.1093/bioinformatics/btm069>.

 preref-StatQuestPCA
Starmer, J. 2018. “StatQuest: Principal Component Analysis (PCA),
Step-by-Step.”[media] Youtube. 2018.
<https://youtu.be/FgakZw6K1QQ?si=ROWqV5NwFNNKCXVv>.

 preref-ade4book
Thioulouse, J., S. Dray, Dufour A.–B, A. Siberchicot, T. Jombart, and
S. Pavoine. 2018. Multivariate Analysis of Ecological Data with
ade4. Springer.

 preref-Venables.2002Masw
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics
with s. 4th ed. Springer.

 preref-ggbiplot_vignette
Vu, V. Q., and M. Friendly. 2024. Ggbiplot: A Grammar of Graphics
Implementation of Biplots.
<https://CRAN.R-project.org/package=ggbiplot>.