---
title: "Advanced workflows: dispersion, cutpoints and model comparison"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Advanced workflows: dispersion, cutpoints and model comparison}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  fig.align = "center"
)
```

This vignette works through three modeling questions: whether the
dispersion depends on covariates, whether the border scores need their own
cutpoints, and how to compare models that differ in link, interval map or
model family. The model and the notation are those of
`vignette("brs-intro")`. The data are simulated and small, so the results
are illustrative.

```{r library}
library(betaregscale)
```

## Variable dispersion

A two-part formula `y ~ x1 + x2 | z1` models the dispersion with its own
linear predictor, $h(\sigma_i) = z_i^\top \gamma$ (logit link by default
under `repar = 2`). `brs_sim()` takes the dispersion coefficients in
`zeta`, on the link scale:

```{r var-sim}
set.seed(11)
d <- data.frame(x1 = rnorm(300), x2 = rbinom(300, 1, 0.5), z1 = rnorm(300))
sv <- brs_sim(~ x1 + x2 | z1, data = d, beta = c(-0.2, 0.5, -0.5),
              zeta = c(qlogis(0.15), 0.6), ncuts = 10)
table(sv$y)
```

The constant-dispersion model is nested in the variable one, so
`anova()` gives a likelihood-ratio test on the added dispersion
coefficients:

```{r var-fit}
fit_const <- brs(y ~ x1 + x2, data = sv)
fit_var <- brs(y ~ x1 + x2 | z1, data = sv)
anova(fit_const, fit_var)
summary(fit_var)
```

The dispersion coefficients carry the `(phi)_` prefix;
`coef(fit, model = "precision")` selects them, and
`predict(type = "precision")` returns $\sigma_i$:

```{r var-coef}
coef(fit_var, model = "precision")
nd <- data.frame(x1 = 0, x2 = 0, z1 = c(-1, 0, 1))
cbind(nd, sigma = predict(fit_var, nd, type = "precision"),
      sd_latent = 10 * sqrt(predict(fit_var, nd, type = "variance")))
```

`brs_marginaleffects(fit, model = "precision")` gives average marginal
effects on the dispersion submodel.

## Boundary cutpoints

Rating data often pile up at 0 and at $K$. With
`cutpoints = "boundary"` the cutpoints between scores 0 and 1 and between
$K - 1$ and $K$ are estimated (see `vignette("brs-intro")` for the
parameterization and the requirements). Scores simulated with border
cutpoints 0.12 and 0.88 instead of the 0.05 and 0.95 of the `"mid"` map:

```{r cut-sim}
set.seed(1)
heap <- brs_sim(~ x1 + x2, data = d, beta = c(0, 0.5, -0.6),
                phi = qlogis(0.2), ncuts = 10, cutpoints = c(0.12, 0.88))
table(heap$y)
```

The fixed map is nested at an interior point of the boundary model, so
the likelihood-ratio test has 2 degrees of freedom ($K \ge 4$):

```{r cut-fit}
fit_fix <- brs(y ~ x1 + x2, data = heap)
fit_bnd <- brs(y ~ x1 + x2, data = heap, cutpoints = "boundary")
anova(fit_fix, fit_bnd)
```

`summary()` reports the estimated cutpoints on $(0, 1)$ and on the score
scale, with delta-method standard errors, next to their fixed values;
`coef()`, `vcov()` and `confint()` accept `model = "cutpoints"`:

```{r cut-summary}
summary(fit_bnd)
confint(fit_bnd, model = "cutpoints")
```

The average predicted score distribution against the observed one shows
what the two cutpoints change:

```{r cut-pmf}
obs <- prop.table(table(factor(heap$y, levels = 0:10)))
round(cbind(observed = as.numeric(obs),
            fixed = colMeans(brs_predict_scoreprob(fit_fix)),
            boundary = colMeans(brs_predict_scoreprob(fit_bnd))), 3)
```

Cross-validation refits the cutpoints in each training fold and scores
the held-out rows with them. A fold without one of the scores 0, 1,
$K - 1$ or $K$ cannot be fitted and is reported in the `error` column.

```{r cut-cv}
set.seed(2)
cv_fix <- brs_cv(y ~ x1 + x2, data = heap, k = 5)
set.seed(2)
cv_bnd <- brs_cv(y ~ x1 + x2, data = heap, k = 5, cutpoints = "boundary")
c(fixed = mean(cv_fix$log_score), boundary = mean(cv_bnd$log_score))
```

Pseudo $R^2$ and the midpoint residuals use the centers of the estimated
border cells, so they are not directly comparable between fixed-map and
boundary fits; compare those fits by likelihood or held-out log score.

## Model comparison

### Links and interval maps

Fits with different links for the same scores are not nested; compare
them by AIC or BIC:

```{r links}
fits <- lapply(c(logit = "logit", probit = "probit", cloglog = "cloglog"),
               function(l) brs(y ~ x1 + x2 | z1, data = sv, link = l))
brs_table(models = fits, sort_by = "AIC")
```

The same holds across interval maps: with `lim = 0.5` each map gives a
pmf over $0, \ldots, K$, so `logLik()`, `AIC()` and `BIC()` are
comparable across `interval = "mid"`, `"right"` and `"left"`, but the maps
are not nested and `anova()` refuses a likelihood-ratio test between them
(see `vignette("brs-intro")`).

### Comparison with an ordinal model

Because the likelihood of `brs()` is the probability of the observed
scores, it is on the same scale as the likelihood of an ordinal
regression for the same scores, and the two can be compared by AIC. A
cumulative logit model $P(S \le s) = \mathrm{logit}^{-1}(\theta_s - x^\top
\beta)$, with $K$ free thresholds, can be fitted in a few lines of base R:

```{r ordinal}
X <- model.matrix(~ x1 + x2, data = heap)[, -1]
s <- heap$y
K <- 10
# Negative log-likelihood; thresholds kept ordered by positive increments
nll_ord <- function(par) {
  theta <- cumsum(c(par[1], exp(par[2:K])))
  eta <- drop(X %*% par[-(1:K)])
  cum <- cbind(0, plogis(outer(-eta, theta, "+")), 1)
  -sum(log(cum[cbind(seq_along(s), s + 2)] - cum[cbind(seq_along(s), s + 1)]))
}
start <- c(qlogis(mean(s == 0)), rep(log(0.4), K - 1), 0, 0)
ord <- optim(start, nll_ord, method = "BFGS", control = list(maxit = 1000))
ord$convergence

ll <- c(brs_fixed = as.numeric(logLik(fit_fix)),
        brs_boundary = as.numeric(logLik(fit_bnd)),
        cumulative_logit = -ord$value)
df <- c(length(coef(fit_fix)), length(coef(fit_bnd)), length(ord$par))
data.frame(logLik = ll, df = df, AIC = -2 * ll + 2 * df)
```

The cumulative logit model spends $K$ parameters on the thresholds; the
beta model with boundary cutpoints uses two plus the dispersion. The same
reasoning applies out of sample: `brs_loglik_contrib(fit, newdata)` and
`brsmm_cv()` return held-out log probabilities of the scores that can be
set against those of any model giving a pmf on $0, \ldots, K$, evaluated
on the same rows or folds.

## A Monte Carlo study

A simulation study of the estimator is a loop of `brs_sim()` and `brs()`
over a fixed design. The code below is a template (not run here): logit
links, `interval = "mid"`, two scale lengths, two sample sizes and three
dispersion levels, with bias, root mean squared error and coverage of the
95% Wald interval per cell. Increase `R` for precise coverages; with
`R = 200` the Monte Carlo standard error of a coverage near 0.95 is about
0.015.

```{r mc-design, eval = FALSE}
R <- 200L
beta <- c(0.5, 0.2, -0.7)
grid <- expand.grid(sigma = c(0.08, 0.2, 0.4), n = c(100L, 500L),
                    K = c(10L, 100L))
one_cell <- function(sigma, n, K) {
  set.seed(1000 + n + K + round(100 * sigma))
  d <- data.frame(x1 = rnorm(n), x2 = rbinom(n, 1, 0.5))   # fixed design
  truth <- c(beta, qlogis(sigma))
  est <- se <- matrix(NA_real_, R, 4)
  for (r in seq_len(R)) {
    s <- brs_sim(~ x1 + x2, data = d, beta = beta, phi = qlogis(sigma),
                 ncuts = K)
    f <- tryCatch(brs(y ~ x1 + x2, data = s), error = function(e) NULL)
    if (is.null(f) || f$convergence != 0L) next
    est[r, ] <- coef(f)
    se[r, ] <- sqrt(diag(vcov(f)))
  }
  ok <- stats::complete.cases(est, se)
  err <- est[ok, ] - rep(truth, each = sum(ok))
  data.frame(sigma, n, K, par = c("b0", "b1", "b2", "logit_sigma"),
             bias = colMeans(err), rmse = sqrt(colMeans(err^2)),
             cover = colMeans(abs(err) <= qnorm(0.975) * se[ok, ]))
}
res <- do.call(rbind, Map(one_cell, grid$sigma, grid$n, grid$K))
```

For mixed models the same loop uses `brsmm_sim()` and `brsmm()`; the
attribute `attr(sim, "truth")$par` holds the true parameters in the order
of `coef()`.

## References

- Lopes, J. E. (2023). *Modelos de regressão beta para dados de escala*.
  Master's dissertation, Universidade Federal do Paraná, Curitiba.
  <https://acervodigital.ufpr.br/xmlui/handle/1884/86624>
- McCullagh, P. (1980). Regression models for ordinal data. *Journal of
  the Royal Statistical Society: Series B*, 42(2), 109–142.