2. Basic statistics of discrete distributional data
This vignette shows how to extract descriptive statistics from a table of discrete distributions, at two levels:
- each cell: the usual statistics of a single distribution (mean, median, quantiles, standard deviation, skewness, …);
- each distributional variable (a column of the table): the Wasserstein barycenter, the Wasserstein variance, covariance and correlation (Irpino and Verde, 2015).
We use the faculty-by-gender groups of the univer data
(see the vignette 1. Building discrete distributional
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)
M
#> a matrix of discrete distributions
#> 26 rows x 5 columns
#> informat willingn officeho
#> F01_M [m=5.35 s=1.73 k=7] [m=6.19 s=1.23 k=7] [m=5.77 s=1.41 k=7]
#> F01_W [m=5.16 s=1.72 k=7] [m=6.42 s=0.674 k=7] [m=5.42 s=1.57 k=7]
#> F02_M [m=6.62 s=0.785 k=7] [m=6.65 s=0.748 k=7] [m=6.3 s=1.19 k=7]
#> F02_W [m=6.47 s=1.07 k=7] [m=6.49 s=0.965 k=7] [m=6.23 s=1.34 k=7]
#> F03_M [m=5.27 s=1.44 k=7] [m=5.97 s=1.11 k=7] [m=5.14 s=1.4 k=7]
#> F03_W [m=5.52 s=1.39 k=7] [m=6.12 s=1.1 k=7] [m=5.4 s=1.28 k=7]
#> F04_M [m=5.95 s=0.74 k=7] [m=6.22 s=0.724 k=7] [m=5.47 s=0.894 k=7]
#> F04_W [m=5.73 s=1.31 k=7] [m=6.03 s=1.43 k=7] [m=5.36 s=1.37 k=7]
#> F05_M [m=6.07 s=0.877 k=7] [m=6.34 s=0.864 k=7] [m=6.17 s=1.05 k=7]
#> F05_W [m=6.06 s=0.942 k=7] [m=6.42 s=0.73 k=7] [m=6.26 s=0.932 k=7]
#> compete global
#> F01_M [m=5.45 s=1.66 k=7] [m=5.87 s=1.5 k=7]
#> F01_W [m=6.26 s=1.02 k=7] [m=6.37 s=0.871 k=7]
#> F02_M [m=6.64 s=0.766 k=7] [m=6.71 s=0.683 k=7]
#> F02_W [m=6.48 s=0.965 k=7] [m=6.41 s=1.07 k=7]
#> F03_M [m=5.35 s=1.39 k=7] [m=5.54 s=1.18 k=7]
#> F03_W [m=5.37 s=1.59 k=7] [m=5.59 s=1.16 k=7]
#> F04_M [m=5.98 s=0.79 k=7] [m=5.98 s=0.724 k=7]
#> F04_W [m=5.73 s=1.39 k=7] [m=5.8 s=1.31 k=7]
#> F05_M [m=6.21 s=0.924 k=7] [m=6.26 s=0.864 k=7]
#> F05_W [m=6.18 s=0.831 k=7] [m=6.27 s=0.809 k=7]
#> ... (only the first rows/columns shown)
#> k = number of support values1. Statistics of a single distribution
Each cell is a distributionD:
d <- get.cell.MatD(M, "F06_M", "global")
d
#> Discrete distribution (7 support values, n = 271)
#> x p cdf
#> 1 0.07011 0.07011
#> 2 0.06273 0.13284
#> 3 0.07749 0.21033
#> 4 0.15129 0.36162
#> 5 0.26199 0.62362
#> 6 0.26937 0.89299
#> 7 0.10701 1.00000
#> mean = 4.70849 std = 1.66134Its moments and quantiles are returned by dedicated functions:
c(mean = meanD(d), std = stdD(d), var = varD(d),
skewness = skewD(d), kurtosis = kurtD(d))
#> mean std var skewness kurtosis
#> 4.708487 1.661336 2.760039 -0.738882 2.761565
compQ(d, c(0.1, 0.25, 0.5, 0.75, 0.9)) # quantiles
#> [1] 2 4 5 6 7
compP(d, c(3, 4, 5)) # cumulative probabilities
#> [1] 0.2103321 0.3616236 0.6236162compQ() uses the left-continuous inverse of the cdf,
\(Q(t) = \min\{x : F(x) \ge t\}\); for
a 1-7 scale it always returns one of the categories.
compP(d, 4) is the share of judgements lower than or equal
to 4 (the “not satisfied” students).
2. Statistics of all the cells
get.MatD.stats() computes a statistic for all the cells
and returns a matrix ($mat). Available statistics are
"mean", "median", "std",
"var", "skewness", "kurtosis",
"min", "max", "quantile" (with
the probability prob), "n" (number of
observations) and "k" (number of support values).
means <- get.MatD.stats(M, "mean")$mat
round(head(means), 2)
#> informat willingn officeho compete global
#> F01_M 5.35 6.19 5.77 5.45 5.87
#> F01_W 5.16 6.42 5.42 6.26 6.37
#> F02_M 6.62 6.65 6.30 6.64 6.71
#> F02_W 6.47 6.49 6.23 6.48 6.41
#> F03_M 5.27 5.97 5.14 5.35 5.54
#> F03_W 5.52 6.12 5.40 5.37 5.59
round(head(get.MatD.stats(M, "std")$mat), 2)
#> informat willingn officeho compete global
#> F01_M 1.73 1.23 1.41 1.66 1.50
#> F01_W 1.72 0.67 1.57 1.02 0.87
#> F02_M 0.79 0.75 1.19 0.77 0.68
#> F02_W 1.07 0.97 1.34 0.96 1.07
#> F03_M 1.44 1.11 1.40 1.39 1.18
#> F03_W 1.39 1.10 1.28 1.59 1.16
round(head(get.MatD.stats(M, "skewness")$mat), 2)
#> informat willingn officeho compete global
#> F01_M -1.56 -1.83 -1.33 -1.41 -1.68
#> F01_W -1.35 -0.74 -1.46 -1.74 -1.26
#> F02_M -2.02 -2.24 -1.56 -2.16 -2.71
#> F02_W -2.20 -2.00 -1.76 -1.82 -2.00
#> F03_M -0.77 -1.68 -0.78 -0.91 -0.86
#> F03_W -1.24 -1.85 -0.86 -1.18 -0.93
head(get.MatD.stats(M, "quantile", prob = 0.1)$mat) # first deciles
#> informat willingn officeho compete global
#> F01_M 2 5 4 3 4
#> F01_W 1 5 2 5 5
#> F02_M 5 6 4 5 6
#> F02_W 5 5 4 5 5
#> F03_M 3 5 3 3 4
#> F03_W 4 5 4 3 4
head(get.MatD.stats(M, "n")$mat) # group sizes
#> informat willingn officeho compete global
#> F01_M 31 31 31 31 31
#> F01_W 19 19 19 19 19
#> F02_M 84 84 84 84 84
#> F02_W 91 91 91 91 91
#> F03_M 127 127 127 127 127
#> F03_W 128 128 128 128 128Almost all the distributions are left-skewed (most judgements are 6 or 7):
These matrices are ordinary R matrices and can be used with any other tool: for instance, the groups with the lowest and highest mean global satisfaction are
round(sort(means[, "global"])[c(1:3, 24:26)], 2)
#> F13_M F13_W F06_M F01_W F02_W F02_M
#> 4.29 4.54 4.71 6.37 6.41 6.71Means alone may hide relevant differences. The share of unsatisfied
students (scores 1-4) can be computed with compP() on all
the cells:
unsat <- apply(M@M, c(1, 2), function(cell) compP(cell[[1]], 4))
round(head(unsat), 3)
#> informat willingn officeho compete global
#> F01_M 0.129 0.065 0.129 0.161 0.161
#> F01_W 0.211 0.000 0.105 0.053 0.053
#> F02_M 0.036 0.036 0.131 0.036 0.036
#> F02_W 0.077 0.066 0.132 0.077 0.066
#> F03_M 0.260 0.087 0.268 0.228 0.157
#> F03_W 0.133 0.047 0.188 0.195 0.1173. The Wasserstein barycenter
The barycenter (Fréchet mean) of a distributional variable is the
distribution whose quantile function is the (weighted) average of the
quantile functions of the cells: it minimizes the sum of squared L2
Wasserstein distances from the data. WD.vec.mean() computes
it exactly:
g <- M[, "global"]
bar <- WD.vec.mean(g)
bar
#> Discrete distribution (107 support values)
#> x p cdf
#> 2.423077 0.00781 0.00781
#> 2.461538 0.00088 0.00870
#> 2.538462 0.00008 0.00877
#> 2.576923 0.00016 0.00893
#> 2.615385 0.00206 0.01099
#> 2.653846 0.00290 0.01389
#> ... (95 rows omitted) ...
#> x p cdf
#> 6.807692 0.00837 0.89299
#> 6.846154 0.00979 0.90278
#> 6.884615 0.02030 0.92308
#> 6.923077 0.00375 0.92683
#> 6.961538 0.02555 0.95238
#> 7.000000 0.04762 1.00000
#> mean = 5.70072 std = 1.01453The barycenter of discrete distributions is discrete, but its support is made of averages of quantiles: it has many values, generally not integers.
A barycenter on the original categories
For presentation it is often useful to express the barycenter on the
categories of the scale. WD.vec.mean.simplified() projects
the exact barycenter on the support of the variable: each value receives
the masses of the points that are nearest to it (for a 1-7 scale, the
value 2 receives the masses in \((1.5,
2.5]\), the value 3 those in \((2.5,
3.5]\), and so on). This nearest-point rule gives the
distribution on the categories closest to the exact barycenter in the
Wasserstein distance. The same projection is available for any
distribution or MatD with simplifyD().
bar_s <- WD.vec.mean.simplified(g, support = 1:7)
bar_s
#> Discrete distribution (7 support values)
#> x p cdf
#> 1 0.00000 0.00000
#> 2 0.00870 0.00870
#> 3 0.04008 0.04878
#> 4 0.10507 0.15385
#> 5 0.15385 0.30769
#> 6 0.45012 0.75781
#> 7 0.24219 1.00000
#> mean = 5.72317 std = 1.11335
WassSqDistD(bar, bar_s) # approximation error
#> [1] 0.09477456The simplified barycenter is an approximation for interpretation and display: all the computations of the package (variances, covariances, k-means prototypes, …) use the exact barycenter.
Barycenter versus pooled distribution
The barycenter is not the distribution obtained by pooling all the respondents. Pooling mixes the groups and adds the differences between them to the variability, while the barycenter describes a “typical group”.
univer$All <- "all"
pooled <- get.cell.MatD(data2MatD(univer, unit = "All", vars = "global", support = 1:7), 1, 1)
n <- get.MatD.stats(g, "n")$mat[, 1]
bar_w <- WD.vec.mean(g, w = n) # weighted by the group sizes
rbind(pooled = c(mean = meanD(pooled), std = stdD(pooled)),
barycenter = c(meanD(bar), stdD(bar)),
weighted_barycenter = c(meanD(bar_w), stdD(bar_w)))
#> mean std
#> pooled 5.71776 1.302630
#> barycenter 5.70072 1.014534
#> weighted_barycenter 5.71776 1.050841Weighting the groups by their sizes, the barycenter has the same mean as the pooled distribution, but a smaller standard deviation.
cmp <- rbind(
data.frame(x = pooled@x, p = pooled@p, type = "pooled respondents"),
data.frame(x = bar_s@x, p = bar_s@p, type = "barycenter (simplified)"),
data.frame(x = WD.vec.mean.simplified(g, w = n, support = 1:7)@x,
p = WD.vec.mean.simplified(g, w = n, support = 1:7)@p,
type = "weighted barycenter (simplified)"))
ggplot(cmp, aes(factor(x), p, fill = type)) +
geom_col(position = "dodge") +
labs(x = "global satisfaction", y = "probability", fill = NULL) +
theme_minimal() + theme(legend.position = "bottom")
With several columns, WD.vec.mean() returns a one-row
MatD with the barycenters of all the variables:
4. Wasserstein variance and its decomposition
The Wasserstein variance of a variable is the mean squared distance of its distributions from the barycenter. Each squared distance can be split into the contributions of the means (position), of the standard deviations (size) and of the shapes: \[d_W^2(x_i, \bar{x}) = (\mu_i - \bar\mu)^2 + (\sigma_i - \bar\sigma)^2 + 2\sigma_i\bar\sigma\,(1 - \rho_{QQ}(x_i, \bar x)).\] Averaging over the groups gives a decomposition of the variance:
WD.var.covar(g)
#> global
#> global 0.5563465
dec <- t(sapply(seq_len(get.MatD.nrows(g)), function(i)
WassSqDistD(get.cell.MatD(g, i, 1), bar, details = TRUE)))
var_dec <- colMeans(dec[, c("POSITION", "SIZE", "SHAPE")])
round(var_dec, 4)
#> POSITION SIZE SHAPE
#> 0.3467 0.0826 0.1271
round(var_dec / sum(var_dec), 3) # shares
#> POSITION SIZE SHAPE
#> 0.623 0.148 0.228About two thirds of the variability of the global satisfaction across
groups is due to differences in the mean judgement; the remaining third
is due to differences in dispersion and shape, which an analysis of the
means alone would ignore. WD.SSQ() returns the total sum of
squares (deviance) of each variable.
5. Covariance, correlation and reliability
Covariances are computed on quantile functions: \(C_{jk} = \sum_i w_i \int_0^1 (Q_{ij}(t) - \bar Q_j(t))(Q_{ik}(t) - \bar Q_k(t))dt\).
round(WD.var.covar(M), 3)
#> informat willingn officeho compete global
#> informat 0.668 0.362 0.366 0.497 0.455
#> willingn 0.362 0.482 0.284 0.361 0.375
#> officeho 0.366 0.284 0.632 0.332 0.342
#> compete 0.497 0.361 0.332 0.614 0.485
#> global 0.455 0.375 0.342 0.485 0.556
R <- WD.correlation(M)
round(R, 3)
#> informat willingn officeho compete global
#> informat 1.000 0.638 0.563 0.776 0.746
#> willingn 0.638 1.000 0.514 0.663 0.725
#> officeho 0.563 0.514 1.000 0.533 0.577
#> compete 0.776 0.663 0.533 1.000 0.830
#> global 0.746 0.725 0.577 0.830 1.000The correlation between two single distributions is measured instead
by the correlation of their quantile functions, rQQ(),
which is the shape component of the distance:
Finally, the internal consistency of a set of items is measured by
the distributional Cronbach’s alpha (GLOB), decomposed into
the part due to the means (AVE) and the part due to the
centred distributions (VAR):
WD.cronbach(M[, c("informat", "willingn", "officeho", "compete")])
#> $alpha
#> GLOB AVE VAR
#> 0.8633627 0.9344928 0.7109092
#>
#> $alpha.if.deleted
#> GLOB AVE VAR
#> informat 0.7955391 0.9066005 0.5693641
#> willingn 0.8328587 0.9015962 0.6854476
#> officeho 0.8704605 0.9543110 0.7109416
#> compete 0.7975606 0.8936676 0.6136929References
- Irpino A., Romano E. (2007). Optimal histogram representation of large data sets: Fisher vs piecewise linear approximation. RNTI E-9, 99-110.
- Irpino A., Verde R. (2015). Basic statistics for distributional symbolic variables: a new metric-based approach. Advances in Data Analysis and Classification, 9(2), 143-175. doi:10.1007/s11634-014-0176-4
- Crocetta C., Irpino A., Antonucci L., Marin C. (2025). A performance indicator and its decomposition according to the impacts of different aspects based on distributional data. Quality & Quantity, 59(3), 2091-2110. doi:10.1007/s11135-024-01956-0