Factor Analysis

Factor analysis was not born in a chemistry laboratory but in psychology. In 1904, the British psychologist Charles Spearman published "'General Intelligence,' Objectively Determined and Measured" [1]. He had noticed that children who did well on one kind of test (say, classics) also tended to do well on apparently unrelated ones (say, discrimination of musical pitch). His explanation was a single hidden quantity, which he called g for "general intelligence", that influenced all the tests at once. The test scores were not the fundamental quantities; they were noisy reflections of something deeper that could not be observed directly. This is the latent variable idea, and it needed new mathematics to make it precise.

In the 1930s and 1940s, the American psychologist L. L. Thurstone expanded that mathematics dramatically. Where Spearman argued for one general factor, Thurstone believed in several distinct abilities (verbal, spatial, numerical and others). He developed multiple factor analysis, introduced factor rotation to reach what he called simple structure, a solution in which each variable depends on as few factors as possible, and turned factor analysis into a general multivariate method [2].

The method entered chemistry through Edmund R. Malinowski, a physical chemist at Stevens Institute of Technology. His book Factor Analysis in Chemistry, written with D. G. Howery and first published in 1980, showed how the same machinery could count the independent chemical components in a set of mixture spectra and help recover their pure spectra [3],[4]. In chemistry, the "latent factors" are not abstract constructs: they are real chemical species whose spectra and concentrations can, at least in principle, be measured on their own. That concrete meaning is what made factor analysis so productive in chemometrics, and it led to a family of methods that includes target factor analysis, evolving factor analysis and self-modelling curve resolution.

The central idea of factor analysis is that the variables you measure are not the variables that matter. A smaller number of hidden factors generate the data, and what you observe is their combined effect plus some noise.

Think of UV-Vis spectra of mixtures of three compounds, measured at 200 wavelengths. You record 200 absorbances per sample, but there are only three independent sources of variation: the three concentrations. The 200 absorbances rise and fall together because they share the same three causes. Factor analysis turns this intuition into a model.

In the factor model, each observed variable $x_j$ is a weighted sum of $m$ common factors $f_1,\dots ,f_m$ plus a term of its own:

The coefficients ${\lambda }_{jk}$ are the factor loadings. They tell how strongly variable $j$ depends on factor $k$. The term ${\epsilon }_j$, the unique factor, collects everything that belongs to variable $j$ alone: measurement noise, a detector artefact, or any variation specific to that wavelength.

This is a generative model: first come the factors, and the data follow from them. Correlations between variables exist because the variables share causes.

Factor analysis and PCA are often confused, and for good reasons: both reduce many variables to a few, and both are often computed from eigenvalues of a covariance or correlation matrix. But they rest on different assumptions and answer different questions.

The practical consequence: if you want to understand the latent structure that produces your data, factor analysis is the natural framework. If you only need a compact summary for later modelling, PCA is simpler and makes fewer assumptions.

For a vector $\mathbf{x}=(x_1,\dots ,x_p{)}^T$ of $p$ observed variables, centred to zero mean, the model is

where:

• $\mathbf{\Lambda }$ is the $p\times m$ loading matrix, linking variables to factors.
• $\mathbf{f}$ is the vector of $m$ common factors, with $m<p$.
• $\mathbf{\epsilon }$ is the vector of $p$ unique factors.

The classical model makes three assumptions:

1. The common factors are uncorrelated and have unit variance: $\text{Cov}(\mathbf{f})=\mathbf{I}$.
2. The unique factors are uncorrelated with each other: $\text{Cov}(\mathbf{\epsilon })=\mathbf{\Psi }$, a diagonal matrix.
3. Common and unique factors are uncorrelated with each other.

Under these assumptions, the covariance matrix of the observed variables splits into two parts:

This is the fundamental equation of factor analysis. $\mathbf{\Lambda }{\mathbf{\Lambda }}^T$ is the covariance produced by the shared factors, and $\mathbf{\Psi }$ adds the variance that belongs to each variable alone. Notice what the equation says: all the covariances between different variables come from the common factors, because $\mathbf{\Psi }$ has nothing off its diagonal.

For each variable, the variance splits in two:

The communality $h_j^2$ is the part of variable $j$'s variance explained by the common factors. The uniqueness ${\psi }_j$ is the rest. For standardised variables, the two add up to 1. In spectroscopy, a wavelength with a high communality is one that responds strongly to the chemical species; a wavelength with a high uniqueness is dominated by noise or by something unrelated to the other wavelengths.

Given a sample covariance (or correlation) matrix $\mathbf{S}$, we look for $\mathbf{\Lambda }$ and $\mathbf{\Psi }$ such that $\mathbf{S}\approx \mathbf{\Lambda }{\mathbf{\Lambda }}^T+\mathbf{\Psi }$. Two methods dominate.

This is the most intuitive approach [5].

The method is simple and widely available, but it does not provide a formal statistical test for the number of factors.

Maximum likelihood factor analysis, introduced by D. N. Lawley in 1940 [6],[7], assumes that the data follow a multivariate normal distribution and looks for the $\mathbf{\Lambda }$ and $\mathbf{\Psi }$ that make the observed data most likely under that model. Up to a constant, it maximises

where $n$ is the number of samples, $|\cdot |$ is the determinant and $\text{tr}$ the trace (the sum of the diagonal). There is no closed-form solution, so it is solved iteratively. Its main advantage is a likelihood ratio test: an approximately chi-squared statistic that tests whether $m$ factors are enough to reproduce the observed covariances. With low-noise spectra, maximum likelihood often drives some uniquenesses to zero (a Heywood case); software then clamps them at a small lower bound and may refuse to compute the test.

The loading matrix is not unique. For any orthogonal matrix $\mathbf{R}$ (a rotation, with $\mathbf{R}{\mathbf{R}}^T=\mathbf{I}$), the loadings ${\mathbf{\Lambda }}^{\ast }=\mathbf{\Lambda }\mathbf{R}$ reproduce the covariance matrix exactly as well, because

(We call the rotation $\mathbf{R}$ to avoid confusion with the scores $\mathbf{T}$ of PCA.)

This rotational indeterminacy means that the first solution an algorithm returns is just one of infinitely many equivalent ones. Rotation uses that freedom to pick the solution that is easiest to interpret.

The unrotated solution tends to produce a first factor that loads on almost everything (a "general" factor) and later factors that contrast groups of variables. It is mathematically natural but hard to read. Rotation redistributes the variance among the factors to approach Thurstone's simple structure: each variable loads strongly on one factor and weakly on the others.

Varimax, introduced by Henry Kaiser in 1958 [8], is the most widely used orthogonal rotation. It maximises the variance of the squared loadings within each factor:

A factor whose squared loadings are all similar has a small variance; a factor with a few large loadings and many near zero has a large one. Maximising $V$ therefore pushes loadings towards either large or near-zero values. (In practice the loadings are often divided by the square root of the communalities before rotation, so that every variable counts equally; this is Kaiser normalisation.) Varimax keeps the factors uncorrelated.

The panel below shows rotation on a simple two-factor example: sixteen wavelengths of a two-component mixture. The data are only mean-centred, as is usual for spectra, so the loadings keep absorbance units and can be compared directly with the bands.

In many real problems, the underlying factors are not independent: in a reaction, for example, the concentration of a product rises as the reactant falls. Promax [9] starts from the varimax solution, raises its loadings to a power (usually 4, keeping the sign) to sharpen the pattern, and then finds the oblique (non-orthogonal) rotation that comes closest to this sharpened target. The result often has simpler structure than varimax, at the price of correlated factors.

In mixture analysis, rotation has a direct physical meaning. The factors obtained from an eigenvalue decomposition of spectral data are abstract: mathematically valid linear combinations of the pure spectra, but not the pure spectra themselves. Rotation can bring them closer to the real spectra, especially when the bands do not overlap too much. Adding chemical constraints, such as non-negative spectra and concentrations, is the basis of target rotation and self-modelling curve resolution, discussed below.

Choosing the number of factors $m$ is one of the most consequential decisions. Several criteria exist, and they do not always agree.

Keep the factors whose eigenvalues of the correlation matrix are larger than 1 [10]. The reasoning is that a factor should explain more variance than a single standardised variable does. The rule is simple and popular, but simulation studies show that it tends to keep too many factors, especially when there are many variables [11].

Plot the eigenvalues in decreasing order and look for the elbow where they change from a steep drop to a flat tail of "rubble" [12]. It is subjective, but often works well when the signal is clearly larger than the noise.

Generate random data sets of the same size as yours (same numbers of samples and variables) but with no factor structure, and compute their eigenvalues. Keep only the factors whose eigenvalues are larger than those of the random data [13]. This accounts for the fact that pure noise also produces eigenvalues larger than zero, and some larger than one. Comparative studies rank parallel analysis among the most accurate rules [11].

With maximum likelihood, fit models with $m=1,2,3,\dots$ factors and choose the smallest $m$ that the chi-squared test does not reject. With very large samples the test rejects even tiny misfits, and with small samples the chi-squared approximation is poor, so use it with judgement.

Leave out part of the data, fit models with different numbers of factors, and check how well the left-out part is reproduced. The number that predicts unseen data best is retained [14]. As with PCA, the left-out values must not be used to compute their own scores.

Chemists brought factor analysis to bear on a concrete problem: given spectra of samples that contain unknown combinations of species, how many species are there, and what are their spectra?

If a set of mixture spectra follows the Beer-Lambert law, the data matrix $\mathbf{D}$ (samples by wavelengths) can be written as concentrations times pure spectra,

and, without noise, its rank (the number of independent rows or columns) equals the number of absorbing species $m$, provided their spectra are linearly independent and their concentrations vary independently. Closure (concentrations that sum to a constant) or species that always appear in fixed ratios lower the rank, and so does mean-centring. With noise, every eigenvalue becomes non-zero, but the first $m$ stay well above a floor of small "noise" eigenvalues. Estimating this chemical rank is the first task of chemical factor analysis.

Malinowski developed several tools for it, including the real error (RE) and the indicator function (IND) [15]. For a matrix with $r$ rows and $c$ columns ($r>c$) and eigenvalues ${\lambda }_j$ of ${\mathbf{D}}^T\mathbf{D}$, they are

RE estimates the size of the noise if $n$ factors were real, and IND typically reaches a minimum at the true number of factors. With the right number of factors, RE comes out close to the true noise level, so it doubles as an estimate of the measurement error. Try it on simulated mixtures:

When spectra are recorded in sequence, as along a chromatographic peak or during a titration or reaction, evolving factor analysis (EFA) computes the eigenvalues of growing windows of the data: the first two spectra, the first three, and so on, and also the same from the end backwards. A new eigenvalue rising out of the noise marks the appearance of a new species; following the backward analysis shows when each one disappears. The result is a "concentration window" for each component, which is valuable information for the next step. EFA was introduced for titrations by Gampp, Maeder and colleagues in 1985 [16] and extended to chromatography by Maeder [17].

The factors from an eigenvalue decomposition are abstract. To reach physically meaningful spectra and concentrations, additional information is needed:

• Target factor analysis (TFA) tests whether a candidate spectrum, for example a library spectrum of a suspected compound, lies in the space spanned by the factors. If it can be reproduced well as a combination of them, it is a plausible component of the mixture [4].
• Self-modelling curve resolution, pioneered by Lawton and Sylvestre in 1971 [18], uses the constraint that spectra and concentrations cannot be negative to find the range of physically acceptable solutions.
• Multivariate curve resolution by alternating least squares (MCR-ALS) generalises this idea with constraints such as non-negativity, unimodality and closure, and is today the most widely used method of this family. It was developed by Romà Tauler and colleagues in the 1990s [19],[20],[21].

The examples simulate spectra of three-component mixtures and fit a three-factor model with varimax rotation. Maximum likelihood (in MATLAB and R) needs more samples than variables, so the simulation uses 200 samples and 30 wavelengths. The noise is set to 0.10: with much less noise the central wavelengths become almost pure signal, maximum likelihood drives many uniquenesses to zero (a Heywood case), and MATLAB then refuses to compute the test.

import numpy as np
from sklearn.decomposition import FactorAnalysis
import matplotlib.pyplot as plt

# --- Simulate a 3-component spectroscopic mixture ---
rng = np.random.default_rng(42)
n_samples, n_wavelengths, n_factors = 200, 30, 3
wavelengths = np.linspace(200, 800, n_wavelengths)
centers = [350, 500, 650]
pure_spectra = np.array([np.exp(-0.5 * ((wavelengths - c) / 60) ** 2) for c in centers])

concentrations = rng.exponential(1.0, (n_samples, n_factors))
spectra = concentrations @ pure_spectra + rng.normal(0, 0.10, (n_samples, n_wavelengths))

# Standardise so that communalities lie between 0 and 1
Z = (spectra - spectra.mean(axis=0)) / spectra.std(axis=0, ddof=1)

# --- Factor analysis with varimax rotation (scikit-learn >= 0.24) ---
fa = FactorAnalysis(n_components=n_factors, rotation='varimax', random_state=0)
fa.fit(Z)
loadings = fa.components_.T                  # (n_wavelengths, n_factors)
communalities = np.sum(loadings ** 2, axis=1)
uniquenesses = fa.noise_variance_

print(f"Mean communality: {communalities.mean():.3f}")
print(f"Mean uniqueness:  {uniquenesses.mean():.3f}")

# Eigenvalues of the correlation matrix (as in PCA), the usual input for Kaiser's rule and the scree test
eigenvalues = np.linalg.eigvalsh(np.corrcoef(Z.T))[::-1]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].plot(range(1, 11), eigenvalues[:10], 'o-')
axes[0].axhline(1, color='r', linestyle='--', label='Kaiser criterion')
axes[0].set_xlabel('Factor number'); axes[0].set_ylabel('Eigenvalue')
axes[0].set_title('Scree plot'); axes[0].legend()

for k in range(n_factors):
    axes[1].plot(wavelengths, loadings[:, k], label=f'Factor {k + 1}')
axes[1].set_xlabel('Wavelength (nm)'); axes[1].set_ylabel('Loading')
axes[1].set_title('Varimax-rotated loadings'); axes[1].legend()

axes[2].plot(wavelengths, communalities, 'g-', linewidth=2)
axes[2].set_ylim(0, 1.05)
axes[2].set_xlabel('Wavelength (nm)'); axes[2].set_ylabel('Communality')
axes[2].set_title('Communalities')

plt.tight_layout()
plt.show()

% --- Simulate a 3-component spectroscopic mixture ---
rng(42);
n_samples = 200; n_wavelengths = 30; n_factors = 3;
wavelengths = linspace(200, 800, n_wavelengths);
centers = [350, 500, 650];
pure_spectra = zeros(n_factors, n_wavelengths);
for i = 1:n_factors
    pure_spectra(i, :) = exp(-0.5 * ((wavelengths - centers(i)) / 60).^2);
end

concentrations = exprnd(1, n_samples, n_factors);
spectra = concentrations * pure_spectra + 0.10 * randn(n_samples, n_wavelengths);

% --- Maximum likelihood factor analysis with varimax rotation ---
% factoran works on the correlation matrix, so communalities = 1 - psi
[loadings, psi, ~, stats] = factoran(spectra, n_factors, 'rotate', 'varimax');
communalities = 1 - psi;

fprintf('Mean communality: %.3f\n', mean(communalities));
if isfield(stats, 'chisq')
    fprintf('Chi-squared = %.1f, df = %d, p = %.3f\n', stats.chisq, stats.dfe, stats.p);
else
    disp('Heywood case: some uniquenesses hit the lower bound, so no test is available');
end

eigenvalues = sort(eig(corrcoef(spectra)), 'descend');

figure('Position', [100 100 1400 400]);
subplot(1, 3, 1);
plot(1:10, eigenvalues(1:10), 'o-'); yline(1, '--r');
xlabel('Factor number'); ylabel('Eigenvalue'); title('Scree plot'); grid on;

subplot(1, 3, 2);
plot(wavelengths, loadings, 'LineWidth', 1.5);
xlabel('Wavelength (nm)'); ylabel('Loading'); title('Varimax-rotated loadings');
legend('Factor 1', 'Factor 2', 'Factor 3'); grid on;

subplot(1, 3, 3);
plot(wavelengths, communalities, 'g-', 'LineWidth', 2); ylim([0 1.05]);
xlabel('Wavelength (nm)'); ylabel('Communality'); title('Communalities'); grid on;

# --- Simulate a 3-component spectroscopic mixture ---
set.seed(42)
n_samples <- 200; n_wavelengths <- 30; n_factors <- 3
wavelengths <- seq(200, 800, length.out = n_wavelengths)
centers <- c(350, 500, 650)
pure_spectra <- t(sapply(centers, function(c) exp(-0.5 * ((wavelengths - c) / 60)^2)))

concentrations <- matrix(rexp(n_samples * n_factors), n_samples, n_factors)
spectra <- concentrations %*% pure_spectra +
           matrix(rnorm(n_samples * n_wavelengths, sd = 0.10), n_samples)

# --- Maximum likelihood factor analysis with varimax rotation ---
# the optimiser needs more than its default 100 iterations on spectra with little noise
fa_result <- factanal(spectra, factors = n_factors, rotation = "varimax",
                      control = list(opt = list(maxit = 1000)))
loadings_mat <- unclass(fa_result$loadings)
communalities <- 1 - fa_result$uniquenesses

cat(sprintf("Mean communality: %.3f\n", mean(communalities)))
cat(sprintf("Chi-squared = %.1f, df = %d, p = %.3f\n", fa_result$STATISTIC, fa_result$dof, fa_result$PVAL))

eigenvalues <- eigen(cor(spectra))$values

par(mfrow = c(1, 3), mar = c(4.5, 4.5, 3, 1))
plot(1:10, eigenvalues[1:10], type = "b", pch = 16,
     xlab = "Factor number", ylab = "Eigenvalue", main = "Scree plot")
abline(h = 1, col = "red", lty = 2)

matplot(wavelengths, loadings_mat, type = "l", lty = 1, lwd = 2,
        xlab = "Wavelength (nm)", ylab = "Loading", main = "Varimax-rotated loadings")
legend("topright", legend = paste("Factor", 1:n_factors), col = 1:n_factors, lwd = 2)

plot(wavelengths, communalities, type = "l", lwd = 2, col = "forestgreen",
     ylim = c(0, 1.05), xlab = "Wavelength (nm)", ylab = "Communality",
     main = "Communalities")

All three find three eigenvalues of the correlation matrix well above 1 and the rest below it, so Kaiser's rule and the scree plot agree on three factors. In MATLAB and R the chi-squared test does not reject the three-factor model (in MATLAB, chi-squared 365.5 with 348 degrees of freedom, p = 0.25). Even so, a few uniquenesses end up at or near the lower bound of the optimiser: with spectra, mild Heywood cases are the rule rather than the exception.

Factor analysis is a good choice for:

• Counting components in mixtures, the fundamental question of spectroscopic mixture analysis.
• Understanding latent structure, when you believe hidden causes drive the correlations, as in environmental, process or sensory data.
• Resolving mixtures, together with constraints or target testing, as the starting point of curve resolution.
• Grouping correlated variables and understanding what drives them.

Other methods are better for:

• Plain dimensionality reduction. If you only need fewer variables for regression or classification, PCA is simpler and makes fewer assumptions.
• Prediction. Factor scores are not built to predict a response; use PLS or PCR.
• Small data sets. Covariance-based factor analysis needs reasonably many samples per variable; with few samples, the estimates are unstable.
• Nonlinear structure. The factor model is linear.

Factor analysis assumes that a few hidden factors produce the correlations you observe, and it separates the variance they share from the variance unique to each variable. Its loadings are only defined up to a rotation; rotation does not change the fit, it only changes how the result is expressed, and methods like varimax choose the expression that is easiest to read. In chemistry, the factors can be real species, which turns counting factors into counting components, and turns the rotation problem into the central question of curve resolution.

The closest relative is principal component analysis, which has no model for the noise and is often the starting point for factor analysis. Two further methods, not yet covered in this library, are worth knowing about: multivariate curve resolution (MCR-ALS), which uses chemical constraints to recover pure spectra and concentration profiles, and independent component analysis (ICA), which looks for factors that are statistically independent rather than merely uncorrelated, and which relies on the factors being non-Gaussian [22].