The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

4. Principal component analysis of discrete distributional data

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

References

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.