DistributionalDDA
Distributional Data Analysis in Bayes Spaces

Package index
Search the DistributionalDDA package
Vignettes
Functions
1086
Source code
56
Man pages
162
Home
/
CRAN
/
DistributionalDDA
/
Principal Component Analysis for Distributional Data

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

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "##",
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE
)

library(DistributionalDDA)

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:

Questions for the Age Data

Applied to the age dataset, we can specifically ask:

Loading and Exploring the Data

Let's start by loading the age distribution dataset:

data("age")
str(age)

The dataset contains:

# Number of observations
nrow(age$data)

# The age midpoints
age$ages

# First few rows of the data
head(age$data[, 1:7])

# Check that distributions sum to 1
summary(rowSums(age$distributions))

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:

  1. Creating a suitable Bayes space with a reference measure
  2. Defining a basis for finite-dimensional approximation
  3. 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

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

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)
dim(ageHist)

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

# Convert to spline representation (needed for PCA)
ageSplines <- as.BDataSplines(ageDensities)
class(ageSplines)
dim(ageSplines)

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

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

  1. Extracts the spline coefficients from the "BDataSplines" object

  2. Centers the coefficient matrix

  3. Computes the covariance matrix and its eigendecomposition

  4. 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:

Summary Statistics

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

summary(pcaAge)

The summary shows:

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:

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:

Detailed interpretation for the age data:

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

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 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:

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:

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:

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

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:

  1. Stacking the clr-representations from all groups into a single vector per observation
  2. Performing a single PCA on this combined representation
  3. 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

The output shows:

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)

How to read this plot:

Interpretation for age data:

When to Use Multivariate vs Separate Grouped PCA

| Approach | Use when... | |:---------|:------------| | Separate (multivariate = FALSE) | You want to compare modes of variation between groups | | Joint (multivariate = TRUE) | You have paired data and want to analyze cross-group covariation |

Example questions:

Predicting Scores for New Data

If you have new age distributions, you can project them onto the existing PC space without re-running PCA:

# Suppose we have the first 5 observations as "training" data
# and want to project observations 6-10 as "new" data

# Create a subset as new BDataSplines (subsetting keeps the class)
newAgeSplines <- BDataSplines(unclass(ageSplines)[6:10, ], basis = basisAge)

# The pcaAge model was fit on all data, but let's demonstrate predict():
newScores <- predict(pcaAge, newAgeSplines)
print(newScores[, 1:3])  # Show first 3 PCs

# Compare with the actual scores
actualScores <- pcaAge$scores[6:10, 1:3]
print(actualScores)

# They should be identical (since we used the same data)

In practice, predict() is useful when:

Practical Workflow Summary

Here's a concise workflow for PCA on distributional data:

# 1. Load data
data("age")

# 2. Create Bayes space and basis
rmAge <- ReferenceMeasureUniform(0, 100)
basisAge <- BayesSpaceBasis(BayesSpace(rmAge), nbasis = 7)

# 3. Convert to BData (adapt based on your data format)
breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50,
            55, 60, 65, 70, 75, 80, 85, 90, 100)
ageSplines <- 
    BDataDensityHistograms(age$distributions, 
                           breaks = breaks, 
                           basis = basisAge) |>
  as.BataDensities() |> 
  as.BDataSplines() # pipe |> instead of nesting 

# 4. Run PCA
pcaResult <- princomp(ageSplines)

# 5. Explore results
summary(pcaResult)                              # Variance explained
plot(pcaResult, type = "screeplot")             # How many components?
plot(pcaResult, type = "loadings", npcs = 2)    # What do components mean?
plot(pcaResult, type = "scores", col = myColors) # Where are my observations?

# 6. Optional: Separate grouped analysis (compare groups)
pcaGrouped <- princomp(ageSplines, group = age$data$sex)
plot(pcaGrouped, type = "loadings", npcs = 2)   # Compare groups

# 7. Optional: Multivariate SFPCA (joint analysis, requires paired data)
pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
plot(pcaMulti, type = "eigenfunctions", npcs = 2)  # Cross-group covariation

Tips and Best Practices

  1. Choose the right basis dimension: Too few dimensions may lose important variation; too many may introduce noise. Start with 5-7 and adjust based on how well reconstructions match originals.

  2. Use the loadings plot: The type = "loadings" plot is the key to interpreting what each PC means in terms of distribution shape.

  3. Consider the reference measure: The reference measure should match your data's support. Use ReferenceMeasureUniform(a, b) for bounded data on [a, b].

  4. Check variance explained: If PC1 explains nearly 100% of variance, your data may vary primarily in one characteristic. This is not a bug but a property of your data.

  5. Use grouped PCA for comparison: When comparing subgroups, grouped PCA reveals whether they have the same or different modes of variation.

Session Info

sessionInfo()


Try the DistributionalDDA package in your browser

Any scripts or data that you put into this service are public.

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

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close