2. Basic statistics of discrete distributional data

This vignette shows how to extract descriptive statistics from a table of discrete distributions, at two levels:

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 values

1. 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.66134

Its 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.6236162

compQ() 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    128

Almost all the distributions are left-skewed (most judgements are 6 or 7):

mean(get.MatD.stats(M, "skewness")$mat < 0)   # share of left-skewed cells
#> [1] 0.9538462

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.71

Means 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.117

3. 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.01453

The 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.09477456

The 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.050841

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

B <- WD.vec.mean.simplified(M, support = 1:7)
round(get.MatD.stats(B, "mean")$mat, 2)
#>      informat willingn officeho compete global
#> Mean     5.56     6.08     5.32    5.63   5.72

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.228

About 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.000

The correlation between two single distributions is measured instead by the correlation of their quantile functions, rQQ(), which is the shape component of the distance:

rQQ(get.cell.MatD(M, "F06_M", "global"), get.cell.MatD(M, "F02_M", "global"))
#> [1] 0.7773405

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.6136929