---
title: "Estimation engines"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimation engines}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(intraclass)
```

Most of the time you never think about the engine: the default just works, and
every example in the other articles uses it. This article is for when you want to
know what that default is doing, or you have a reason to switch. An
[**engine**](glossary.html#engine) is the computational backend `icc()` uses to
estimate the [variance components](glossary.html#variance-component), chosen with
the `engine` argument. Some engine choices are purely computational: the same
estimator, a different solver. Others compute a genuinely different estimator, one
that is asymptotically equivalent (it converges to the same answer as the sample
grows), or move to a fully Bayesian fit. This article covers the mixed-model
engines (glmmTMB, lme4), the structural-equation engine (lavaan), and the Bayesian
engine (brms), and when the distinction matters. Any unfamiliar term is defined in
the [*Glossary*](glossary.html).

## The mixed-model engines: glmmTMB and lme4

By default `icc()` fits the variance components with **glmmTMB**. You can instead
request **lme4** with `engine = "lme4"` for the random two-way design. The lme4
package itself arrives with the installation, because glmmTMB names lme4 in its
own `Imports:`. But the engine needs one piece more: **merDeriv**. It supplies the
parameter covariance the default Monte-Carlo interval is built from, and every
lme4 fit checks for it on entry whatever interval method you ask for. merDeriv
sits in `Suggests:` and is fetched only on request, so a plain install leaves you
with the lme4 package but not the lme4 engine. Both engines are
[REML](glossary.html#reml) mixed-model fits of the same model. (REML, restricted
maximum likelihood, is the standard way to estimate variance components.) On a
given dataset they therefore return the same coefficients to numerical tolerance:
the choice is about the fitting backend, not the estimand.

```{r engine, eval = requireNamespace("lme4", quietly = TRUE) && requireNamespace("merDeriv", quietly = TRUE)}
glmmtmb <- tidy(icc(ratings, score, subject, rater, engine = "glmmTMB", seed = 1))
lme4 <- tidy(icc(ratings, score, subject, rater, engine = "lme4", seed = 1))
data.frame(
  term = glmmtmb$term,
  glmmTMB = round(glmmtmb$estimate, 4),
  lme4 = round(lme4$estimate, 4)
)
```

The two point estimates agree to well within rounding, and their Monte-Carlo
intervals coincide to about `0.01`. The lme4 interval is built from the parameter
covariance supplied by the **merDeriv** package, transformed onto the same
boundary-aware log-scale glmmTMB uses. glmmTMB remains the recommended default. It
is the engine this package declares in `Imports:`, and it is robust when a variance
component sits exactly at the [zero boundary](glossary.html#zero-variance-boundary).
At that boundary the lme4 route cannot form an interval, and directs you back to
glmmTMB. lme4 otherwise has full design parity with glmmTMB, covering the
fixed-rater and every multilevel design, on balanced **and** incomplete/ragged
data. It degrades to glmmTMB only at that variance boundary.

## A structural-equation engine (`lavaan`)

`engine = "lavaan"` fits the same design as a **structural equation model**, a
common-factor generalizability model in the sense of Jorgensen (2021), for the
random two-way design. Unlike lme4, this is not just a different backend for the
*same* estimator: it matters which coefficient you ask for.

```{r lavaan, eval = requireNamespace("lavaan", quietly = TRUE)}
glmmtmb <- tidy(icc(ratings, score, subject, rater, engine = "glmmTMB", seed = 1))
lavaan <- tidy(icc(ratings, score, subject, rater, engine = "lavaan", seed = 1))
data.frame(
  term = glmmtmb$term,
  glmmTMB = round(glmmtmb$estimate, 4),
  lavaan = round(lavaan$estimate, 4)
)
```

**Consistency** coefficients are a ratio of the subject and residual variances, so
the SEM returns them identically to the mixed model. **Absolute agreement** is
different. The SEM has no random rater effect to estimate: a rater is a single
column, so its effect lives in the column *means*. Following Jorgensen (2021), the
rater variance is recovered from the mean structure as the variance of the estimated
indicator intercepts. This [**indicator-mean
estimator**](glossary.html#indicator-mean-estimator) is a genuinely different
estimator of the rater variance than the mixed model's random effect. The two are
asymptotically equivalent and match conventional generalizability-theory software
(GENOVA, and formerly the `gtheory` package) closely on real data
[Vispoel et al. 2022]. But on a small design they differ by a modest amount: here
`ICC(A,1)` is about `0.284` from lavaan versus `0.290` from the mixed model. The gap
here comes from the raw variance of only four estimated rater means, which carries
small-sample noise the mixed model shrinks away.

Which is "right"? Neither is wrong: they are two defensible estimators of the same
population quantity. Use `"glmmTMB"` (the default) if you want the mixed-model
random-rater estimate and its wider, generalize-to-new-raters interval. Reach for
`"lavaan"` if you are working inside an SEM generalizability-theory workflow and want
results comparable to that literature. The SEM engine covers the random **and
fixed-rater** two-way design, on complete **and incomplete**
([FIML](glossary.html#fiml)) data, with the parametric bootstrap available on
complete data and Monte-Carlo throughout. It also covers the crossed (Design 1)
**multilevel** design as a two-level SEM, reporting the subject- and cluster-level
ICCs off one five-component fit. That route takes random raters, on balanced data
and on incomplete or unbalanced data alike, and gives a Monte-Carlo interval. Two-level
SEM estimation is full-information ML with no REML analog. So with few clusters its
cluster-level components sit slightly below the REML estimates, and its agreement
rater term slightly above. Both are documented differences that shrink as the
cluster count grows, while consistency ICCs agree essentially exactly. One-way
designs are still directed to the mixed-model engines.

## A Bayesian engine (`brms`)

`engine = "brms"` fits the variance components in a fully **Bayesian** framework
(Stan, via the **brms** package). This is the approach ten Hove, Jorgensen & van der
Ark (2020) developed for interrater reliability. Instead of a single REML point with
a Monte-Carlo interval, it samples the posterior of every variance component and
reads the ICC off the draws. The point estimate is the [posterior **mode**
(MAP)](glossary.html#posterior-mode-map) and the interval is a **credible**
interval, covered in [*Confidence-interval
methods*](interval-methods.html#bayesian-credible-intervals-ci_method-posterior).
Because it samples a Stan model, this engine is slower than the others and needs the
`brms` package (an optional `Suggests` dependency).

The examples below are shown with pre-computed output: fitting a Stan model needs a
toolchain not available when this site is built, so these chunks are not evaluated at
knit time. What they show is `icc()`'s own printed output. On a design this small the
sampler may additionally report divergent transitions, and whether it does depends on
the seed and the platform, so those warnings are not reproduced here.

```{r brms, eval = FALSE}
bayes <- icc(ratings, score, subject, rater, engine = "brms", type = "agreement", seed = 1)
bayes
```

```
#> ── Intraclass correlation: two-way random, absolute agreement ──────────────────
#> Subjects: 6 | Raters: 4 (random) | Observations: 24 of 24 cells (complete)
#> Engine: brms (MCMC) | CI: 95% posterior credible (4000 draws)
#>
#>   index     estimate   95% CI
#>   ICC(A,1)     0.241   [0.066, 0.649]
#>   ICC(A,k)     0.679   [0.221, 0.881]
#>
#> Variance components: subject 1.522, rater 2.653, residual 0.962
#> Shrout & Fleiss equivalent: ICC(A,1) = ICC(2,1), ICC(A,k) = ICC(2,k)
```

The header names a **brms (MCMC)** engine and a **posterior credible** interval, and
`ci_method = "posterior"` is automatic. On this tiny six-subject design the MAP
`ICC(A,1)` is about `0.24`, a little below the glmmTMB REML value of `0.29`. The MAP
is the mode of a wide, right-skewed posterior, which the small sample pulls down (ten
Hove et al. 2020 note this at small rater counts). Sampler settings pass through
`brm_args`: chains, iterations, backend, and parallel `cores`, as in
`brm_args = list(chains = 4, cores = 4)`.

### The prior, and overriding it

By default the engine places a weakly-informative **half-*t*(4, 0, 1)**
[prior](glossary.html#prior) on every random-effect standard deviation. That is the
*sourced* prior (ten Hove et al. 2020, §3.3/§4.1) that every coverage result in this
package depends on. You can supply your own prior for prior-sensitivity or
method-comparison work with the `prior` argument, which takes any brms prior object.
But `icc()` warns loudly, because leaving the sourced prior **voids the coverage
guarantees** and a poorly chosen SD prior can *worsen* the small-sample boundary
behavior:

```{r brms-prior, eval = FALSE}
library(brms)
icc(ratings, score, subject, rater, engine = "brms",
  prior = set_prior("normal(0, 0.1)", class = "sd"), seed = 1)
#> Warning message:
#> Using a custom `prior` instead of the sourced half-t(4, 0, 1).
#> ! This VOIDS the package's coverage guarantees: the credible-interval coverage
#>   results (ten Hove et al. 2020) hold only for the sourced prior.
#> ℹ A vague or flat SD prior can WORSEN small-`k` boundary bias -- the half-t is
#>   weakly informative on purpose (Principle #3's regime).
#> ℹ Leave `prior` unset for the sourced default unless you are running
#>   prior-sensitivity or method-comparison work.
```

Here the deliberately over-tight `normal(0, 0.1)` prior squeezes every standard
deviation toward zero, collapsing the ICC to nearly nothing. It is a vivid reminder
that the prior is load-bearing, not a casual knob. Leave `prior` unset unless you
specifically intend to depart from the sourced default. The brms engine covers the
same design family as the mixed-model engines (two-way random and fixed, one-way,
and the multilevel designs), on balanced and incomplete data.