--- title: "Principal Component Analysis for Distributional Data" author: Matthias Templ date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Principal Component Analysis for Distributional Data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "##", fig.width = 7, fig.height = 5, warning = FALSE, message = FALSE ) ``` ```{r setup} 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: - 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: ```{r load-data} data("age") str(age) ``` 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) ```{r explore-data} # 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: ```{r plot-raw, fig.height=6} # 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: ```{r plot-by-sex, fig.height=6} # 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: ```{r create-basis} # 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 ```{r convert-to-bdata} # 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: ```{r plot-bdata, fig.height=6} # 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)) ``` ```{r plot-bdata-sex, fig.height=6} # 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: ```{r run-pca} 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: - **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: ```{r pca-summary} summary(pcaAge) ``` 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: ```{r screeplot} 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 ```{r 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: ```{r score-plot, fig.height=6} # Basic score plot plot(pcaAge, type = "scores") ``` ### Coloring by Sex A key question is whether male and female distributions differ systematically: ```{r scores-by-sex, fig.height=6} # 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: ```{r scores-with-labels, fig.height=6} # 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. ```{r loadings-plot, fig.height=6} # 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`: ```{r loadings-scaled, fig.height=6} # 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: ```{r loadings-clr, fig.height=6} # 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: ```{r covariance, fig.height=6} 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? ```{r reconstruction, fig.height=6} # 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: ```{r boxplot} 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: ```{r biplot, fig.height=6} 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: ```{r grouped-pca} pcaGrouped <- princomp(ageSplines, group = age$data$sex) pcaGrouped ``` ## Comparing Screeplots Do males and females have similar dimensionality of variation? ```{r grouped-screeplot, fig.width=10} plot(pcaGrouped, type = "screeplot") ``` ## Comparing Loadings Are the modes of variation the same for both groups? ```{r grouped-loadings, fig.width=10} 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 ```{r grouped-scores, fig.width=10} 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: ```{r multivariate-pca} pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE) pcaMulti ``` 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 ```{r multi-screeplot} 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 ```{r multi-scores} 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: ```{r multi-eigenfunctions, fig.width=10, fig.height=6} plot(pcaMulti, type = "eigenfunctions", npcs = 2) ``` **How to read this plot**: - Each row is one principal component - Each column is one group (male, female) - Black line: mean distribution for that group - Blue line: positive PC direction (mean + 1 SD along PC) - Red dashed line: negative PC direction (mean - 1 SD along PC) **Interpretation for age data**: - If PC1 shows similar patterns for both males and females (e.g., both shift toward older ages), it captures **joint aging** - communes where both sexes have older populations. - If patterns differ between groups for a PC, that component captures **differential** variation - e.g., communes where male and female age structures diverge. ## 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**: - "Do males and females vary in similar ways?" → Use separate grouped PCA - "Which communes have unusual joint male-female age patterns?" → Use multivariate SFPCA # Predicting Scores for New Data If you have new age distributions, you can project them onto the existing PC space without re-running PCA: ```{r predict-example} # 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: - You have a PCA model from historical data - New observations arrive and you want to position them in the existing PC space - You want to classify or cluster new data based on the learned PC coordinates # Practical Workflow Summary Here's a concise workflow for PCA on distributional data: ```{r workflow-summary, eval=FALSE} # 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 ```{r session-info} sessionInfo() ```