princomp.BData: Principal Component Analysis for Distributional Data
In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

princomp.BData

R Documentation

Principal Component Analysis for Distributional Data

Description

Performs Simplicial Functional Principal Component Analysis (SFPCA) on distributional data where each observation is an entire probability distribution. This function extends the classical 'princomp' to work with BData objects, enabling dimensionality reduction and exploratory analysis of datasets of distributions.

Usage

## S3 method for class 'BData'
princomp(
  x,
  cor = FALSE,
  scores = TRUE,
  covmat = NULL,
  center = NULL,
  group = NULL,
  multivariate = FALSE,
  ...
)

## S3 method for class 'princomp.BData'
print(x, ...)

Arguments

`x`	A princomp.BData object to print.
`cor`	Logical. If 'TRUE', the analysis uses the correlation matrix instead of the covariance matrix. This standardizes each coefficient to unit variance before analysis, which can be useful when coefficients have very different scales. Default is 'FALSE'.
`scores`	Logical. If 'TRUE' (default), the principal component scores are computed for each observation. Scores represent the coordinates of each distribution in the reduced PC space.
`covmat`	An optional pre-computed covariance matrix. If provided, this matrix is used instead of computing the covariance from 'x'. Useful for robust or regularized covariance estimates. Note: When provided with a non-orthonormal basis, this should be the covariance in the transformed (orthonormalized) space.
`center`	An optional centering vector. If 'NULL' (default), the column means of the coefficient matrix are used. Providing a custom center allows for centering around a specific reference distribution.
`group`	An optional factor or vector specifying group membership for each observation. If provided, PCA is performed separately for each group level (unless 'multivariate = TRUE'), and the result is an object of class '"princomp.BData.grouped"' containing a list of PCA results. This is useful for comparing principal components across different subgroups (e.g., males vs females, treatment vs control).
`multivariate`	Logical. Only relevant when 'group' is provided. If 'FALSE' (default), PCA is performed separately for each group. If 'TRUE', performs multivariate SFPCA (2D-SFPCA per Hron et al. 2016) where clr-representations are stacked across groups and a single joint PCA is performed. This captures the covariation between groups and yields common eigenfunctions spanning all groups. Requires that observations are paired across groups (same number of observations per group, same observation index = same subject).
`...`	Additional arguments (currently unused).

Details

## Conceptual Background

In distributional data analysis, each observation is not a single number but an entire probability distribution. Examples include:

Particle size distributions from different samples
Income distributions across different regions
Age distributions of populations over time
Spectral distributions in chemistry or physics

PCA for such data aims to find the main "modes of variation" - the principal ways in which distributions in the dataset differ from each other.

## Method

The function works by:

Converting all input distributions to a common spline basis representation using as.BDataSplines. This yields a matrix where each row is an observation and each column is a spline basis coefficient.
**Handling non-orthonormal basis:** If the basis has a non-identity Gram matrix S, the coefficients are transformed via C' = C \cdot S^{1/2} to orthonormalize the representation before PCA.
Computing the covariance (or correlation) matrix of these (transformed) coefficients.
Performing eigendecomposition to find the principal directions.
Projecting the centered data onto these directions to obtain scores.
**Back-transforming loadings:** If transformation was applied, loadings are back-transformed via L = S^{-1/2} \cdot L' so they remain in the original coefficient space.

## Non-Orthonormal Basis Handling

**Important:** Many spline bases, including ZB-splines (zero-integral B-splines), are **not orthonormal**. For such bases, standard PCA on raw coefficients gives incorrect results because the Euclidean inner product in coefficient space does not correspond to the L2 inner product in function space.

This function automatically detects non-orthonormal bases by checking the Gram matrix S stored in the BayesSpaceBasis. If S \neq I, the function applies the correction described in Hron et al. (2016, Eq. 6):

N^{-1} C' C M b = \rho b

where M is the Gram matrix. This is implemented via Cholesky transformation:

Compute S^{1/2} via Cholesky decomposition: S = L L^T
Transform coefficients: \tilde{C} = C \cdot L
Perform standard PCA on \tilde{C}
Back-transform loadings: V = L^{-1} \cdot \tilde{V}

The scores remain unchanged since they are computed with the orthonormalized representation. The loadings are transformed back to the original coefficient space for interpretability.

**For users:** If you want to use an orthonormal basis, you can either:

Use an explicitly orthonormal spline basis
Provide the correct Gram matrix S when constructing the BayesSpaceBasis (the function will handle the correction)

The Gram matrix S defines the inner product in coefficient space: \langle c_1, c_2 \rangle_S = c_1^T S c_2. If S = I (identity), the basis is orthonormal and no correction is needed.

This approach is analogous to princomp.acomp from the compositions package, which performs PCA on compositional data using log-ratio coordinates.

## Interpretation

Eigenvalues/sdev: The variance explained by each component. Large values indicate important modes of variation.
Eigenfunctions: Directions in coefficient space representing functional perturbations. When added to/subtracted from the mean, they show how distributions change along each PC.
Scores: Coordinates of each observation in the PC space. Similar scores indicate similar distribution shapes.

Use plot.princomp.BData with 'type="eigenfunctions"' to visualize how the distribution shape changes along each principal component direction.

## Important Considerations

The number of meaningful components is limited by both the number of observations and the basis dimension.
If your data shows that PC1 explains nearly 100% of variance, the distributions may vary primarily in one characteristic (e.g., only in scale or only in location). This is not a bug but a property of your data.
For meaningful PCA, ensure your dataset has genuine multivariate variation. Use B2GaussianDistribution with an appropriate covariance structure for simulation studies.

Value

An object of class 'c("princomp.BData", "princomp")' containing:

sdev: Numeric vector of standard deviations of the principal components, i.e., the square roots of the eigenvalues of the covariance matrix. These quantify the amount of variation captured by each component.
loadings: Matrix of eigenfunctions in coefficient form (class '"loadings"'). Each column is an eigenvector representing a principal component direction in the space of spline coefficients. Rows correspond to basis coefficients. **Note:** For non-orthonormal bases, these are back-transformed to the original coefficient space. Kept for compatibility with 'stats::princomp'.
eigenfunctions: Alias for 'loadings'. This is the proper SFPCA terminology (Hron et al. 2016).
center: Numeric vector containing the mean of each spline coefficient, used for centering the data before PCA. This represents the "average distribution" in coefficient space.
scale: Numeric vector of scalings applied to each coefficient. All ones unless 'cor=TRUE', in which case these are the standard deviations used for standardization.
n.obs: Integer giving the number of observations (distributions) in the dataset.
scores: Matrix of principal component scores (if 'scores=TRUE'). Each row corresponds to an observation, each column to a principal component. The score indicates how far each distribution deviates from the mean along that PC direction.
call: The matched function call.
Center: The mean distribution as a BDataSplines object. This can be plotted or converted to other representations (PDF, CDF, etc.).
Loadings: The eigenfunctions as a BDataSplines object. Each row represents one principal component as a functional perturbation. Kept for compatibility.
Eigenfunctions: Alias for 'Loadings'. This is the proper SFPCA terminology.
basis: The BayesSpaceBasis used for the spline representation.
OrigClass: Character string recording the original class of the input data, useful for back-conversion.
orthonormal: Logical indicating whether the basis was orthonormal ('TRUE') or required Gram matrix correction ('FALSE').
GramMatrix: The Gram matrix S from the basis. Identity matrix if the basis is orthonormal.

Methods (by generic)

print(princomp.BData): Print method displaying key PCA results

Author(s)

Matthias Templ

References

Hron, K., Menafoglio, A., Templ, M., Hruzova, K., & Filzmoser, P. (2016). Simplicial principal component analysis for density functions in Bayes spaces. Computational Statistics & Data Analysis, 94, 330-350.

Examples

## Example 1: Basic PCA on simulated distributional data
## -------------------------------------------------------
## We use B2GaussianDistribution to generate distributions with
## genuine multivariate variation (not just scaling differences).

# Define a Gaussian process on Bayes space with exponential covariance
meanDist <- BDataPDFs(list(dunif), basis = B2Unif01Basis7)
covFunction <- function(x, y) exp(-4 * abs(x - y)^1.7)
gaussProc <- B2GaussianDistribution(meanDist, covFunction)

# Simulate 30 distributions from this process
simData <- simulate(gaussProc, nsim = 30, seed = 42)
splineData <- as.BDataSplines(simData)

# Perform PCA
pcaResult <- princomp(splineData)

# View results
pcaResult
summary(pcaResult)

## Example 2: Visualizing and interpreting the PCA results
## ---------------------------------------------------------

# Screeplot: How much variance does each PC explain?
plot(pcaResult, type = "screeplot")
# Interpretation: The first two components (PC1 and PC2) capture most of the
# variance (around 70%), with a clear "elbow" after PC3. This suggests that
# the main patterns in the distributional data can be described by 2-3 modes
# of variation.

# Score plot: How do observations cluster in PC space?
plot(pcaResult, type = "scores")
# Interpretation: Each point represents one distribution. Distributions that
# are close together in this plot have similar shapes. Look for clusters
# (groups of similar distributions), outliers (isolated points), or gradients
# (continuous variation along a direction).

# To identify individual distributions, use labels:
plot(pcaResult, type = "scores", labels = TRUE)
# Now each distribution is labeled with its observation number, making it
# easy to trace outliers or interesting cases back to the original data.

# Loadings plot: How do distributions change along each PC?
# This is THE KEY PLOT for interpretation!
plot(pcaResult, type = "loadings", npcs = 2)
# Interpretation: The black line shows the mean distribution. For each PC,
# the solid colored line shows what a distribution looks like when it has a
# HIGH score on that PC (positive direction), while the dashed line shows a
# LOW score (negative direction). For example, if PC1+ is shifted right and
# PC1- is shifted left compared to the mean, then PC1 captures location
# variation. If PC1+ is wider and PC1- is narrower, then PC1 captures spread.

# Biplot: Combine scores and loadings in one plot
biplot(pcaResult)
# Interpretation: Numbers show observations, red arrows show how spline
# coefficients contribute to the PCs. Observations in similar positions have
# similar distributions. The arrows are harder to interpret directly since
# they represent abstract basis coefficients - use the loadings plot above
# for functional interpretation.

## Example 3: Predicting scores for new data
## -------------------------------------------
# Generate new distributions from the same process
newData <- as.BDataSplines(simulate(gaussProc, nsim = 5, seed = 123))

# Project onto the existing PC space
newScores <- predict(pcaResult, newData)
print(newScores)

## Example 4: PCA on real age distribution data
## ----------------------------------------------
# Load age distribution data from Upper Austrian communes
data("age")

# Create a basis for age distributions (support 0-100 years)
rmAge <- ReferenceMeasureUniform(0, 100)
B2Age <- BayesSpace(rmAge)
basisAge <- BayesSpaceBasis(B2Age, nbasis = 7)

# Convert proportions to BDataDensityHistograms
# The data are proportions in 5-year age bins (last bin is 10 years: 90-100)
breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 100)
ageHist <- BDataDensityHistograms(age$distributions, breaks = breaks, basis = basisAge)

# Convert: histograms -> densities -> splines
ageDensities <- as.BDataDensities(ageHist)
ageSplines <- as.BDataSplines(ageDensities)

# Perform PCA
pcaAge <- princomp(ageSplines)
summary(pcaAge)

# Screeplot: Most variance is in PC1 and PC2
plot(pcaAge, type = "screeplot")

# Score plot colored by sex
sexColors <- ifelse(age$data$sex == "male", "blue", "red")
plot(pcaAge, type = "scores", col = sexColors)
legend("topright", c("Male", "Female"), col = c("blue", "red"), pch = 19)
# Interpretation: Males and females separate along PC1, indicating that the
# primary mode of variation captures sex differences in age structure
# (e.g., females tend to live longer, creating different elderly proportions).

# Loadings plot: What do the PCs represent?
plot(pcaAge, type = "loadings", npcs = 2)
# Interpretation: PC1 captures the balance between younger and older ages.
# PC2 may capture differences in the working-age population bulge.

## Example 5: Grouped PCA - compare males vs females
## ---------------------------------------------------
# Perform PCA separately for each sex group
pcaGrouped <- princomp(ageSplines, group = age$data$sex)
pcaGrouped

# Compare screeplots: Do males and females have similar dimensionality?
plot(pcaGrouped, type = "screeplot")

# Compare loadings: Are the modes of variation the same?
plot(pcaGrouped, type = "loadings", npcs = 2)
# Interpretation: If the loadings look similar across groups, the same
# types of variation exist in both. If they differ, each group has
# distinct patterns of age distribution variation.

# Compare score distributions with boxplots
plot(pcaGrouped, type = "boxplot", npcs = 3)

## Example 6: Multivariate SFPCA (2D-SFPCA) - joint analysis across groups
## -------------------------------------------------------------------------
# Unlike grouped PCA which analyzes groups separately, multivariate SFPCA
# stacks the representations and finds COMMON principal directions across groups.
# This requires paired data: same observation index = same subject across groups.

pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
pcaMulti

# Screeplot: Total variance across both male and female distributions
plot(pcaMulti, type = "screeplot")

# Score plot: Each point is a commune, scored on joint PCs
plot(pcaMulti, type = "scores")

# Eigenfunctions: Shows how each group component of the joint eigenfunction
# perturbs the mean. Left column = male, right column = female.
plot(pcaMulti, type = "eigenfunctions", npcs = 2)
# Interpretation: If PC1 shows similar patterns for both groups, they vary
# together. If patterns differ, the PC captures group-specific variation
# or interactions between groups.

DistributionalDDA documentation built on Oct. 20, 2026, 5:06 p.m.