## ----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)
xv <- c("informat", "willingn", "officeho", "compete")

## -----------------------------------------------------------------------------
pars <- WD.regression.two.components(M, Yvar = "global", Xvars = xv)
round(pars, 3)

## -----------------------------------------------------------------------------
means <- as.data.frame(get.MatD.stats(M, "mean")$mat)
round(coef(lm(global ~ informat + willingn + officeho + compete, data = means)), 3)

## -----------------------------------------------------------------------------
pred <- WD.regression.two.components.predict(M[, xv], pars)
gof <- WD.regression.GOF(M[, "global"], pred)
gof$RMSE_W
gof$OMEGA
gof$PSEUDOR2$index

## ----fig.height = 5-----------------------------------------------------------
sel <- c("F02_M", "F06_W", "F11_M", "F13_W")
obs <- MatD2long(M[sel, "global"]);   obs$type <- "observed"
prd <- MatD2long(pred[sel, ]);        prd$type <- "predicted"
both <- rbind(obs, prd)
ggplot(both, aes(x, cdf, colour = type)) +
  geom_step(direction = "hv", linewidth = 0.8) +
  facet_wrap(~ row) +
  scale_colour_manual(values = c(observed = "grey20", predicted = "#c62320")) +
  labs(x = "global satisfaction", y = "cumulative distribution", colour = NULL) +
  theme_minimal()

## -----------------------------------------------------------------------------
q_obs <- sapply(c(0.25, 0.5, 0.75), function(p) get.MatD.stats(M[, "global"], "quantile", p)$mat)
q_prd <- sapply(c(0.25, 0.5, 0.75), function(p) get.MatD.stats(pred, "quantile", p)$mat)
df <- data.frame(group = get.MatD.rownames(M), obs = q_obs[, 2], prd = q_prd[, 2],
                 obs_lo = q_obs[, 1], obs_hi = q_obs[, 3])
ggplot(df, aes(prd, obs)) +
  geom_abline(colour = "grey60") +
  geom_errorbar(aes(ymin = obs_lo, ymax = obs_hi), width = 0, colour = "grey60") +
  geom_point(size = 2) +
  labs(x = "predicted median", y = "observed median (bars: observed IQR)") +
  theme_minimal()

## -----------------------------------------------------------------------------
set.seed(42)
B <- 200
boot <- t(replicate(B, {
  idx <- sample(get.MatD.nrows(M), replace = TRUE)
  WD.regression.two.components(M[idx, ], Yvar = "global", Xvars = xv)
}))
ci <- t(apply(boot, 2, quantile, probs = c(0.025, 0.975)))
round(cbind(estimate = pars, ci), 3)