In the L2 Wasserstein space a distribution is represented by its
quantile function. Principal component analysis of distributional data
can therefore be performed on the matrix of quantiles:
each distribution is described by its quantiles at a grid of probability
levels, and a non-standardized PCA of this matrix approximates the
Wasserstein PCA (Verde, Irpino and Balzanella, 2016). With several
distributional variables, each one is a group of quantile columns and
the groups are balanced by a multiple factor analysis
(MFA), as proposed by Crocetta et al. (2025) to build a performance
indicator from customer satisfaction data.
library(DDAWASS)
library(ggplot2)
data(univer, package = "CUB")
univer$Faculty <- factor(sprintf("F%02d", univer$faculty))
univer$Gender <- factor(univer$gender, levels = 0:1, labels = c("M", "W"))
items <- c("informat", "willingn", "officeho", "compete", "global")
M <- data2MatD(univer, unit = c("Faculty", "Gender"), vars = items,
support = 1:7, min.n = 10)
gender <- sub(".*_", "", get.MatD.rownames(M))
The matrix of quantiles
WD.quantile.matrix() computes the quantiles of each
distribution at the levels \(0, 1/q, \dots,
1\). For a 1-7 scale the minimum and maximum are often 1 and 7
for all groups, and these columns carry no information; the PCA
functions add a negligible noise to constant columns.
Q <- WD.quantile.matrix(M, quantiles = 10, vars = "global")
Q[1:5, ]
#> global.Min global.Q(0.10) global.Q(0.20) global.Q(0.30) global.Q(0.40)
#> F01_M 1 4 5 6 6
#> F01_W 1 5 6 6 6
#> F02_M 1 6 7 7 7
#> F02_W 1 5 6 6 7
#> F03_M 1 4 5 5 5
#> global.Q(0.50) global.Q(0.60) global.Q(0.70) global.Q(0.80)
#> F01_M 6 7 7 7
#> F01_W 7 7 7 7
#> F02_M 7 7 7 7
#> F02_W 7 7 7 7
#> F03_M 6 6 6 7
#> global.Q(0.90) global.Max
#> F01_M 7 7
#> F01_W 7 7
#> F02_M 7 7
#> F02_W 7 7
#> F03_M 7 7
PCA of one distributional variable
WD.1d.PCA() analyses a single variable; the means and
the standard deviations of the distributions are projected as
supplementary variables to help the interpretation of the axes.
pca_g <- WD.1d.PCA(M, var = "global", quantiles = 10, plots = FALSE)
round(pca_g$eig[1:4, ], 2)
#> eigenvalue percentage of variance cumulative percentage of variance
#> comp 1 0.35 81.68 81.68
#> comp 2 0.02 4.94 86.62
#> comp 3 0.02 3.56 90.18
#> comp 4 0.01 2.66 92.84
round(pca_g$quanti.sup$cor[, 1:2], 3)
#> Dim.1 Dim.2
#> Mean 0.990 0.010
#> Std -0.721 0.301
The first component summarizes most of the variability; it is almost
perfectly correlated with the mean of the distributions and negatively
correlated with their standard deviation: it opposes satisfied groups
with concentrated judgements to less satisfied groups with dispersed
judgements. The second component explains a small share of variability
and describes minor differences in the shape of the distributions.
The Spanish-fan plot shows the quantile variables as
consecutive triangles, from the minimum to the maximum (red: extreme
quantiles, yellow: median):
WD.plot.Spanish.fan(pca_g)
Groups can be plotted on the first plane, coloured by their mean
global satisfaction:
ind <- data.frame(pca_g$ind$coord[, 1:2], label = rownames(pca_g$ind$coord),
mean = get.MatD.stats(M, "mean")$mat[, "global"], gender = gender)
ggplot(ind, aes(Dim.1, Dim.2, colour = mean, shape = gender, label = label)) +
geom_hline(yintercept = 0, colour = "grey70") + geom_vline(xintercept = 0, colour = "grey70") +
geom_point(size = 2.5) + geom_text(vjust = -0.8, size = 3, show.legend = FALSE) +
scale_colour_gradient(low = "#c62320", high = "#275024") +
labs(title = "Global satisfaction: groups on the first factorial plane",
x = sprintf("Dim 1 (%.1f%%)", pca_g$eig[1, 2]),
y = sprintf("Dim 2 (%.1f%%)", pca_g$eig[2, 2])) +
theme_minimal()
Multiple factor analysis of several items
WD.MultiplePCA() analyses several distributional
variables together. Here we use the four specific items, keeping
global aside (it is the response of the regression
vignette).
mfa <- WD.MultiplePCA(M, list.of.vars = c("informat", "willingn", "officeho", "compete"),
quantiles = 6, plots = FALSE)
round(mfa$eig[1:4, ], 2)
#> eigenvalue percentage of variance cumulative percentage of variance
#> comp 1 3.36 69.41 69.41
#> comp 2 0.46 9.39 78.80
#> comp 3 0.21 4.31 83.12
#> comp 4 0.18 3.69 86.81
The first factor is common to all the items; the contributions and
the coordinates of the groups of quantiles (one group per item) show
which items drive each factor:
round(mfa$group$contrib[, 1:2], 1)
#> Dim.1 Dim.2
#> informat 26.4 9.5
#> willingn 27.4 11.3
#> officeho 20.1 71.5
#> compete 26.1 7.7
round(mfa$group$coord[, 1:2], 3)
#> Dim.1 Dim.2
#> informat 0.890 0.043
#> willingn 0.922 0.052
#> officeho 0.675 0.325
#> compete 0.878 0.035
The first factor is a general satisfaction factor shared by the four
items, while the second one is mostly specific to officeho:
the evaluation of the office hours follows a different pattern from the
evaluation of the staff and of the information.
ind <- data.frame(mfa$ind$coord[, 1:2], label = rownames(mfa$ind$coord), gender = gender)
ggplot(ind, aes(Dim.1, Dim.2, colour = gender, label = label)) +
geom_hline(yintercept = 0, colour = "grey70") + geom_vline(xintercept = 0, colour = "grey70") +
geom_point(size = 2.5) + geom_text(vjust = -0.8, size = 3, show.legend = FALSE) +
labs(title = "MFA of the four items", x = sprintf("Dim 1 (%.1f%%)", mfa$eig[1, 2]),
y = sprintf("Dim 2 (%.1f%%)", mfa$eig[2, 2])) +
theme_minimal()
The coordinates of the groups on the first factor can be used as a
synthetic performance indicator of the orientation
services, which takes into account the whole distributions of the
judgements and not only their means (Crocetta et al., 2025):
perf <- sort(mfa$ind$coord[, 1], decreasing = TRUE)
round(head(perf, 5), 2)
#> F02_M F02_W F08_W F05_W F05_M
#> 3.19 3.19 1.83 1.82 1.69
round(tail(perf, 5), 2)
#> F09_W F06_W F13_W F06_M F13_M
#> -1.65 -2.65 -3.38 -3.56 -3.91