Introduction

Principal Component Analysis (PCA) is a fundamental technique for exploring multivariate data. In distributional data analysis, each observation is not a single number but an entire probability distribution. The DistributionalDDA package extends R’s classical princomp() function to work seamlessly with “BData” objects - the package’s representation of distributional data. This means you can use the familiar princomp(), plot(), summary(), predict(), biplot(), and screeplot() functions on your distributional data, with specialized methods that understand the functional nature of distributions.

This vignette demonstrates PCA for distributional data using the age dataset included in the package. This dataset contains age distributions from 57 communes in Upper Austria, separated by sex - a real-world example where we want to understand the main modes of variation in how populations are distributed across age groups.

What PCA for Distributions Reveals

In general, PCA for distributional data identifies the main “modes of variation” - the principal ways in which distributions in a dataset differ from each other. It answers questions like:

How many independent dimensions of variation exist in the data?
What does each dimension represent in terms of distribution shape?
Which distributions are similar or outliers?
How much of the total variation is captured by a reduced representation?

Questions for the Age Data

Applied to the age dataset, we can specifically ask:

What are the main ways in which age distributions differ across communes?
Are there systematic differences between male and female age distributions?
Which communes have similar or unusual age structures?
Can we characterize communes by just a few principal components?

Loading and Exploring the Data

Let’s start by loading the age distribution dataset:

data("age")
str(age)
## List of 3
##  $ data         :'data.frame':   114 obs. of  22 variables:
##   ..$ commune     : chr [1:114] "4010" "4020" "4030" "4040" ...
##   ..$ municipality: chr [1:114] "Linz" "Steyr" "Wels" "Altheim" ...
##   ..$ sex         : Factor w/ 2 levels "male","female": 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ AGE_2       : num [1:114] 0.0517 0.0515 0.0532 0.0479 0.0512 ...
##   ..$ AGE_7       : num [1:114] 0.0446 0.0476 0.0539 0.0464 0.0532 ...
##   ..$ AGE_12      : num [1:114] 0.0452 0.0478 0.0548 0.0543 0.0552 ...
##   ..$ AGE_17      : num [1:114] 0.0507 0.0612 0.0642 0.0651 0.0649 ...
##   ..$ AGE_22      : num [1:114] 0.0689 0.0706 0.0661 0.0654 0.0638 ...
##   ..$ AGE_27      : num [1:114] 0.0824 0.0696 0.0726 0.0624 0.061 ...
##   ..$ AGE_32      : num [1:114] 0.0822 0.0644 0.0713 0.0639 0.0666 ...
##   ..$ AGE_37      : num [1:114] 0.0676 0.0607 0.0616 0.059 0.0704 ...
##   ..$ AGE_42      : num [1:114] 0.0729 0.0693 0.0796 0.0762 0.0799 ...
##   ..$ AGE_47      : num [1:114] 0.0816 0.0794 0.0825 0.0911 0.0934 ...
##   ..$ AGE_52      : num [1:114] 0.0756 0.0811 0.0774 0.0824 0.0771 ...
##   ..$ AGE_57      : num [1:114] 0.0615 0.0662 0.0623 0.0698 0.0716 ...
##   ..$ AGE_62      : num [1:114] 0.0541 0.0577 0.0549 0.0546 0.0552 ...
##   ..$ AGE_67      : num [1:114] 0.0481 0.0529 0.0456 0.0459 0.0404 ...
##   ..$ AGE_72      : num [1:114] 0.0473 0.0506 0.045 0.0451 0.0404 ...
##   ..$ AGE_77      : num [1:114] 0.0265 0.0283 0.0234 0.0319 0.026 ...
##   ..$ AGE_82      : num [1:114] 0.0222 0.0222 0.0181 0.0244 0.0204 ...
##   ..$ AGE_87      : num [1:114] 0.01236 0.01352 0.00991 0.01099 0.00687 ...
##   ..$ AGE_92      : num [1:114] 0.00448 0.00543 0.00362 0.00341 0.00242 ...
##  $ ages         : num [1:19] 2 7 12 17 22 27 32 37 42 47 ...
##  $ distributions: num [1:114, 1:19] 0.0517 0.0515 0.0532 0.0479 0.0512 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:114] "Linz_male" "Steyr_male" "Wels_male" "Altheim_male" ...
##   .. ..$ : chr [1:19] "AGE_2" "AGE_7" "AGE_12" "AGE_17" ...

The dataset contains:

data: A data frame with commune codes, sex, and proportions for 19 age groups
ages: The midpoint ages of each 5-year age group (2, 7, 12, …, 92)
distributions: A matrix where each row is an age distribution (sums to 1)

# Number of observations
nrow(age$data)
## [1] 114

# The age midpoints
age$ages
##  [1]  2  7 12 17 22 27 32 37 42 47 52 57 62 67 72 77 82 87 92

# First few rows of the data
head(age$data[, 1:7])
##   commune           municipality  sex      AGE_2      AGE_7     AGE_12
## 1    4010                   Linz male 0.05167582 0.04464853 0.04522866
## 2    4020                  Steyr male 0.05147977 0.04762422 0.04778713
## 3    4030                   Wels male 0.05318663 0.05388969 0.05480367
## 4    4040                Altheim male 0.04785847 0.04643079 0.05425202
## 5    4041 Gilgenberg am Weilhart male 0.05120178 0.05322157 0.05524137
## 6    4042         Maria Schmolln male 0.04995852 0.05198636 0.05281593
##       AGE_17
## 1 0.05071258
## 2 0.06120011
## 3 0.06422470
## 4 0.06511484
## 5 0.06493638
## 6 0.06728731

# Check that distributions sum to 1
summary(rowSums(age$distributions))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1

Visualizing the Raw Data

Before performing PCA, let’s visualize the age distributions:

# Plot all distributions with transparency, plus the mean
matplot(age$ages, t(age$distributions), type = "l", lty = 1,
        col = adjustcolor("gray", alpha.f = 0.4),
        xlab = "Age", ylab = "Proportion",
        main = "All Age Distributions with Mean")

# Add the mean distribution
meanDist <- colMeans(age$distributions)
lines(age$ages, meanDist, col = "black", lwd = 3)
legend("topright", c("Individual", "Mean"),
       col = c("gray", "black"), lwd = c(1, 3))

We can also visualize the distributions separately by sex:

# Separate by sex
maleIdx <- which(age$data$sex == "male")
femaleIdx <- which(age$data$sex == "female")

# Plot all distributions colored by sex
matplot(age$ages, t(age$distributions[maleIdx, ]), type = "l",
        lty = 1, col = adjustcolor("blue", alpha.f = 0.3),
        ylim = c(0, 0.1), xlab = "Age", ylab = "Proportion",
        main = "Age Distributions by Sex")
matlines(age$ages, t(age$distributions[femaleIdx, ]),
         lty = 1, col = adjustcolor("red", alpha.f = 0.3))

# Add means for each sex
lines(age$ages, colMeans(age$distributions[maleIdx, ]), col = "blue", lwd = 3)
lines(age$ages, colMeans(age$distributions[femaleIdx, ]), col = "red", lwd = 3)
legend("topright", c("Male", "Female"), col = c("blue", "red"), lwd = 3)

Converting to BData Format

To perform PCA, we need to convert our data into the “BData” format used by the package. This involves:

Creating a suitable Bayes space with a reference measure
Defining a basis for finite-dimensional approximation
Converting the histogram data to “BDataSplines”

Step 1: Create the Bayes Space

Age distributions have support on [0, 100] (approximately), so we use a uniform reference measure on this interval:

# Reference measure: Uniform on [0, 100]
rmAge <- ReferenceMeasureUniform(0, 100)

# Create the Bayes space
B2Age <- BayesSpace(rmAge)
B2Age
## BayesSpace with reference measure:ReferenceMeasureUniform(a = 0, b = 100)

# Create a basis with 7 dimensions
# (this provides a good balance between flexibility and parsimony)
basisAge <- BayesSpaceBasis(B2Age, nbasis = 7)
basisAge
## BayesSpaceBasis.ReferenceMeasureUniform(basis = basis$referenceMeasure, 
##     bayesSpace = basis, nbasis = 7)

Step 2: Convert Histogram Data to BData

The age distributions are given as proportions in 5-year age bins. We treat these as density histograms and convert them through the pipeline: histograms -> densities -> splines

# Define the histogram breaks
# Note: The last bin (90-100) is 10 years wide, others are 5 years
breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50,
            55, 60, 65, 70, 75, 80, 85, 90, 100)

# Create density histograms from the proportions
ageHist <- BDataDensityHistograms(age$distributions,
                                   breaks = breaks,
                                   basis = basisAge)
class(ageHist)
## [1] "BDataDensityHistograms" "BDataHistograms"        "BData"
dim(ageHist)
## [1] 114  19

# Convert to continuous densities
ageDensities <- as.BDataDensities(ageHist)

# Convert to spline representation (needed for PCA)
ageSplines <- as.BDataSplines(ageDensities)
class(ageSplines)
## [1] "BDataSplines" "BData"
dim(ageSplines)
## [1] 114   7

Step 3: Visualize the BData

Let’s check that our conversion worked correctly:

# Plot all distributions as PDFs with transparency
plot(ageSplines, as = "PDF", col = adjustcolor("gray", alpha.f = 0.4),
     main = "All Age Distributions (as PDFs) with Mean")

# Add the mean distribution
meanAge <- mean(ageSplines)
plot(meanAge, as = "PDF", add = TRUE, col = "black", lwd = 3)
legend("topright", c("Individual", "Mean"), col = c("gray", "black"), lwd = c(1, 3))

# Plot by sex with transparency
maleSplines <- BDataSplines(unclass(ageSplines)[maleIdx, ], basis = basisAge)
femaleSplines <- BDataSplines(unclass(ageSplines)[femaleIdx, ], basis = basisAge)

plot(maleSplines, as = "PDF", col = adjustcolor("blue", alpha.f = 0.3),
     main = "Age Distributions by Sex (as PDFs)")
plot(femaleSplines, as = "PDF", col = adjustcolor("red", alpha.f = 0.3), add = TRUE)

# Add sex-specific means
plot(mean(maleSplines), as = "PDF", add = TRUE, col = "blue", lwd = 3)
plot(mean(femaleSplines), as = "PDF", add = TRUE, col = "red", lwd = 3)
legend("topright", c("Male", "Female"), col = c("blue", "red"), lwd = 3)

Performing PCA

Now we can perform PCA on our distributional data. The syntax is identical to R’s classical princomp() function - simply pass the “BData” object:

pcaAge <- princomp(ageSplines)
pcaAge
## Call:
## princomp(x = ageSplines)
## 
## Standard deviations:
##     Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6     Comp.7 
## 0.23002255 0.03460215 0.02770435 0.02165770 0.01681350 0.01309618 0.01272954 
## 
## PC1 explains 94.5 % of the variance, PC1-2 explain 96.7 %
## Basis dimension: 7 
## Number of observations: 114 
## Basis: non-orthonormal (Gram matrix correction applied)

Behind the scenes, R dispatches to the princomp.BData() method, which:

Extracts the spline coefficients from the “BDataSplines” object
Centers the coefficient matrix
Computes the covariance matrix and its eigendecomposition
Returns a “princomp.BData” object that extends the standard “princomp” class

The output shows familiar PCA information (standard deviations, variance explained) plus additional details specific to distributional data:

Basis dimension: The number of spline coefficients (here 7), which determines the maximum number of principal components
Number of observations: The number of distributions in the dataset

Summary Statistics

The summary() method works just like for standard PCA:

summary(pcaAge)
## Importance of components:
##                        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## Standard deviation      0.230  0.035  0.028  0.022  0.017  0.013  0.013
## Proportion of Variance  0.945  0.021  0.014  0.008  0.005  0.003  0.003
## Cumulative Proportion   0.945  0.967  0.981  0.989  0.994  0.997  1.000

The summary shows:

Standard deviation: The square root of each eigenvalue - how much variation each PC captures
Proportion of Variance: The fraction of total variance explained by each PC
Cumulative Proportion: Running total of variance explained

Visualizing PCA Results

The plot() method for “princomp.BData” objects provides several plot types for exploring and interpreting PCA results. The following table summarizes all available options:

Type	Description
`"screeplot"`	Bar chart of variance explained by each PC; helps decide how many components to retain
`"variance"`	Cumulative variance plot with reference lines at 80%, 90%, 95%
`"scores"`	Scatter plot of observations in PC space; reveals clusters and outliers
`"loadings"`	Shows how distribution shape changes along each PC direction
`"biplot"`	Combines scores and loadings in one plot
`"covariance"`	Heatmap of functional covariance surface across the domain
`"reconstruction"`	Distributions reconstructed from a reduced number of PCs
`"boxplot"`	Boxplots of score distributions for each PC

Screeplot: How Many Components?

The screeplot helps decide how many principal components are meaningful:

plot(pcaAge, type = "screeplot")

Interpretation: Look for an “elbow” where the variance drops off. Components before the elbow capture the main structure in the data. Here, PC1 captures most of the variance, with PC2 adding some additional information.

Cumulative Variance Plot

plot(pcaAge, type = "variance")

Interpretation: The dashed lines show common thresholds (80%, 90%, 95%). This helps decide how many components to retain based on a target variance explained.

Score Plot: Exploring Structure

The score plot shows each distribution as a point in PC space:

# Basic score plot
plot(pcaAge, type = "scores")

Coloring by Sex

A key question is whether male and female distributions differ systematically:

# Create color vector
sexColors <- ifelse(age$data$sex == "male", "blue", "red")

# Score plot colored by sex
plot(pcaAge, type = "scores", col = sexColors)
legend("topright", c("Male", "Female"), col = c("blue", "red"), pch = 19)

Interpretation: If males and females separate along PC1, this suggests that the primary mode of variation in age distributions captures sex differences. This makes sense because women tend to live longer, creating different proportions in older age groups.

Using Labels

To identify specific distributions, we can display labels instead of points. With labels = TRUE, the row names are used - in this case, the commune codes (e.g., “4101”, “4102”, etc.) combined with sex:

# Use row names as labels (commune codes)
plot(pcaAge, type = "scores", labels = TRUE, col = sexColors, cex = 0.7)

Interpretation: The labels allow us to identify which communes have unusual age structures. Communes appearing far from the center or separated from the main cluster may have distinctive demographic characteristics - for example, urban centers with younger populations or rural areas with aging communities. Pairs of labels (same commune code, different colors) that are far apart indicate communes where male and female age distributions differ substantially.

Loadings Plot: The Key to Interpretation

This is the most important plot for understanding what PCA reveals.

The loadings plot shows how the distribution shape changes along each principal component direction. For each PC, two perturbations of the mean distribution are shown:

PC+: mean + scale.sdev × sdev × loading (positive direction)
PC-: mean − scale.sdev × sdev × loading (negative direction)

where sdev is the standard deviation (square root of the eigenvalue) for that PC, and loading is the unit-length eigenvector. By default, scale.sdev = 1, so the perturbations represent distributions that are ±1 standard deviation away from the mean along each PC direction.

# Use colorblind-friendly colors: black (mean), blue (PC1), orange (PC2)
cbPalette <- c("black", "blue", "orange")
plot(pcaAge, type = "loadings", npcs = 2, col = cbPalette)

How to read this plot:

Black line: The mean (average) age distribution
Blue lines: PC1 - solid for PC1+, dashed for PC1-
Orange lines: PC2 - solid for PC2+, dashed for PC2-

Detailed interpretation for the age data:

PC1 (blue lines) captures the dominant mode of variation:

PC1+ (solid blue): Higher proportions for young children (ages 0-12) and elderly (ages 60+), but lower proportions for working-age adults (roughly ages 15-45). This pattern is typical for female populations: women live longer (more elderly) and the male-female ratio shifts with age.
PC1- (dashed blue): Close to the mean for children, but higher proportions in the working-age range (15-45) and lower proportions among the elderly. This pattern is typical for male populations.

Since males and females separate along PC1 (as seen in the score plot), PC1 essentially captures the sex difference in age structure - primarily driven by women’s longer life expectancy creating higher proportions in older age groups.

PC2 (orange lines) captures the second mode of variation:

PC2+ (solid orange): Higher proportions in middle-aged groups (roughly 35-55) and lower proportions for young adults (20-35) and the very old (80+).
PC2- (dashed orange): The opposite pattern - more young adults and very old, fewer middle-aged.

PC2 likely captures demographic differences between communes - for example, urban areas attracting young adults (negative PC2) versus suburban/rural areas with established families and middle-aged residents (positive PC2).

Adjusting Scale

The scale.sdev parameter controls how far from the mean the perturbed distributions are shown. By default (scale.sdev = 1), the plot displays distributions at ±1 standard deviation along each PC direction - this represents a “typical” deviation from the mean.

To see more extreme perturbations, increase scale.sdev:

# Show +/- 2 standard deviations from mean
plot(pcaAge, type = "loadings", npcs = 2, scale.sdev = 2, col = cbPalette)

When to adjust the scale:

Use scale.sdev = 2 or higher when the differences between solid and dashed lines are subtle and hard to see at the default scale
Use scale.sdev = 0.5 if the perturbations look unrealistically extreme
The appropriate scale depends on your data - choose a value that shows meaningful differences without producing implausible distributions

Alternative Representations

The loadings can be displayed in different forms using the as parameter. Each representation emphasizes different aspects of the distributions:

`as`	Description
`"PDF"`	Probability density function (default) - shows proper densities
`"logPDF"`	Log of PDF - emphasizes differences in the tails
`"clr"`	Centered log-ratio - the natural Bayes space representation
`"density"`	Density with respect to the reference measure
`"logDensity"`	Log density with respect to reference measure

The CLR (centered log-ratio) representation is particularly useful because it is the representation in the subspace \(L^2_0\) of zero-integral square integrable log-density functions, a subspace isometric to the Bayes space - and hence the set where PCA is actually performed:

# Centered log-ratio (CLR) representation
plot(pcaAge, type = "loadings", npcs = 2, as = "clr", col = cbPalette)

Interpretation: In the CLR representation, values are centered around zero. Positive values indicate age groups with higher-than-geometric-mean proportions; negative values indicate lower-than-geometric-mean proportions. This makes it easy to see which parts of the age range increase or decrease along each PC direction.

Covariance Surface

The covariance surface shows which parts of the age range co-vary:

plot(pcaAge, type = "covariance")

Interpretation:

High values (red) on the diagonal show high variance at that age
Off-diagonal patterns show how different age groups co-vary
Negative values indicate that when one age group is above average, another tends to be below average

Reconstruction: Checking Approximation Quality

How well do the first few PCs approximate the original data?

# Reconstruct using only 2 PCs with transparency
sexColorsTransp <- ifelse(age$data$sex == "male",
                          adjustcolor("blue", alpha.f = 0.3),
                          adjustcolor("red", alpha.f = 0.3))
plot(pcaAge, type = "reconstruction", npcs = 2, col = sexColorsTransp)

Compare this to the original distributions. If the reconstructed curves look similar to the originals, 2 PCs capture most of the structure. With only 2 PCs, some detail is lost, but the main patterns (especially the male/female separation) are preserved.

Boxplot of Scores

To see how the PC scores vary across all age distributions:

plot(pcaAge, type = "boxplot", npcs = 4)

This helps identify potential outliers (points beyond the whiskers) and compare the spread along each PC direction.

Biplot

The biplot combines scores and loadings in one visualization:

biplot(pcaAge)

Interpretation:

Numbers represent observations (distributions)
Red arrows represent spline coefficients
Observations close together have similar distributions
The arrows are harder to interpret directly (they are abstract basis functions )

For functional interpretation, the loadings plot (type = "loadings") is more informative.

Grouped PCA: Comparing Subgroups

When comparing groups (e.g., males vs females), we can perform PCA separately for each group:

pcaGrouped <- princomp(ageSplines, group = age$data$sex)
pcaGrouped
## Grouped PCA on distributional data
## Groups: male, female 
## 
## --- Group: male ---
##   Observations: 57 
##   PC1 explains: 79.8 % of variance
## --- Group: female ---
##   Observations: 57 
##   PC1 explains: 80.2 % of variance

Comparing Screeplots

Do males and females have similar dimensionality of variation?

plot(pcaGrouped, type = "screeplot")

Comparing Loadings

Are the modes of variation the same for both groups?

plot(pcaGrouped, type = "loadings", npcs = 2)

Interpretation: If the loadings look similar across groups, both sexes vary in similar ways. If they differ, each group has distinct patterns of age distribution variation.

Comparing Score Distributions

plot(pcaGrouped, type = "scores")

Multivariate SFPCA: Joint Analysis Across Groups

The grouped PCA above performs separate PCA for each group. But what if we want to analyze how male and female age distributions jointly vary across communes? For example, communes where males have older age structures might also have females with older age structures - this covariation between groups is not captured by separate PCAs.

Multivariate SFPCA (also called 2D-SFPCA, following Hron et al. 2016) solves this by:

Stacking the clr-representations from all groups into a single vector per observation
Performing a single PCA on this combined representation
Finding common principal directions that span all groups

This requires paired observations: the same index must correspond to the same subject (commune) across groups. In our data, each commune has both a male and female age distribution - these are paired.

Running Multivariate SFPCA

Simply add multivariate = TRUE to the princomp() call:

pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
pcaMulti
## Multivariate SFPCA (2D-SFPCA) on distributional data
## Groups: male, female 
## Observations per group: 57 
## Coefficients per group: 7 
## Total dimensions: 14 
## Note: Non-orthonormal basis, Gram matrix correction applied
## 
## Importance of components:
##                PC1       PC2       PC3       PC4       PC5        
##   Std. Dev.    0.1352    0.0642    0.0419    0.0303    0.0256     
##   Prop. Var.   0.6703    0.1510    0.0644    0.0337    0.0240     
##   Cum. Prop.   0.6703    0.8213    0.8857    0.9194    0.9434

The output shows:

Groups: The factor levels being analyzed jointly (male, female)
Observations per group: Number of paired observations (57 communes)
Total dimensions: Coefficients per group × number of groups (7 × 2 = 14)

Screeplot: Variance Across All Groups

plot(pcaMulti, type = "screeplot")

Interpretation: This shows the total variance explained across both male and female distributions. Compare this to the separate grouped PCA - the joint analysis typically explains variance differently because it captures cross-group covariation.

Score Plot: Communes in Joint PC Space

plot(pcaMulti, type = "scores")

Interpretation: Each point represents a commune (not a single distribution). Communes close together have similar age structures for both males and females. Outliers are communes with unusual joint age patterns.

Eigenfunctions: How Groups Vary Together

This is the key plot for interpretation. It shows how distributions change along each PC, with separate panels for each group:

plot(pcaMulti, type = "eigenfunctions", npcs = 2)