A second opinion for a mixed model

Why

An analyst has just fitted a mixed model – a random intercept for the blocks in a trial, or two crossed intercepts for variety and trial in a multi-environment trial (MET) – and wants an independent second opinion: not another run of the same sampler on the same posterior, but an opinion built from inference methods and model-structure variants that are genuinely different from whatever produced the first fit. Does a heavy-tailed alternative agree with a Gaussian fit? Does a closed-form treatment of the random effects agree with what a sampler-based fit reported? koine answers by stacking over a curated battery of inference lenses in predictive space and returning one synthesised opinion, or by abstaining outright – declining to guess – when the battery cannot offer enough independent perspectives to say anything. This vignette walks the full spine an analyst would use to get that opinion: the request, the battery that resolves against it, a single lens, the diversity that makes a synthesis worth having, and the two mixed-model shapes koine supports – a single random intercept and a multi-environment trial’s two crossed intercepts – ending with what koine does when it should not answer at all.

What

A koine_pass() is an engine-neutral, declarative request: it carries the formula, the data, an evaluable family, the canonical observation axis, and a content hash, and it is the one object every lens and every synthesis call reads. resolve_battery() enrols only the lenses whose family, feature support and tier fit that request, and records why every other lens was excluded, so absence is always explained rather than silent. A lens’s own fit_*() function – fit_dram(), fit_robust_t(), fit_glmm_marginal(), and others in the registry – takes a pass and returns draws, or a typed abstention if the lens cannot field the model. synthesise() is the front door: it resolves the battery, fits every opinion lens, gates on PSIS reliability, and stacks the survivors into one predictive opinion, reporting the stacking weights, the effective number of distinct components, and an over-dispersed envelope that reflects between-lens disagreement as well as within-lens uncertainty. abstention() and is_abstention() make a decline a typed, first-class value rather than an error or a silently wrong number.

Do

set.seed(42)
n_pop <- 40L
x_pop <- rnorm(n_pop)
d_pop <- data.frame(
  y = 1.5 + 0.8 * x_pop + rnorm(n_pop, sd = 0.7), x = x_pop
)
pass <- koine_pass(y ~ x, data = d_pop, family = gaussian())
pass
#> <koine pass>
#>   formula : y ~ x
#>   family  : gaussian / identity
#>   N       : 40
#>   outputs : posterior, log_likelihood
#>   ir_hash : c708f46cca12...

fs <- as_family_spec(gaussian())
fs$loglik_pointwise(y = c(0, 1), eta = c(0, 1), aux = list(sigma = 1))
#> [1] -0.9189385 -0.9189385

resolved <- resolve_battery(pass)
vapply(resolved$enrolled, function(lens) lens$id, character(1))
#> [1] "dram"         "nimble_gibbs" "dream_zs"     "robust_t"     "horseshoe"   
#> [6] "lasso"        "ridge"

fit_pop <- fit_dram(pass, n_iter = 1500L, n_burn = 750L, n_chains = 2L)
fit_pop
#> <koine fb_fit>
#>   status   : contributed
#>   method   : dram  exactness: exact
#>   family/kernel: random_walk / dram
#>   draws    : 750 x 2 x 3 (chain x draw x var)
#>   log_lik  : NULL (unassembled)
#>   maturity : koine-implemented-validated

has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
has_lme4 <- requireNamespace("lme4", quietly = TRUE)

set.seed(19)
n_ht <- 60L
x_ht <- rnorm(n_ht)
y_ht <- 1 + 0.8 * x_ht + rt(n_ht, df = 5) * 0.6      # moderately heavy-tailed
pass_ht <- koine_pass(y_ht ~ x_ht, data.frame(y_ht = y_ht, x_ht = x_ht),
                      gaussian())

reg <- koine_stage1_registry()
lens_id <- function(id) which(vapply(reg, function(z) z$id, "") == id)
reg_ht <- list(reg[[lens_id("dram")]], reg[[lens_id("robust_t")]])
reg_ht[[1]]$fit_fn <- function(p, ls) fit_dram(p, ls, 1500L, 750L, 2L)
reg_ht[[2]]$fit_fn <- function(p, ls) {
  fit_robust_t(p, ls, nu = 4, n_iter = 1500L, n_burn = 750L, n_chains = 2L)
}
res_ht <- synthesise(pass_ht, registry = reg_ht)
res_ht
#> <koine result --- synthesised fourth opinion>
#>   components: 2
#>   weights   :
#>     dram           0.449
#>     robust_t       0.551
#>   elpd (stacked LOO): -65.83
#>   PSIS k threshold  : 0.685
#>   diversity         : 1.98 effective lens(es)
#>   pool envelope     : 2 lenses, 0 disagreeing on >=1 obs (95%, 2x inflated)
#>   evidence axis     : 1 lens(es), log-ml [prior-conditional]

Figure: the stacking weights

weight_df <- data.frame(
  lens = names(w_ht), weight = as.numeric(w_ht), stringsAsFactors = FALSE
)
ggplot2::ggplot(
  weight_df, ggplot2::aes(x = stats::reorder(lens, weight), y = weight,
                          fill = lens)
) +
  ggplot2::geom_col() +
  ggplot2::coord_flip() +
  ggplot2::scale_fill_viridis_d(guide = "none") +
  ggplot2::labs(
    x = "lens", y = "stacking weight",
    title = sprintf("n_eff = %.2f effective lens(es)", n_eff_ht)
  ) +
  ggplot2::theme_minimal()

Stacking weights for the two-lens heavy-tailed battery, with the effective-component count n_eff annotated. Neither lens is at zero, so the mixture is a genuine blend, not a corroboration collapse onto one structure.

Figure: the pooled predictive envelope

ord_ht <- order(x_ht)
points_df <- data.frame(x = x_ht[ord_ht], y = y_ht[ord_ht])
band_df <- data.frame(
  x = x_ht[ord_ht], lower = lpe_ht$lower[ord_ht], upper = lpe_ht$upper[ord_ht],
  pool_mean = lpe_ht$pool_mean[ord_ht]
)
ggplot2::ggplot() +
  ggplot2::geom_ribbon(
    data = band_df, ggplot2::aes(x = x, ymin = lower, ymax = upper),
    fill = "#0072B2", alpha = 0.25
  ) +
  ggplot2::geom_line(
    data = band_df, ggplot2::aes(x = x, y = pool_mean), colour = "#0072B2",
    linewidth = 1
  ) +
  ggplot2::geom_point(
    data = points_df, ggplot2::aes(x = x, y = y), colour = "grey35", size = 1.2
  ) +
  ggplot2::labs(
    x = "x", y = "y",
    title = "Pool envelope: pooled mean response and its 95% band"
  ) +
  ggplot2::theme_minimal()

The over-dispersed linear-pool envelope for the heavy-tailed battery: the pooled posterior-mean response (line) and its 95% envelope (band), built from the two lenses’ means and their between-lens disagreement, against the fitted x range.

set.seed(20)
n_re <- 60L
n_groups_re <- 10L
g_re <- factor(rep(seq_len(n_groups_re), each = n_re / n_groups_re))
x_re <- rnorm(n_re)
b_re <- rnorm(n_groups_re, 0, 1.1)                  # group random intercepts
y_re <- 1 + 0.5 * x_re + b_re[g_re] + rnorm(n_re, 0, 0.7)
pass_re <- koine_pass(
  y_re ~ x_re + (1 | g_re), data.frame(y_re = y_re, x_re = x_re, g_re = g_re),
  gaussian()
)
fit_re <- fit_glmm_marginal(pass_re, n_iter = 1500L, n_burn = 750L,
                            n_chains = 2L)

set.seed(404)
n_variety <- 10L                                    # varieties
n_trial <- 5L                                       # trials (environments)
u_variety <- rnorm(n_variety, 0, 1.0)                # variety effects
v_trial <- rnorm(n_trial, 0, 2.5)                    # trial effects (dominant)
rows_met <- vector("list", n_variety)
for (i in seq_len(n_variety)) {
  trials_i <- sample(seq_len(n_trial), 3L)          # incomplete, unbalanced
  obs_i <- lapply(trials_i, function(t) {
    data.frame(
      variety = i, trial = t,
      yield = 5 + u_variety[i] + v_trial[t] + rnorm(1L, 0, 0.6)
    )
  })
  rows_met[[i]] <- do.call(rbind, obs_i)
} # ends i, over the varieties
met <- do.call(rbind, rows_met)
met$variety <- factor(met$variety)
met$trial <- factor(met$trial)
nrow(met)
#> [1] 30

pass_met <- koine_pass(
  yield ~ 1 + (1 | variety) + (1 | trial), met, gaussian(), seed = 5L
)
pass_met                                            # tagged [MET], crossed
#> <koine pass>
#>   formula : yield ~ 1 + (1 | variety) + (1 | trial)
#>   random  : (1 | variety) + (1 | trial)  [MET]
#>   family  : gaussian / identity
#>   N       : 30
#>   outputs : posterior, log_likelihood
#>   ir_hash : b84c3ccf9161...
fit_met <- fit_glmm_marginal(pass_met, n_iter = 1200L, n_burn = 600L,
                             n_chains = 2L)

prod_fit <- lme4::lmer(
  yield ~ 1 + (1 | variety) + (1 | trial), met, REML = TRUE
)
blup_prod <- lme4::ranef(prod_fit)$variety[, 1]
single_fit <- lme4::lmer(yield ~ 1 + (1 | variety), met, REML = TRUE)
blup_single <- lme4::ranef(single_fit)$variety[, 1]

# koine MET variety effects (GLS BLUP at the posterior-mean variance components)
mv_met <- stats::setNames(sm_met$mean, sm_met$variable)
design_met <- koine:::.pass_design(pass_met)
red_met <- koine:::.pass_re_design(pass_met)
taus_met <- c(mv_met["tau_variety"], mv_met["tau_trial"])
vmat_met <- koine:::.re_marginal_cov(
  design_met$N, red_met$factors, taus_met, mv_met["sigma"]
)
b_names <- paste0("b_", colnames(design_met$X))
resid_met <- design_met$y - as.numeric(design_met$X %*% mv_met[b_names])
z_variety <- koine:::.re_indicator(
  red_met$factors[[1L]]$group, red_met$factors[[1L]]$n_groups
)
u_koine <- as.numeric(
  taus_met[1L]^2 * t(z_variety) %*% solve(vmat_met, resid_met)
)

cor_met <- cor(u_koine, blup_prod)
cor_single <- cor(blup_single, blup_prod)

counts_met <- transform(met, n = rpois(nrow(met), 4))
fit_glmm_marginal(
  koine_pass(n ~ 1 + (1 | variety) + (1 | trial), counts_met, poisson())
)$reason
#> [1] "feature_unsupported"

pass_bad <- koine_pass(
  y_re ~ x_re + (x_re | g_re),
  data.frame(y_re = y_re, x_re = x_re, g_re = g_re), gaussian()
)
synthesise(pass_bad)$reason
#> [1] "feature_unsupported"

out_abstain <- synthesise(pass, registry = list(koine_stage1_registry()[[1]]))
is_abstention(out_abstain)
#> [1] TRUE
out_abstain
#> <koine abstention>
#>   scope : synthesis
#>   reason: no_diverse_battery
#>   detail: only 1 opinion lens(es) resolved; >= 2 independent lenses needed.

Read

The DRAM lens recovers the data-generating intercept and slope closely: 1.57 against a true 1.5, and 0.79 against a true 0.8.

On the heavy-tailed toy, the two-lens battery does not collapse onto a single structure: DRAM carries 0.45 of the stacking weight and robust-t carries 0.55, giving 1.98 effective lens(es) – above the package’s degeneracy threshold of 1.5, so the diagnostic records genuine diversity rather than a corroboration collapse. The weights figure shows that split directly, and the envelope figure shows what it buys: a pooled band wide enough that, on this particular draw, 0 of the 60 points fall outside the 95% band even though the band is inflated for between-lens disagreement as well as within-lens uncertainty.

The single-random-intercept marginal lens recovers the data-generating values in the same qualitative way: intercept 0.9 (true 1), slope 0.61 (true 0.5), residual SD 0.69 (true 0.7) and group SD 1.34 (true 1.1), on 60 observations in 10 groups – fewer groups than a production analysis would carry, which is why the variance component is recovered only qualitatively rather than to two decimal places.

On the multi-environment toy, koine’s MET marginal variety effects correlate 1.0000 with the production (lme4) MET fit, against 0.9280 for the single-random-intercept reduction that ignores the trial structure. The gap between the two reproduces, at this much smaller scale, the qualitative point the section makes: confounding a variety with the trials it happened to be grown in costs real agreement with the model that actually generated the data, even on an incomplete, unbalanced design of only 30 records.

Every abstention in this vignette is grounded in the same registry the battery table used above: the Poisson MET call, the random-slope call, and the single-lens synthesise() call all decline through the typed koine_abstention path (is_abstention() returns TRUE for the last of the three) rather than returning a silently wrong number.

Limits

The random-effects lenses carry a fixed weakly-informative prior – N(0, 5) on every fixed-effect coefficient and lognormal(0, 2) on every variance-component standard deviation for the Gaussian and binomial families – and koine_pass()’s prior argument, though present in the signature, does not reach them: it is reserved for lenses that do not yet exist. On a response measured in the tens or hundreds this reference prior pulls the intercept toward zero and lets the group standard deviation absorb the resulting offset, because the prior’s own scale, not the data’s, starts to dominate the posterior. The examples in this vignette are engineered on roughly unit scale for exactly this reason: they are not a demonstration that the prior is harmless, only that it is invisible when the design already matches its assumption. A caller with a response in its raw, natural units – yield in kilograms per hectare, a count in the hundreds, anything not already standardised – should centre and scale the response and every non-intercept covariate before building the pass, and map any recovered coefficients back to the original scale afterwards; the random-intercept and MET marginal lenses in this vignette would return a materially different, prior-dominated opinion on the same data left in its natural units.

Every fitted example in this vignette uses at most 60 observations and a few hundred to two thousand MCMC iterations over two chains, chosen so the whole vignette renders in under two minutes rather than to reflect the scale koine is meant to run at; the MET section is the one exception, because demonstrating a crossed two-grouping structure at all needs two grouping factors, so it shrinks every other dimension of the design instead. The two-lens diversity split shown above is itself sensitive to this small a sample: at n = 60, most seeds tried while preparing this vignette gave a degenerate, single-lens-dominated split rather than the genuine two-lens blend shown here, whereas the same scenario at a few hundred observations is far less seed-sensitive. A reader who wants a diversity diagnostic that does not depend on which random draw they happened to get should run the same battery at production scale, not at this vignette’s render-time scale.

The MET marginal is exact only for the Gaussian identity-link case; a non-Gaussian MET routes through a marginalised sampler that can itself abstain on convergence grounds, and that sampler is not exercised in this vignette (the Poisson MET example above abstains before any sampling is attempted, because the structure itself, not a convergence failure, is what is declined). koine’s random-effects corner is limited throughout to a single random intercept or the two-grouping MET structure: random slopes, correlated random effects, and more than two crossed or nested groupings are all refused rather than approximated, as the random-slope example above shows.

What to read next

The package’s own INTEGRATION.md (installed alongside the package, system.file("INTEGRATION.md", package = "koine")) specifies the contract a parent triangulation package drives to use koine as a fourth, synthesised witness, including the exact supported-model corner this vignette exercises. The quorum package’s Getting started with quorum vignette shows koine consumed as one lens in a wider quorum of independent estimates, through its emit_koine() wrapper; that wrapper standardises the response and covariates before calling koine, which is the practical answer to the prior-scale limit stated above.

Reproduce

Seeds are set per section (42 for the population regression, 19 for the heavy-tailed diversity toy, 20 for the random-intercept toy, 404 for the MET toy) so each fitted example is reproducible on its own; package versions follow.

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] koine_0.3.1.9000 rmarkdown_2.32  
#> 
#> loaded via a namespace (and not attached):
#>  [1] bridgesampling_1.2-1 tensorA_0.36.2.1     sass_0.4.10         
#>  [4] generics_0.1.4       stringi_1.8.9        lattice_0.23-1      
#>  [7] lme4_2.0-6           digest_0.6.39        magrittr_2.0.5      
#> [10] evaluate_1.0.5       grid_4.6.1           RColorBrewer_1.1-3  
#> [13] mvtnorm_1.4-2        fastmap_1.2.0        jsonlite_2.0.0      
#> [16] Matrix_1.7-6         backports_1.5.1      Brobdingnag_1.3-1   
#> [19] viridisLite_0.4.3    scales_1.4.0         jquerylib_0.1.4     
#> [22] reformulas_0.4.4     abind_1.4-8          Rdpack_2.6.6        
#> [25] cli_3.6.6            rlang_1.3.0          rbibutils_2.4.1     
#> [28] splines_4.6.1        withr_3.0.3          cachem_1.1.0        
#> [31] yaml_2.3.12          otel_0.2.0           parallel_4.6.1      
#> [34] tools_4.6.1          coda_0.19-4.1        nloptr_2.2.1        
#> [37] checkmate_2.3.4      minqa_1.2.8          dplyr_1.2.1         
#> [40] ggplot2_4.0.3        boot_1.3-32          buildtools_1.0.0    
#> [43] vctrs_0.7.3          posterior_1.7.0      R6_2.6.1            
#> [46] matrixStats_1.5.0    lifecycle_1.0.5      stringr_1.6.0       
#> [49] MASS_7.3-66          pkgconfig_2.0.3      pillar_1.11.1       
#> [52] bslib_0.12.0         gtable_0.3.6         loo_2.10.1          
#> [55] glue_1.8.1           Rcpp_1.1.2           xfun_0.61           
#> [58] tibble_3.3.1         tidyselect_1.2.1     sys_3.4.3           
#> [61] knitr_1.52           farver_2.1.2         htmltools_0.5.9     
#> [64] nlme_3.1-171         labeling_0.4.3       maketools_1.3.2     
#> [67] compiler_4.6.1       S7_0.2.2             distributional_0.9.0