---
title: GLMM Fitting and Model Comparison
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{GLMM Fitting and Model Comparison}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
has_lme4 <- requireNamespace("lme4", quietly = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  eval = has_lme4
)
```

```{r load-package, eval = TRUE}
library(mixeff)
```

```{r lme4-missing, eval = !has_lme4, echo = FALSE, results = "asis"}
cat("The executable cbpp examples in this vignette require the `lme4` package.\n")
```

This vignette shows the current GLMM contract in `mixeff`. The default path
is `glmm(..., method = "pirls_profiled")`, which delegates to the upstream
profiled PIRLS GLMM fitter. `method = "joint_laplace"` requests the labelled
joint route (`fast = FALSE`, `nAGQ = 1`) and is available in the dependency-light
vendored build through the native optimizer.

The example uses the `lme4::cbpp` data. `glmm()` does not yet expose binomial
case weights, so the grouped counts are expanded to Bernoulli rows before
fitting.

```{r cbpp-data}
env <- new.env(parent = emptyenv())
utils::data("cbpp", package = "lme4", envir = env)
cbpp <- get("cbpp", envir = env, inherits = FALSE)

expand_cbpp <- function(data) {
  rows <- lapply(seq_len(nrow(data)), function(i) {
    successes <- data$incidence[[i]]
    failures <- data$size[[i]] - successes
    data.frame(
      herd = data$herd[[i]],
      period = data$period[[i]],
      y = c(rep.int(1L, successes), rep.int(0L, failures))
    )
  })
  droplevels(do.call(rbind, rows))
}

cbpp_binary <- expand_cbpp(cbpp)
c(n_grouped_rows = nrow(cbpp), n_binary_rows = nrow(cbpp_binary))
```

## Family and Audit

The statistical intent is a binomial-logit GLMM with a herd random intercept.
Before fitting, the same formula can be compiled and explained so the fixed and
random-effect structure is visible before numerical optimization starts.

```{r audit-first}
cbpp_family <- binomial(link = "logit")
cbpp_formula <- y ~ period + (1 | herd)

cbpp_spec <- compile_model(cbpp_formula, cbpp_binary)
explain_model(cbpp_spec)
```

## Fit

```{r fit-laplace}
glmm_fit <- glmm(
  cbpp_formula,
  cbpp_binary,
  family = cbpp_family,
  method = "pirls_profiled",
  nAGQ = 1L,
  control = mm_control(verbose = -1)
)

glmm_fit
fixef(glmm_fit)
VarCorr(glmm_fit)
```

The summary prints fitted quantities and deliberately avoids coefficient tests
unless the inference contract certifies them.

```{r summary}
summary(glmm_fit)
```

Core extractors are available from the durable R object.

```{r extractors}
head(fitted(glmm_fit))
head(residuals(glmm_fit))
head(ranef(glmm_fit)[[1L]])
c(
  logLik = as.numeric(logLik(glmm_fit)),
  deviance = deviance(glmm_fit),
  AIC = AIC(glmm_fit),
  BIC = BIC(glmm_fit)
)
```

## Quadrature Sensitivity

`nAGQ` is part of the fit request and is recorded on the object. Values above
one are a profiled-path sensitivity check; the labelled joint-Laplace path is
restricted to `nAGQ = 1`.

```{r agq-sensitivity}
glmm_fit_agq3 <- glmm(
  cbpp_formula,
  cbpp_binary,
  family = cbpp_family,
  method = "pirls_profiled",
  nAGQ = 3L,
  control = mm_control(verbose = -1)
)

data.frame(
  nAGQ = c(glmm_fit$nAGQ, glmm_fit_agq3$nAGQ),
  logLik = c(as.numeric(logLik(glmm_fit)),
             as.numeric(logLik(glmm_fit_agq3))),
  AIC = c(AIC(glmm_fit), AIC(glmm_fit_agq3)),
  check.names = FALSE
)
```

## Current Boundaries

In-sample `fitted()` and `residuals()` are available, and
`summary(fit, tests = "coefficients")` returns a Wald-z fixed-effect
table built from the upstream PIRLS/Laplace working-Hessian covariance
payload (reliability label: `moderate`). GLMM `predict()` (population and
conditional, on the link or response scale), `contrast()`, `drop1()`, and
multi-model `anova()` are available as R-side computations on top of that
certified covariance. `confint()` (Wald) requires a joint-Laplace fit: the
fast-PIRLS/profiled covariance is a working-Hessian payload and the engine
does not certify Wald intervals for it. `simulate()` and `refit()` for GLMMs
are not yet implemented and are exposed as ordinary R conditions rather than
silently filled in. The joint-Laplace fit route is available and labelled
separately from the profiled fast-PIRLS default.

```{r boundaries}
glmm_boundary <- function(expr) {
  cnd <- tryCatch({
    force(expr)
    NULL
  }, error = function(cnd) cnd)
  if (is.null(cnd)) {
    return(data.frame(status = "available", class = NA_character_,
                      message = "", check.names = FALSE))
  }
  data.frame(
    status = "unavailable",
    class = class(cnd)[[1L]],
    message = conditionMessage(cnd),
    check.names = FALSE
  )
}

rbind(
  predict = glmm_boundary(predict(glmm_fit)),
  confint = glmm_boundary(confint(glmm_fit)),
  coefficient_tests = glmm_boundary(summary(glmm_fit, tests = "coefficients")),
  simulate = glmm_boundary(stats::simulate(glmm_fit, nsim = 1L)),
  refit = glmm_boundary(refit(glmm_fit, fitted(glmm_fit))),
  joint_laplace = glmm_boundary(glmm(
    cbpp_formula,
    cbpp_binary,
    family = cbpp_family,
    method = "joint_laplace",
    control = mm_control(verbose = -1)
  ))
)
```

```{r checks, include = FALSE}
stopifnot(inherits(glmm_fit, "mm_glmm"))
stopifnot(identical(glmm_fit$nAGQ, 1L))
stopifnot(identical(glmm_fit_agq3$nAGQ, 3L))
stopifnot(nobs(glmm_fit) == nrow(cbpp_binary))
stopifnot(all(is.finite(fixef(glmm_fit))))
stopifnot(is.finite(as.numeric(logLik(glmm_fit))))
stopifnot(is.finite(AIC(glmm_fit)))
stopifnot(inherits(summary(glmm_fit), "summary.mm_glmm"))
stopifnot(length(fitted(glmm_fit)) == nrow(cbpp_binary))
stopifnot(length(residuals(glmm_fit)) == nrow(cbpp_binary))

boundaries <- rbind(
  predict = glmm_boundary(predict(glmm_fit)),
  confint = glmm_boundary(confint(glmm_fit)),
  coefficient_tests = glmm_boundary(summary(glmm_fit, tests = "coefficients")),
  simulate = glmm_boundary(stats::simulate(glmm_fit, nsim = 1L)),
  joint_laplace = glmm_boundary(glmm(
    cbpp_formula,
    cbpp_binary,
    family = cbpp_family,
    method = "joint_laplace",
    control = mm_control(verbose = -1)
  ))
)
# predict / coefficient tests are available for PIRLS-profiled fits;
# confint requires joint-Laplace (PIRLS working-Hessian not certified for Wald);
# simulate remains unavailable; joint-Laplace route is available.
stopifnot(identical(boundaries["predict", "status"], "available"))
stopifnot(identical(boundaries["confint", "status"], "unavailable"))
stopifnot(identical(boundaries["coefficient_tests", "status"], "available"))
stopifnot(identical(boundaries["simulate", "status"], "unavailable"))
stopifnot(identical(boundaries["joint_laplace", "status"], "available"))
```