## ----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)
# covariates of the groups, recovered from the labels
groups <- data.frame(label = get.MatD.rownames(M))
groups$Faculty <- sub("_.*", "", groups$label)
groups$Gender  <- sub(".*_", "", groups$label)

## -----------------------------------------------------------------------------
D <- WD_MAT_DIST(M)
round(D[1:6, 1:6], 3)

## ----fig.height = 5-----------------------------------------------------------
hc <- WD_hclust(M, method = "ward.D2")
plot(hc, main = "Ward clustering of faculty x gender groups", xlab = "", sub = "")
rect.hclust(hc, k = 3)

## -----------------------------------------------------------------------------
cl_h <- cutree(hc, k = 3)
table(cluster = cl_h, gender = groups$Gender)
same <- tapply(cl_h, groups$Faculty, function(z) length(unique(z)) == 1)
mean(same)   # share of faculties whose two groups are in the same cluster

## -----------------------------------------------------------------------------
set.seed(123)
qual <- sapply(2:7, function(k) WD_kmeans(M, k, rep = 10)$quality)
ggplot(data.frame(k = 2:7, quality = qual), aes(k, quality)) +
  geom_line() + geom_point() +
  labs(y = "between SS / total SS", title = "Quality of the partition") +
  theme_minimal()

## -----------------------------------------------------------------------------
set.seed(123)
km <- WD_kmeans(M, k = 3, rep = 20)
km$quality
km$solution$cardinality
split(names(km$solution$IDX), km$solution$IDX)

## -----------------------------------------------------------------------------
table(kmeans = km$solution$IDX, hclust = cl_h)

## ----fig.height = 4-----------------------------------------------------------
proto <- km$solution$centers
ord <- order(get.MatD.stats(proto, "mean")$mat[, "global"])
proto <- proto[ord, ]
rownames(proto@M) <- c("Low", "Medium", "High")
round(get.MatD.stats(proto, "mean")$mat, 2)
round(get.MatD.stats(proto, "quantile", prob = 0.1)$mat, 2)
plot(proto)

## -----------------------------------------------------------------------------
comp <- t(sapply(items, function(v)
  WassSqDistD(get.cell.MatD(proto, 1, v), get.cell.MatD(proto, 3, v), details = TRUE)))
round(comp, 3)

## -----------------------------------------------------------------------------
set.seed(123)
km_s <- WD_kmeans(M, k = 3, rep = 20, standardize = TRUE)
table(raw = km$solution$IDX, standardized = km_s$solution$IDX)