## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
has_cub <- requireNamespace("CUB", quietly = TRUE)
knitr::opts_chunk$set(eval = has_cub)

## ----eval = !has_cub, echo = FALSE, results = "asis"--------------------------
# cat("**This vignette needs the 'CUB' package (for the `univer` data), which is not installed: code is shown but not run.**")

## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
d <- get.cell.MatD(M, "F06_M", "global")
d

## -----------------------------------------------------------------------------
c(mean = meanD(d), std = stdD(d), var = varD(d),
  skewness = skewD(d), kurtosis = kurtD(d))
compQ(d, c(0.1, 0.25, 0.5, 0.75, 0.9))   # quantiles
compP(d, c(3, 4, 5))                     # cumulative probabilities

## -----------------------------------------------------------------------------
means <- get.MatD.stats(M, "mean")$mat
round(head(means), 2)
round(head(get.MatD.stats(M, "std")$mat), 2)
round(head(get.MatD.stats(M, "skewness")$mat), 2)
head(get.MatD.stats(M, "quantile", prob = 0.1)$mat)   # first deciles
head(get.MatD.stats(M, "n")$mat)                      # group sizes

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

## -----------------------------------------------------------------------------
round(sort(means[, "global"])[c(1:3, 24:26)], 2)

## -----------------------------------------------------------------------------
unsat <- apply(M@M, c(1, 2), function(cell) compP(cell[[1]], 4))
round(head(unsat), 3)

## -----------------------------------------------------------------------------
g <- M[, "global"]
bar <- WD.vec.mean(g)
bar

## -----------------------------------------------------------------------------
bar_s <- WD.vec.mean.simplified(g, support = 1:7)
bar_s
WassSqDistD(bar, bar_s)       # approximation error

## -----------------------------------------------------------------------------
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)))

## ----fig.height = 3.5---------------------------------------------------------
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")

## -----------------------------------------------------------------------------
B <- WD.vec.mean.simplified(M, support = 1:7)
round(get.MatD.stats(B, "mean")$mat, 2)

## -----------------------------------------------------------------------------
WD.var.covar(g)
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)
round(var_dec / sum(var_dec), 3)       # shares

## -----------------------------------------------------------------------------
round(WD.var.covar(M), 3)
R <- WD.correlation(M)
round(R, 3)

## -----------------------------------------------------------------------------
rQQ(get.cell.MatD(M, "F06_M", "global"), get.cell.MatD(M, "F02_M", "global"))

## -----------------------------------------------------------------------------
WD.cronbach(M[, c("informat", "willingn", "officeho", "compete")])