---
title: Why mixeff?
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Why mixeff?}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
```

```{r load-package}
library(mixeff)
```

If you fit mixed-effects models in R, you most likely use
`lme4::lmer()`. `mixeff::lmm()` aims to be *functionally equivalent*:
the formula language is the same, and the extractors —
`fixef()`, `ranef()`, `VarCorr()`, `predict()`, `simulate()`,
`anova()`, `summary()`, `update()` — do what you expect. Statistical
answers agree with `lme4` within documented tolerances on the parity
datasets shipped with the package. It is not a literal *drop-in*: you
call `lmm()` rather than `lmer()`, results are not bit-exact, and the
package is audit-first — it reports or refuses rather than silently
transforming a model.

The reason to switch is what `mixeff` does around the fit. It is
faster on most designs — typically by a factor of two to five — and
it makes four things explicit that `lme4` leaves implicit:

1.  **The formula stays familiar.** Anything you would hand to
    `lmer()` you can hand to `lmm()`.
2.  **Singular fits become labelled facts.** A reduced-rank
    random-effect covariance is reported with codes, severity, and
    effective rank, instead of a single warning that scrolls off the
    screen.
3.  **Inference labels what asymptotics can and cannot do.** Each
    p-value carries the method that produced it; where Wald,
    Satterthwaite, and Kenward-Roger are unavailable in principle,
    a parametric bootstrap is offered as a labelled first-class peer.
4.  **The fitted object is a record.** A model saved today reopens
    six months later with the same coefficients, the same audit
    trail, and the same method labels — without depending on the
    original Rust handle.

This page demonstrates each of the four on one small dataset.

## The dataset

A small repeated-measures study: 18 subjects, 10 daily reaction-time
measurements each, intercepts and slopes that are nearly perfectly
correlated by construction. This is the kind of design that produces a
singular fit in any modern engine.

```{r make-data}
set.seed(3)
n_subj <- 18L
days   <- 0:9
b0     <- rnorm(n_subj, sd = 30)
b1     <- 0.5 * b0 / 30 * 10 + rnorm(n_subj, sd = 0.5)

sleep_like <- do.call(rbind, lapply(seq_len(n_subj), function(i) {
  data.frame(
    subj = factor(i),
    days = days,
    rt   = 250 + b0[i] + (10 + b1[i]) * days +
           rnorm(length(days), sd = 20)
  )
}))

head(sleep_like, 4)
```

## A. The formula stays familiar

`mixeff` keeps the lme4 random-effects syntax. If you can read `(x | g)`, you
can read `mixeff`.

```{r fit-model}
fit <- lmm(
  rt ~ days + (1 + days | subj),
  sleep_like,
  control = mm_control(verbose = -1)
)
fit
```

```{r familiar-checks, include = FALSE}
stopifnot(inherits(fit, "mm_lmm"))
stopifnot(all(is.finite(fixef(fit))))
stopifnot(length(fitted(fit)) == nrow(sleep_like))
```

`fixef()`, `ranef()`, `VarCorr()`, `predict()`, `simulate()`, `anova()`, and
`summary()` all do what you expect.

## B. When a fit is degenerate, you find out *which* part

`lme4::lmer()` fits this model and reports the situation in one
parsimonious line: `boundary (singular) fit`. The fact is correct.
What it leaves implicit is *which* variance component reached the
boundary, what the effective rank of the random-effect covariance
is, and which downstream inference methods are no longer defined.

```{r lme4-side, eval = requireNamespace("lme4", quietly = TRUE)}
m <- suppressMessages(lme4::lmer(
  rt ~ days + (1 + days | subj),
  data = sleep_like
))
m
lme4::isSingular(m)
```

`mixeff` reports the same fact, and then unpacks it. `fit_status()` names
the convergence outcome, `changes()` shows the requested-to-effective
transition, and `diagnostics()` returns stable codes.

```{r mixeff-side}
fit_status(fit)
is_singular(fit)
changes(fit)
diagnostics(fit)$table[, c("code", "severity", "stage", "message")]
```

```{r diag-checks, include = FALSE}
stopifnot(identical(fit_status(fit), "converged_reduced_rank"))
stopifnot(isTRUE(is_singular(fit)))
diag_codes <- diagnostics(fit)$table$code
stopifnot("covariance_reduced" %in% diag_codes)
ch <- changes(fit)$table
stopifnot("design_time_reduction" %in% ch$stage ||
          "covariance_transition" %in% ch$stage)
```

A reduced-rank covariance is now a *labelled* fact about the fit, not a
warning that scrolls off the screen.

## C. Expose when asymptotics are weak, then offer bootstrap

`mixeff` has two jobs on the inference side. First, it tells you when
the asymptotic methods (Wald z, Satterthwaite, Kenward-Roger) are
unreliable on your fit. Then, when bootstrap is the defensible route,
it offers it as a labelled, first-class peer to the asymptotic methods
— not a fallback trick.

`inference_options()` is the audit verb for that judgment. It
enumerates the inference methods available on the current fit, gives
each one a [closed-enum *reason*](inference-method-glossary.html)
for its status, and names the verb
you would call to invoke it.

```{r inference-options}
opt <- inference_options(fit, "days", nsim = 200)
opt$table[, c("method", "expected_status",
              "expected_reliability_reason", "current")]
```

```{r options-checks, include = FALSE}
stopifnot(inherits(opt, "mm_inference_options"))
satt_row <- opt$table[opt$table$method == "satterthwaite", ]
stopifnot(identical(satt_row$expected_status, "not_assessed"))
stopifnot(identical(satt_row$expected_reliability_reason,
                    "satterthwaite_unavailable_at_boundary"))
```

Two routes are available on this fit: asymptotic Wald z (immediate, but
labelled `low` reliability), and bootstrap (~seconds, labelled by
replicate count and Monte-Carlo SE). Satterthwaite and Kenward-Roger
refuse with a stable reason — `*_unavailable_at_boundary` — because at
a boundary fit the variance-parameter derivative they need is not
defined. That is a fact about the math, not a bug.

The asymptotic Wald row carries its own warrant. `summary()` now
prints `reliability_reason` next to `reliability`:

```{r summary-status}
inf <- inference_table(fit)$table
inf[, c("term", "method", "status", "reliability", "reliability_reason")]
```

```{r summary-status-checks, include = FALSE}
stopifnot("reliability_reason" %in% names(inf))
stopifnot(all(inf$reliability[inf$status == "available"] == "low"))
stopifnot(all(nzchar(inf$reliability_reason)))
```

`degrees_of_freedom_unavailable_so_z_substituted` is the closed-enum
warrant: a t reference distribution was the requested target but the
df could not be computed at this boundary fit, so a standard normal
was substituted. The number is real; the *grade* is calibrated.

For a defensible p-value on this same fit, route through `contrast()`
with `method = "bootstrap"`. The Rust engine simulates from the
constrained null, refits each replicate, and returns a labelled
inference row plus a run payload (boundary rate, MCSE, replicate
count) for audit.

```{r bootstrap-answer}
ct <- contrast(fit, c(0, 1), method = "bootstrap",
               bootstrap = bootstrap_control(nsim = 200, seed = 1))
ct$table[, c("contrast", "estimate", "p_value",
             "method", "status", "reliability")]

run <- ct$table$details[[1]]$bootstrap
data.frame(
  successful_replicates = run$successful_replicates,
  boundary_rate         = round(run$boundary_rate, 3),
  mcse                  = round(run$mcse, 4)
)
```

```{r bootstrap-checks, include = FALSE}
stopifnot(identical(ct$table$status, "available"))
stopifnot(is.finite(ct$table$p_value))
stopifnot(run$successful_replicates >= 100)
```

The bootstrap p-value is `available`, the method is named, the run
payload makes the simulation provenance explicit. The boundary rate is
visible because singular fits propagate boundary behaviour into their
own bootstrap replicates — that is honest, not a bug. `mcse` quantifies
the Monte-Carlo uncertainty of the p-value estimate; raise `nsim` for
a tighter MCSE.

There are three states for any reported quantity: *available with a
named method and a closed-enum warrant*, *unavailable with a stable
reason code*, or *typed error*. There is no fourth state where the
package guesses.

## D. The fit is the record

The fitted object is a serialisable record. `saveRDS()` followed by
`readRDS()` and `revive()` reproduces the audit trail and the extractors
without depending on the original Rust handle.

```{r round-trip}
path <- tempfile(fileext = ".rds")
saveRDS(fit, path)
restored <- revive(readRDS(path))

identical(fixef(restored),               fixef(fit))
identical(changes(restored)$table,       changes(fit)$table)
identical(diagnostics(restored)$table,   diagnostics(fit)$table)
```

```{r round-trip-checks, include = FALSE}
stopifnot(identical(fixef(restored), fixef(fit)))
stopifnot(identical(changes(restored)$table, changes(fit)$table))
stopifnot(identical(diagnostics(restored)$table, diagnostics(fit)$table))
```

A reviewer reading the `.rds` six months from now sees the same
convergence status, the same reduced-rank diagnostic, the same method
labels on the same coefficients.

## What this page did not show

`mixeff` does not (yet, by design):

- match `lme4` numerics bit-for-bit; statistical equivalence within
  documented tolerances on parity datasets is the bar.
- provide the joint-Laplace / AGQ GLMM backend; `glmm()` currently ships
  the labelled profiled-PIRLS path and refuses unavailable joint methods
  explicitly (`vignette("glmm", package = "mixeff")`).
- ship Kenward-Roger or profile-likelihood confidence intervals in v0.

Each of those is a stable boundary with a name, not a missing feature
hidden behind a fallback.

## Where to read next

- `vignette("lmm-basics", package = "mixeff")` — fitting and the standard
  extractors at a slower pace.
- `vignette("inference", package = "mixeff")` — coefficient tests,
  contrasts, term tests, and model comparisons.
- `vignette("demystifying-formulas", package = "mixeff")` — what
  `(1 | g)`, `(x | g)`, split blocks, and `||` actually mean.
- `vignette("saving-and-reviving", package = "mixeff")` — the round-trip
  story in detail.
