From 36c7b282f11b467b1c2f664acdbd2f5d8c1c857e Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 17:49:31 -0400
Subject: [PATCH 01/19] Add mixed models draft metadata

---
 mixed-models/DESCRIPTION | 4 ++++
 1 file changed, 4 insertions(+)
 create mode 100644 mixed-models/DESCRIPTION

diff --git a/mixed-models/DESCRIPTION b/mixed-models/DESCRIPTION
new file mode 100644
index 0000000..18bfcba
--- /dev/null
+++ b/mixed-models/DESCRIPTION
@@ -0,0 +1,4 @@
+Title: Mixed Models Explorer
+Description: Explore grouped data, random effects, partial pooling, shrinkage, and model diagnostics.
+Categories: statistics, simulation, regression
+Status: draft

From 1a5a9e76d3e9078412cd723b43cad9b8f84d8388 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 17:49:46 -0400
Subject: [PATCH 02/19] Document mixed models learning flow

---
 mixed-models/readme.md | 47 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 mixed-models/readme.md

diff --git a/mixed-models/readme.md b/mixed-models/readme.md
new file mode 100644
index 0000000..e2a26a6
--- /dev/null
+++ b/mixed-models/readme.md
@@ -0,0 +1,47 @@
+## How it works
+
+This experiment separates two ideas that are easy to mix up:
+
+1. **Data structure** controls how the simulated data are generated.
+2. **Fitted model** controls how we choose to explain those same data.
+
+The simulated relationship is
+
+\[
+y_{ij} = \beta_0 + b_{0j} + (\beta_1 + b_{1j})x_{ij} + \varepsilon_{ij},
+\]
+
+where the group deviations are centered at zero and the observation noise is normal:
+
+\[
+\begin{pmatrix}b_{0j}\\b_{1j}\end{pmatrix} \sim N(0, \Sigma),
+\qquad
+\varepsilon_{ij} \sim N(0, \sigma^2).
+\]
+
+Depending on the selected data structure, the random-intercept or random-slope variance can be exactly zero.
+
+### Pooling
+
+**Complete pooling** ignores group differences and estimates one relationship for everyone.
+
+**No pooling** estimates a separate relationship for each group.
+
+**Partial pooling** is the mixed-model middle ground: groups have their own effects, but those effects are estimated together and share information.
+
+### Shrinkage
+
+Shrinkage is the visible result of partial pooling. A group estimate is pulled toward the population relationship when its own data are uncertain. In random-effect notation, the estimated group deviation is pulled toward zero because zero means "no deviation from the population effect."
+
+Groups with less information generally shrink more. Groups with more observations, less noise, or stronger evidence of genuine between-group differences generally shrink less.
+
+### What to inspect
+
+- **Main view:** simulated truth and fitted group relationships.
+- **Residuals vs fitted:** remaining structure or changing residual spread.
+- **Normal Q-Q:** whether residuals look compatible with a normal-error assumption.
+- **Random effects:** estimated group deviations around zero.
+- **Shrinkage:** no-pooling estimates compared with the estimates from the selected model.
+- **Variance components:** true simulation standard deviations compared with those estimated by the fitted mixed model.
+
+A singular mixed-model fit is informative here: it often means the fitted random-effect covariance has reached a boundary, commonly because one random-effect variance is estimated close to zero.

From d0f433f21d33b8e0f2f3065af39a52dc0f9d20ac Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 17:49:57 -0400
Subject: [PATCH 03/19] Add mixed models credits

---
 mixed-models/credits.md | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 mixed-models/credits.md

diff --git a/mixed-models/credits.md b/mixed-models/credits.md
new file mode 100644
index 0000000..d2e18a3
--- /dev/null
+++ b/mixed-models/credits.md
@@ -0,0 +1 @@
+App made by [Joshua Kunst](https://jkunst.com) with ❤️ and ☕ using Shiny for R ✨. Code [here](https://github.com/jbkunst/visual-data-lab).

From d9893acb37690b53192e25629f01b42d72f12428 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 17:51:28 -0400
Subject: [PATCH 04/19] Prototype mixed models explorer with base R plots

---
 mixed-models/app.R | 631 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 631 insertions(+)
 create mode 100644 mixed-models/app.R

diff --git a/mixed-models/app.R b/mixed-models/app.R
new file mode 100644
index 0000000..8ac2569
--- /dev/null
+++ b/mixed-models/app.R
@@ -0,0 +1,631 @@
+# packages ----------------------------------------------------------------
+library(shiny)
+library(bslib)
+library(lme4)
+library(vdltheme)
+
+# theme -------------------------------------------------------------------
+apptheme <- theme_vdl()
+
+thematic::thematic_shiny(font = "auto")
+
+sidebar <- purrr::partial(bslib::sidebar, width = 320)
+
+card <- purrr::partial(
+  bslib::card,
+  full_screen = TRUE,
+  wrapper = purrr::partial(bslib::card_body, padding = 0)
+)
+
+primary_color <- unname(bs_get_variables(apptheme, "primary"))
+
+# helpers -----------------------------------------------------------------
+simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
+  set.seed(seed)
+
+  beta0 <- 2
+  beta1 <- 1
+  sigma_b0 <- if (structure %in% c("intercept", "both")) 1.5 else 0
+  sigma_b1 <- if (structure %in% c("slope", "both")) 0.65 else 0
+  sigma_e <- 0.8
+  rho <- if (structure == "both") 0.35 else 0
+
+  groups <- paste0("G", seq_len(n_groups))
+  z0 <- rnorm(n_groups)
+  z1 <- rnorm(n_groups)
+
+  b0 <- sigma_b0 * z0
+  b1 <- sigma_b1 * (rho * z0 + sqrt(1 - rho^2) * z1)
+
+  effects <- data.frame(
+    group = factor(groups, levels = groups),
+    b0 = b0,
+    b1 = b1
+  )
+
+  dat <- do.call(
+    rbind,
+    lapply(seq_len(n_groups), function(i) {
+      x <- sort(runif(n_per_group, -2, 2))
+      epsilon <- rnorm(n_per_group, 0, sigma_e)
+      y <- beta0 + b0[i] + (beta1 + b1[i]) * x + epsilon
+
+      data.frame(
+        group = factor(groups[i], levels = groups),
+        x = x,
+        y = y
+      )
+    })
+  )
+
+  list(
+    data = dat,
+    effects = effects,
+    truth = list(
+      beta0 = beta0,
+      beta1 = beta1,
+      sigma_b0 = sigma_b0,
+      sigma_b1 = sigma_b1,
+      sigma_e = sigma_e,
+      rho = rho
+    )
+  )
+}
+
+fit_selected_model <- function(dat, model) {
+  switch(
+    model,
+    pooled = lm(y ~ x, data = dat),
+    separate = lm(y ~ x * group, data = dat),
+    ri = lmer(y ~ x + (1 | group), data = dat, REML = TRUE),
+    rs = lmer(y ~ x + (0 + x | group), data = dat, REML = TRUE),
+    ris = lmer(y ~ x + (1 + x | group), data = dat, REML = TRUE)
+  )
+}
+
+truth_lines <- function(sim) {
+  groups <- levels(sim$data$group)
+  x_grid <- seq(-2.1, 2.1, length.out = 100)
+
+  do.call(
+    rbind,
+    lapply(groups, function(g) {
+      effect <- sim$effects[sim$effects$group == g, ]
+
+      data.frame(
+        group = factor(g, levels = groups),
+        x = x_grid,
+        y = sim$truth$beta0 + effect$b0 +
+          (sim$truth$beta1 + effect$b1) * x_grid
+      )
+    })
+  )
+}
+
+fitted_lines <- function(dat, fit) {
+  groups <- levels(dat$group)
+  x_grid <- seq(-2.1, 2.1, length.out = 100)
+  grid <- expand.grid(x = x_grid, group = groups)
+  grid$group <- factor(grid$group, levels = groups)
+  grid$y <- predict(fit, newdata = grid)
+  grid
+}
+
+no_pool_coefficients <- function(dat) {
+  fits <- lapply(split(dat, dat$group), function(d) coef(lm(y ~ x, data = d)))
+  out <- do.call(rbind, fits)
+  colnames(out) <- c("intercept", "slope")
+  out
+}
+
+current_group_coefficients <- function(dat, fit, model) {
+  groups <- levels(dat$group)
+
+  if (model == "pooled") {
+    coefs <- coef(fit)[c("(Intercept)", "x")]
+    out <- matrix(rep(coefs, each = length(groups)), ncol = 2)
+    rownames(out) <- groups
+    colnames(out) <- c("intercept", "slope")
+    return(out)
+  }
+
+  if (model == "separate") {
+    return(no_pool_coefficients(dat))
+  }
+
+  out <- as.matrix(coef(fit)$group[, c("(Intercept)", "x"), drop = FALSE])
+  colnames(out) <- c("intercept", "slope")
+  out
+}
+
+population_coefficients <- function(dat, fit, model) {
+  if (inherits(fit, "merMod")) {
+    out <- fixef(fit)[c("(Intercept)", "x")]
+  } else if (model == "pooled") {
+    out <- coef(fit)[c("(Intercept)", "x")]
+  } else {
+    out <- coef(lm(y ~ x, data = dat))[c("(Intercept)", "x")]
+  }
+
+  unname(out)
+}
+
+estimated_sds <- function(fit) {
+  out <- c(intercept = NA_real_, slope = NA_real_, residual = sigma(fit))
+
+  if (!inherits(fit, "merMod")) {
+    return(out)
+  }
+
+  vc <- VarCorr(fit)$group
+  sd_re <- attr(vc, "stddev")
+
+  if ("(Intercept)" %in% names(sd_re)) {
+    out["intercept"] <- unname(sd_re["(Intercept)"])
+  }
+
+  if ("x" %in% names(sd_re)) {
+    out["slope"] <- unname(sd_re["x"])
+  }
+
+  out
+}
+
+model_labels <- c(
+  pooled = "Complete pooling",
+  separate = "No pooling",
+  ri = "Random intercept",
+  rs = "Random slope",
+  ris = "Random intercept + slope"
+)
+
+structure_labels <- c(
+  none = "No group differences",
+  intercept = "Different intercepts",
+  slope = "Different slopes",
+  both = "Different intercepts + slopes"
+)
+
+# ui ----------------------------------------------------------------------
+ui <- page_fillable(
+  theme = apptheme,
+  padding = 0,
+  layout_sidebar(
+    fillable = TRUE,
+    padding = "0.75rem",
+    sidebar = sidebar(
+      title = "Mixed Models Explorer",
+      radioButtons(
+        "data_structure",
+        input_label_vdl(
+          "1. Data structure",
+          "Changes the process that generates the points."
+        ),
+        choices = structure_labels,
+        selected = "both"
+      ),
+      radioButtons(
+        "model",
+        input_label_vdl(
+          "2. Fitted model",
+          "Keeps the data fixed and changes how group structure is modeled."
+        ),
+        choices = model_labels,
+        selected = "ri"
+      ),
+      sliderInput(
+        "n_groups",
+        tags$small("Groups"),
+        min = 4,
+        max = 10,
+        value = 6,
+        step = 1
+      ),
+      sliderInput(
+        "n_per_group",
+        tags$small("Observations per group"),
+        min = 5,
+        max = 40,
+        value = 18,
+        step = 1
+      ),
+      actionButton("resimulate", "Resimulate data", width = "100%"),
+      tags$hr(),
+      tags$small(tags$strong("Underlying model")),
+      uiOutput("truth_formula"),
+      tags$small(tags$strong("Fitted model")),
+      uiOutput("fit_formula"),
+      uiOutput("pooling_note"),
+      uiOutput("model_check"),
+      accordion(
+        open = FALSE,
+        accordion_panel(
+          "How it works",
+          tags$small(htmltools::includeMarkdown("readme.md"))
+        )
+      ),
+      tags$small(htmltools::includeMarkdown("credits.md"))
+    ),
+    layout_columns(
+      col_widths = c(12, 6, 6, 6, 6),
+      gap = "0.75rem",
+      card(
+        card_header(uiOutput("main_title")),
+        plotOutput("main_plot", width = "100%", height = "45vh")
+      ),
+      card(
+        card_header("Residuals vs fitted"),
+        plotOutput("residual_plot", width = "100%", height = "27vh")
+      ),
+      card(
+        card_header("Normal Q-Q"),
+        plotOutput("qq_plot", width = "100%", height = "27vh")
+      ),
+      card(
+        card_header("Random effects"),
+        plotOutput("random_effects_plot", width = "100%", height = "27vh")
+      ),
+      card(
+        card_header("Pooling / shrinkage"),
+        plotOutput("shrinkage_plot", width = "100%", height = "27vh")
+      ),
+      card(
+        card_header("Variance components"),
+        plotOutput("variance_plot", width = "100%", height = "27vh")
+      )
+    )
+  )
+)
+
+# server ------------------------------------------------------------------
+server <- function(input, output, session) {
+  seed <- reactiveVal(100L)
+
+  observeEvent(input$resimulate, {
+    seed(seed() + 1L)
+  })
+
+  sim <- reactive({
+    simulate_grouped_data(
+      structure = input$data_structure,
+      n_groups = input$n_groups,
+      n_per_group = input$n_per_group,
+      seed = seed()
+    )
+  })
+
+  fit <- reactive({
+    fit_selected_model(sim()$data, input$model)
+  })
+
+  output$truth_formula <- renderUI({
+    truth <- sim()$truth
+
+    formula <- switch(
+      input$data_structure,
+      none = "y_{ij} = \\beta_0 + \\beta_1 x_{ij} + \\varepsilon_{ij}",
+      intercept = "y_{ij} = \\beta_0 + b_{0j} + \\beta_1 x_{ij} + \\varepsilon_{ij}",
+      slope = "y_{ij} = \\beta_0 + (\\beta_1 + b_{1j})x_{ij} + \\varepsilon_{ij}",
+      both = "y_{ij} = \\beta_0 + b_{0j} + (\\beta_1 + b_{1j})x_{ij} + \\varepsilon_{ij}"
+    )
+
+    random_assumption <- switch(
+      input$data_structure,
+      none = "",
+      intercept = "<div>\\(b_{0j} \\sim N(0, \\sigma^2_{b0})\\)</div>",
+      slope = "<div>\\(b_{1j} \\sim N(0, \\sigma^2_{b1})\\)</div>",
+      both = "<div>\\(\\mathbf b_j \\sim N(\\mathbf 0, \\Sigma)\\)</div>"
+    )
+
+    withMathJax(
+      HTML(
+        paste0(
+          "<div>\\(", formula, "\\)</div>",
+          random_assumption,
+          "<div>\\(\\varepsilon_{ij} \\sim N(0, \\sigma^2)\\)</div>",
+          sprintf(
+            "<div><small>True SDs: intercept %.2f · slope %.2f · residual %.2f</small></div>",
+            truth$sigma_b0,
+            truth$sigma_b1,
+            truth$sigma_e
+          )
+        )
+      )
+    )
+  })
+
+  output$fit_formula <- renderUI({
+    math <- switch(
+      input$model,
+      pooled = "y_{ij} = \\beta_0 + \\beta_1x_{ij} + \\varepsilon_{ij}",
+      separate = "y_{ij} = \\alpha_{0j} + \\alpha_{1j}x_{ij} + \\varepsilon_{ij}",
+      ri = "y_{ij} = \\beta_0 + b_{0j} + \\beta_1x_{ij} + \\varepsilon_{ij}",
+      rs = "y_{ij} = \\beta_0 + (\\beta_1 + b_{1j})x_{ij} + \\varepsilon_{ij}",
+      ris = "y_{ij} = \\beta_0 + b_{0j} + (\\beta_1 + b_{1j})x_{ij} + \\varepsilon_{ij}"
+    )
+
+    code <- switch(
+      input$model,
+      pooled = "lm(y ~ x, data = dat)",
+      separate = "lm(y ~ x * group, data = dat)",
+      ri = "lmer(y ~ x + (1 | group), data = dat)",
+      rs = "lmer(y ~ x + (0 + x | group), data = dat)",
+      ris = "lmer(y ~ x + (1 + x | group), data = dat)"
+    )
+
+    tagList(
+      withMathJax(HTML(paste0("<div>\\(", math, "\\)</div>"))),
+      tags$code(code)
+    )
+  })
+
+  output$pooling_note <- renderUI({
+    text <- switch(
+      input$model,
+      pooled = "Complete pooling · one relationship is shared by every group.",
+      separate = "No pooling · each group gets its own OLS relationship.",
+      "Partial pooling · group deviations are estimated jointly and shrink toward zero."
+    )
+
+    tags$div(class = "small text-muted mt-2", tags$strong("Information sharing: "), text)
+  })
+
+  output$model_check <- renderUI({
+    mod <- fit()
+    residual_mean <- mean(residuals(mod))
+
+    if (!inherits(mod, "merMod")) {
+      return(
+        tags$div(
+          class = "small mt-2",
+          tags$strong("Model check"),
+          tags$div("Random effects: not modeled"),
+          tags$div(sprintf("Mean residual: %.3f", residual_mean))
+        )
+      )
+    }
+
+    convergence_messages <- mod@optinfo$conv$lme4$messages
+    converged <- is.null(convergence_messages)
+    singular <- isSingular(mod, tol = 1e-4)
+    sds <- estimated_sds(mod)
+
+    tags$div(
+      class = "small mt-2",
+      tags$strong("Model check"),
+      tags$div(if (converged) "✓ optimizer converged" else "⚠ convergence warning"),
+      tags$div(
+        if (singular) {
+          "⚠ singular fit: random-effect covariance is on a boundary"
+        } else {
+          "✓ random-effect covariance is full rank"
+        }
+      ),
+      tags$div(
+        sprintf(
+          "Estimated SDs: intercept %s · slope %s",
+          ifelse(is.na(sds["intercept"]), "—", sprintf("%.3f", sds["intercept"])),
+          ifelse(is.na(sds["slope"]), "—", sprintf("%.3f", sds["slope"]))
+        )
+      ),
+      tags$div(sprintf("Mean residual: %.3f", residual_mean))
+    )
+  })
+
+  output$main_title <- renderUI({
+    tags$span(
+      structure_labels[[input$data_structure]],
+      tags$span(class = "text-muted", "  →  "),
+      model_labels[[input$model]]
+    )
+  })
+
+  output$main_plot <- renderPlot({
+    simulation <- sim()
+    dat <- simulation$data
+    mod <- fit()
+    truth <- truth_lines(simulation)
+    fitted <- fitted_lines(dat, mod)
+    groups <- levels(dat$group)
+    cols <- setNames(hcl.colors(length(groups), "Dark 3"), groups)
+
+    plot(
+      dat$x,
+      dat$y,
+      col = adjustcolor(cols[as.character(dat$group)], alpha.f = 0.65),
+      pch = 16,
+      xlab = "x",
+      ylab = "y",
+      main = ""
+    )
+
+    for (g in groups) {
+      dtruth <- truth[truth$group == g, ]
+      lines(dtruth$x, dtruth$y, col = adjustcolor(cols[g], alpha.f = 0.45), lty = 2, lwd = 2)
+    }
+
+    if (input$model == "pooled") {
+      abline(mod, col = primary_color, lwd = 3)
+    } else {
+      for (g in groups) {
+        dfit <- fitted[fitted$group == g, ]
+        lines(dfit$x, dfit$y, col = cols[g], lwd = 2.5)
+      }
+    }
+
+    legend(
+      "topleft",
+      legend = c("Simulated truth", "Fitted model"),
+      lty = c(2, 1),
+      lwd = c(2, 2.5),
+      bty = "n",
+      cex = 0.85
+    )
+  })
+
+  output$residual_plot <- renderPlot({
+    mod <- fit()
+    x <- fitted(mod)
+    y <- residuals(mod)
+
+    plot(
+      x,
+      y,
+      pch = 16,
+      col = adjustcolor(primary_color, alpha.f = 0.55),
+      xlab = "Fitted",
+      ylab = "Residual",
+      main = ""
+    )
+    abline(h = 0, lty = 2)
+    lines(lowess(x, y), lwd = 2)
+  })
+
+  output$qq_plot <- renderPlot({
+    r <- residuals(fit())
+    qqnorm(r, pch = 16, col = adjustcolor(primary_color, alpha.f = 0.55), main = "")
+    qqline(r, lwd = 2)
+  })
+
+  output$random_effects_plot <- renderPlot({
+    mod <- fit()
+
+    if (!inherits(mod, "merMod")) {
+      plot.new()
+      text(0.5, 0.55, "No random effects in this model", cex = 1.05)
+      text(0.5, 0.43, "Complete/no pooling do not estimate b_j", cex = 0.85)
+      return(invisible())
+    }
+
+    re <- ranef(mod)$group
+    groups <- rownames(re)
+    values <- as.matrix(re)
+    xr <- range(c(0, values))
+    pad <- max(diff(xr) * 0.08, 0.1)
+
+    plot(
+      c(xr[1] - pad, xr[2] + pad),
+      c(0.5, length(groups) + 0.5),
+      type = "n",
+      yaxt = "n",
+      xlab = "Deviation from population effect",
+      ylab = "",
+      main = ""
+    )
+    axis(2, at = seq_along(groups), labels = groups, las = 1, cex.axis = 0.8)
+    abline(v = 0, lty = 2)
+
+    pchs <- c(16, 1)
+    for (j in seq_len(ncol(values))) {
+      points(values[, j], seq_along(groups), pch = pchs[j], cex = 1.15)
+    }
+
+    legend(
+      "topright",
+      legend = colnames(values),
+      pch = pchs[seq_len(ncol(values))],
+      bty = "n",
+      cex = 0.8
+    )
+    mtext(expression(b[j] %~% N(0, Sigma)), side = 3, line = -1.2, adj = 0, cex = 0.8)
+  })
+
+  output$shrinkage_plot <- renderPlot({
+    dat <- sim()$data
+    mod <- fit()
+    groups <- levels(dat$group)
+    no_pool <- no_pool_coefficients(dat)
+    current <- current_group_coefficients(dat, mod, input$model)
+    population <- population_coefficients(dat, mod, input$model)
+
+    xr <- range(c(no_pool[, "intercept"], current[, "intercept"], population[1]))
+    yr <- range(c(no_pool[, "slope"], current[, "slope"], population[2]))
+    xpad <- max(diff(xr) * 0.12, 0.2)
+    ypad <- max(diff(yr) * 0.12, 0.15)
+
+    plot(
+      c(xr[1] - xpad, xr[2] + xpad),
+      c(yr[1] - ypad, yr[2] + ypad),
+      type = "n",
+      xlab = "Group intercept",
+      ylab = "Group slope",
+      main = ""
+    )
+
+    segments(
+      no_pool[, "intercept"],
+      no_pool[, "slope"],
+      current[, "intercept"],
+      current[, "slope"],
+      col = "grey70"
+    )
+    points(no_pool[, "intercept"], no_pool[, "slope"], pch = 1, cex = 1.1)
+    points(current[, "intercept"], current[, "slope"], pch = 16, cex = 1.1)
+    points(population[1], population[2], pch = 8, cex = 1.5, lwd = 2)
+
+    if (input$model %in% c("ri", "rs", "ris")) {
+      text(
+        current[, "intercept"],
+        current[, "slope"],
+        labels = groups,
+        pos = 3,
+        cex = 0.7
+      )
+    }
+
+    legend(
+      "topright",
+      legend = c("No-pooling estimate", "Selected model", "Population effect"),
+      pch = c(1, 16, 8),
+      bty = "n",
+      cex = 0.78
+    )
+  })
+
+  output$variance_plot <- renderPlot({
+    truth <- sim()$truth
+    mod <- fit()
+
+    true_sd <- c(
+      intercept = truth$sigma_b0,
+      slope = truth$sigma_b1,
+      residual = truth$sigma_e
+    )
+    estimated_sd <- estimated_sds(mod)
+
+    values <- c(true_sd, estimated_sd)
+    ymax <- max(values, na.rm = TRUE) * 1.2
+    if (!is.finite(ymax) || ymax == 0) ymax <- 1
+
+    plot(
+      c(0.7, 3.3),
+      c(0, ymax),
+      type = "n",
+      xaxt = "n",
+      xlab = "",
+      ylab = "Standard deviation",
+      main = ""
+    )
+    axis(1, at = 1:3, labels = c("Random intercept", "Random slope", "Residual"))
+
+    for (i in seq_along(true_sd)) {
+      if (!is.na(estimated_sd[i])) {
+        segments(i, true_sd[i], i, estimated_sd[i], col = "grey70")
+      }
+    }
+
+    points(1:3, true_sd, pch = 1, cex = 1.35, lwd = 2)
+    keep <- !is.na(estimated_sd)
+    points((1:3)[keep], estimated_sd[keep], pch = 16, cex = 1.15)
+
+    legend(
+      "topright",
+      legend = c("True", "Estimated"),
+      pch = c(1, 16),
+      bty = "n",
+      cex = 0.82
+    )
+  })
+}
+
+shinyApp(ui, server)

From 59f670a3b4e56e712a85ba9f913975cfa08275c8 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 17:53:32 -0400
Subject: [PATCH 05/19] Refine mixed models diagnostic layout

---
 mixed-models/app.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index 8ac2569..5a6aeb1 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -247,7 +247,7 @@ ui <- page_fillable(
       tags$small(htmltools::includeMarkdown("credits.md"))
     ),
     layout_columns(
-      col_widths = c(12, 6, 6, 6, 6),
+      col_widths = c(12, 6, 6, 4, 4, 4),
       gap = "0.75rem",
       card(
         card_header(uiOutput("main_title")),

From f7114b494a9e334edada6c03b84ded87d86023c9 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 22:37:39 -0400
Subject: [PATCH 06/19] Fix mixed model selector values

---
 mixed-models/app.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index 5a6aeb1..20be465 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -201,7 +201,7 @@ ui <- page_fillable(
           "1. Data structure",
           "Changes the process that generates the points."
         ),
-        choices = structure_labels,
+        choices = setNames(names(structure_labels), structure_labels),
         selected = "both"
       ),
       radioButtons(
@@ -210,7 +210,7 @@ ui <- page_fillable(
           "2. Fitted model",
           "Keeps the data fixed and changes how group structure is modeled."
         ),
-        choices = model_labels,
+        choices = setNames(names(model_labels), model_labels),
         selected = "ri"
       ),
       sliderInput(
@@ -628,4 +628,4 @@ server <- function(input, output, session) {
   })
 }
 
-shinyApp(ui, server)
+shinyApp(ui, server)
\ No newline at end of file

From 34c21d624bdf1039bbc0947a83adebf840478ed4 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 22:42:51 -0400
Subject: [PATCH 07/19] Fix MathJax escaping in mixed models readme

---
 mixed-models/readme.md | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/mixed-models/readme.md b/mixed-models/readme.md
index e2a26a6..c44dbaf 100644
--- a/mixed-models/readme.md
+++ b/mixed-models/readme.md
@@ -7,17 +7,17 @@ This experiment separates two ideas that are easy to mix up:
 
 The simulated relationship is
 
-\[
-y_{ij} = \beta_0 + b_{0j} + (\beta_1 + b_{1j})x_{ij} + \varepsilon_{ij},
-\]
+\\[
+y_{ij} = \\beta_0 + b_{0j} + (\\beta_1 + b_{1j})x_{ij} + \\varepsilon_{ij},
+\\]
 
 where the group deviations are centered at zero and the observation noise is normal:
 
-\[
-\begin{pmatrix}b_{0j}\\b_{1j}\end{pmatrix} \sim N(0, \Sigma),
-\qquad
-\varepsilon_{ij} \sim N(0, \sigma^2).
-\]
+\\[
+\\begin{pmatrix}b_{0j}\\\\b_{1j}\\end{pmatrix} \\sim N(0, \\Sigma),
+\\qquad
+\\varepsilon_{ij} \\sim N(0, \\sigma^2).
+\\]
 
 Depending on the selected data structure, the random-intercept or random-slope variance can be exactly zero.
 

From 99b2474b9100688e649853cf571a42ea763df6d7 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 22:47:07 -0400
Subject: [PATCH 08/19] Refine mixed models explorer layout

---
 mixed-models/app.R | 103 ++++++++++++++++++++++++++++++++++++---------
 1 file changed, 83 insertions(+), 20 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index 20be465..d8762ce 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -33,7 +33,6 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   groups <- paste0("G", seq_len(n_groups))
   z0 <- rnorm(n_groups)
   z1 <- rnorm(n_groups)
-
   b0 <- sigma_b0 * z0
   b1 <- sigma_b1 * (rho * z0 + sqrt(1 - rho^2) * z1)
 
@@ -190,6 +189,60 @@ structure_labels <- c(
 ui <- page_fillable(
   theme = apptheme,
   padding = 0,
+  tags$head(
+    tags$style(HTML("
+      .formula-block {
+        margin: -0.15rem 0 0.75rem 0;
+        padding: 0.55rem 0.65rem;
+        border-left: 3px solid var(--bs-primary);
+        background: color-mix(in srgb, var(--bs-primary) 5%, transparent);
+        font-size: 0.82rem;
+      }
+
+      .formula-block code {
+        display: block;
+        margin-top: 0.25rem;
+        white-space: normal;
+      }
+
+      .mixed-grid {
+        display: grid;
+        grid-template-columns: minmax(0, 1fr) minmax(0, 1fr) minmax(230px, 0.55fr);
+        grid-template-rows: minmax(0, 1fr) minmax(0, 1fr) minmax(170px, 0.55fr);
+        gap: 0.75rem;
+        width: 100%;
+        height: calc(100vh - 1.5rem);
+        min-height: 760px;
+      }
+
+      .mixed-main { grid-column: 1 / span 2; grid-row: 1 / span 2; }
+      .mixed-residuals { grid-column: 3; grid-row: 1; }
+      .mixed-qq { grid-column: 3; grid-row: 2; }
+      .mixed-random { grid-column: 1; grid-row: 3; }
+      .mixed-shrinkage { grid-column: 2; grid-row: 3; }
+      .mixed-variance { grid-column: 3; grid-row: 3; }
+
+      .mixed-grid .card { min-width: 0; min-height: 0; }
+      .mixed-grid .shiny-plot-output { height: 100% !important; min-height: 0; }
+
+      @media (max-width: 1100px) {
+        .mixed-grid {
+          display: grid;
+          grid-template-columns: 1fr 1fr;
+          grid-template-rows: auto;
+          height: auto;
+          min-height: 0;
+        }
+
+        .mixed-main { grid-column: 1 / -1; grid-row: auto; min-height: 65vh; }
+        .mixed-residuals, .mixed-qq, .mixed-random, .mixed-shrinkage, .mixed-variance {
+          grid-column: auto;
+          grid-row: auto;
+          min-height: 260px;
+        }
+      }
+    "))
+  ),
   layout_sidebar(
     fillable = TRUE,
     padding = "0.75rem",
@@ -204,6 +257,11 @@ ui <- page_fillable(
         choices = setNames(names(structure_labels), structure_labels),
         selected = "both"
       ),
+      tags$div(
+        class = "formula-block",
+        tags$strong("Underlying model"),
+        uiOutput("truth_formula")
+      ),
       radioButtons(
         "model",
         input_label_vdl(
@@ -213,6 +271,12 @@ ui <- page_fillable(
         choices = setNames(names(model_labels), model_labels),
         selected = "ri"
       ),
+      tags$div(
+        class = "formula-block",
+        tags$strong("Fitted model"),
+        uiOutput("fit_formula"),
+        uiOutput("pooling_note")
+      ),
       sliderInput(
         "n_groups",
         tags$small("Groups"),
@@ -230,12 +294,6 @@ ui <- page_fillable(
         step = 1
       ),
       actionButton("resimulate", "Resimulate data", width = "100%"),
-      tags$hr(),
-      tags$small(tags$strong("Underlying model")),
-      uiOutput("truth_formula"),
-      tags$small(tags$strong("Fitted model")),
-      uiOutput("fit_formula"),
-      uiOutput("pooling_note"),
       uiOutput("model_check"),
       accordion(
         open = FALSE,
@@ -246,32 +304,37 @@ ui <- page_fillable(
       ),
       tags$small(htmltools::includeMarkdown("credits.md"))
     ),
-    layout_columns(
-      col_widths = c(12, 6, 6, 4, 4, 4),
-      gap = "0.75rem",
+    tags$div(
+      class = "mixed-grid",
       card(
+        class = "mixed-main",
         card_header(uiOutput("main_title")),
-        plotOutput("main_plot", width = "100%", height = "45vh")
+        plotOutput("main_plot", width = "100%", height = "100%")
       ),
       card(
+        class = "mixed-residuals",
         card_header("Residuals vs fitted"),
-        plotOutput("residual_plot", width = "100%", height = "27vh")
+        plotOutput("residual_plot", width = "100%", height = "100%")
       ),
       card(
+        class = "mixed-qq",
         card_header("Normal Q-Q"),
-        plotOutput("qq_plot", width = "100%", height = "27vh")
+        plotOutput("qq_plot", width = "100%", height = "100%")
       ),
       card(
+        class = "mixed-random",
         card_header("Random effects"),
-        plotOutput("random_effects_plot", width = "100%", height = "27vh")
+        plotOutput("random_effects_plot", width = "100%", height = "100%")
       ),
       card(
+        class = "mixed-shrinkage",
         card_header("Pooling / shrinkage"),
-        plotOutput("shrinkage_plot", width = "100%", height = "27vh")
+        plotOutput("shrinkage_plot", width = "100%", height = "100%")
       ),
       card(
+        class = "mixed-variance",
         card_header("Variance components"),
-        plotOutput("variance_plot", width = "100%", height = "27vh")
+        plotOutput("variance_plot", width = "100%", height = "100%")
       )
     )
   )
@@ -367,7 +430,7 @@ server <- function(input, output, session) {
       "Partial pooling · group deviations are estimated jointly and shrink toward zero."
     )
 
-    tags$div(class = "small text-muted mt-2", tags$strong("Information sharing: "), text)
+    tags$div(class = "small text-muted mt-1", text)
   })
 
   output$model_check <- renderUI({
@@ -377,7 +440,7 @@ server <- function(input, output, session) {
     if (!inherits(mod, "merMod")) {
       return(
         tags$div(
-          class = "small mt-2",
+          class = "small mt-2 mb-2",
           tags$strong("Model check"),
           tags$div("Random effects: not modeled"),
           tags$div(sprintf("Mean residual: %.3f", residual_mean))
@@ -391,7 +454,7 @@ server <- function(input, output, session) {
     sds <- estimated_sds(mod)
 
     tags$div(
-      class = "small mt-2",
+      class = "small mt-2 mb-2",
       tags$strong("Model check"),
       tags$div(if (converged) "✓ optimizer converged" else "⚠ convergence warning"),
       tags$div(
@@ -628,4 +691,4 @@ server <- function(input, output, session) {
   })
 }
 
-shinyApp(ui, server)
\ No newline at end of file
+shinyApp(ui, server)

From 03185a6926c5465d9dafaf8fd47df631d62712c3 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 22:53:44 -0400
Subject: [PATCH 09/19] Resize mixed models main panel to 60 percent

---
 mixed-models/app.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index d8762ce..dd851c9 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -207,8 +207,8 @@ ui <- page_fillable(
 
       .mixed-grid {
         display: grid;
-        grid-template-columns: minmax(0, 1fr) minmax(0, 1fr) minmax(230px, 0.55fr);
-        grid-template-rows: minmax(0, 1fr) minmax(0, 1fr) minmax(170px, 0.55fr);
+        grid-template-columns: minmax(0, 1fr) minmax(0, 1fr) minmax(230px, 1.33fr);
+        grid-template-rows: minmax(0, 1fr) minmax(0, 1fr) minmax(170px, 1.33fr);
         gap: 0.75rem;
         width: 100%;
         height: calc(100vh - 1.5rem);

From bc9f05d778dd214a932c9dff71be13ee1afe565f Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 23:10:32 -0400
Subject: [PATCH 10/19] Make mixed models dashboard grid square

---
 mixed-models/app.R | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index dd851c9..515ccd8 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -207,12 +207,12 @@ ui <- page_fillable(
 
       .mixed-grid {
         display: grid;
-        grid-template-columns: minmax(0, 1fr) minmax(0, 1fr) minmax(230px, 1.33fr);
-        grid-template-rows: minmax(0, 1fr) minmax(0, 1fr) minmax(170px, 1.33fr);
+        grid-template-columns: repeat(3, minmax(0, 1fr));
+        grid-template-rows: repeat(3, minmax(0, 1fr));
         gap: 0.75rem;
-        width: 100%;
-        height: calc(100vh - 1.5rem);
-        min-height: 760px;
+        width: min(100%, calc(100vh - 1.5rem));
+        aspect-ratio: 1 / 1;
+        margin-inline: auto;
       }
 
       .mixed-main { grid-column: 1 / span 2; grid-row: 1 / span 2; }
@@ -230,8 +230,9 @@ ui <- page_fillable(
           display: grid;
           grid-template-columns: 1fr 1fr;
           grid-template-rows: auto;
-          height: auto;
-          min-height: 0;
+          width: 100%;
+          aspect-ratio: auto;
+          margin-inline: 0;
         }
 
         .mixed-main { grid-column: 1 / -1; grid-row: auto; min-height: 65vh; }
@@ -691,4 +692,4 @@ server <- function(input, output, session) {
   })
 }
 
-shinyApp(ui, server)
+shinyApp(ui, server)
\ No newline at end of file

From 6e3ed0e7c1ca1c4c03076e98de711bc58ecce1e7 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Mon, 17 Aug 2026 23:17:37 -0400
Subject: [PATCH 11/19] Keep mixed models grid at two-thirds without square
 aspect

---
 mixed-models/app.R | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index 515ccd8..d2e5b49 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -210,9 +210,9 @@ ui <- page_fillable(
         grid-template-columns: repeat(3, minmax(0, 1fr));
         grid-template-rows: repeat(3, minmax(0, 1fr));
         gap: 0.75rem;
-        width: min(100%, calc(100vh - 1.5rem));
-        aspect-ratio: 1 / 1;
-        margin-inline: auto;
+        width: 100%;
+        height: calc(100vh - 1.5rem);
+        min-height: 720px;
       }
 
       .mixed-main { grid-column: 1 / span 2; grid-row: 1 / span 2; }
@@ -231,8 +231,8 @@ ui <- page_fillable(
           grid-template-columns: 1fr 1fr;
           grid-template-rows: auto;
           width: 100%;
-          aspect-ratio: auto;
-          margin-inline: 0;
+          height: auto;
+          min-height: 0;
         }
 
         .mixed-main { grid-column: 1 / -1; grid-row: auto; min-height: 65vh; }

From 99f06863d3a6ee19cf0ca65d56c66cc6ed8d9b30 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 02:03:54 -0400
Subject: [PATCH 12/19] Clarify grouped structure in mixed models simulation

---
 mixed-models/app.R | 35 +++++++++++++++++++++++++++--------
 1 file changed, 27 insertions(+), 8 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index d2e5b49..34407bc 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -36,6 +36,12 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   b0 <- sigma_b0 * z0
   b1 <- sigma_b1 * (rho * z0 + sqrt(1 - rho^2) * z1)
 
+  x_centers <- if (structure %in% c("intercept", "both")) {
+    seq(-1.4, 1.4, length.out = n_groups)
+  } else {
+    rep(0, n_groups)
+  }
+
   effects <- data.frame(
     group = factor(groups, levels = groups),
     b0 = b0,
@@ -45,7 +51,11 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   dat <- do.call(
     rbind,
     lapply(seq_len(n_groups), function(i) {
-      x <- sort(runif(n_per_group, -2, 2))
+      x <- if (structure %in% c("intercept", "both")) {
+        sort(x_centers[i] + runif(n_per_group, -0.85, 0.85))
+      } else {
+        sort(runif(n_per_group, -2, 2))
+      }
       epsilon <- rnorm(n_per_group, 0, sigma_e)
       y <- beta0 + b0[i] + (beta1 + b1[i]) * x + epsilon
 
@@ -84,12 +94,13 @@ fit_selected_model <- function(dat, model) {
 
 truth_lines <- function(sim) {
   groups <- levels(sim$data$group)
-  x_grid <- seq(-2.1, 2.1, length.out = 100)
 
   do.call(
     rbind,
     lapply(groups, function(g) {
       effect <- sim$effects[sim$effects$group == g, ]
+      observed <- sim$data[sim$data$group == g, ]
+      x_grid <- seq(min(observed$x), max(observed$x), length.out = 100)
 
       data.frame(
         group = factor(g, levels = groups),
@@ -103,11 +114,19 @@ truth_lines <- function(sim) {
 
 fitted_lines <- function(dat, fit) {
   groups <- levels(dat$group)
-  x_grid <- seq(-2.1, 2.1, length.out = 100)
-  grid <- expand.grid(x = x_grid, group = groups)
-  grid$group <- factor(grid$group, levels = groups)
-  grid$y <- predict(fit, newdata = grid)
-  grid
+
+  do.call(
+    rbind,
+    lapply(groups, function(g) {
+      observed <- dat[dat$group == g, ]
+      grid <- data.frame(
+        x = seq(min(observed$x), max(observed$x), length.out = 100),
+        group = factor(g, levels = groups)
+      )
+      grid$y <- predict(fit, newdata = grid)
+      grid
+    })
+  )
 }
 
 no_pool_coefficients <- function(dat) {
@@ -692,4 +711,4 @@ server <- function(input, output, session) {
   })
 }
 
-shinyApp(ui, server)
\ No newline at end of file
+shinyApp(ui, server)

From 3a6003dfac6a1616e5b337fcdba165751d2867ab Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 02:04:11 -0400
Subject: [PATCH 13/19] Explain staggered group ranges in mixed models app

---
 mixed-models/readme.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/mixed-models/readme.md b/mixed-models/readme.md
index c44dbaf..8865ab2 100644
--- a/mixed-models/readme.md
+++ b/mixed-models/readme.md
@@ -21,6 +21,8 @@ where the group deviations are centered at zero and the observation noise is nor
 
 Depending on the selected data structure, the random-intercept or random-slope variance can be exactly zero.
 
+For scenarios with intercept differences, groups are observed over partly different ranges of \\(x\\). Those ranges are generated independently of the random effects; they simply make the contrast between a pooled relationship and within-group relationships easier to see.
+
 ### Pooling
 
 **Complete pooling** ignores group differences and estimates one relationship for everyone.

From 4aafef64c36184141cf275c4e20d959d517ca7d7 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 02:16:44 -0400
Subject: [PATCH 14/19] Show group colors and coefficient recovery

---
 mixed-models/app.R | 201 ++++++++++++++++++++++++++++++---------------
 1 file changed, 137 insertions(+), 64 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index 34407bc..b7ca0b4 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -20,6 +20,12 @@ card <- purrr::partial(
 primary_color <- unname(bs_get_variables(apptheme, "primary"))
 
 # helpers -----------------------------------------------------------------
+group_colors <- function(groups) {
+  palette <- hcl.colors(256, "viridis")
+  index <- round(seq(1 + 0.12 * 255, 1 + 0.88 * 255, length.out = length(groups)))
+  setNames(palette[index], groups)
+}
+
 simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   set.seed(seed)
 
@@ -56,6 +62,7 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
       } else {
         sort(runif(n_per_group, -2, 2))
       }
+
       epsilon <- rnorm(n_per_group, 0, sigma_e)
       y <- beta0 + b0[i] + (beta1 + b1[i]) * x + epsilon
 
@@ -136,6 +143,15 @@ no_pool_coefficients <- function(dat) {
   out
 }
 
+truth_coefficients <- function(sim) {
+  out <- cbind(
+    intercept = sim$truth$beta0 + sim$effects$b0,
+    slope = sim$truth$beta1 + sim$effects$b1
+  )
+  rownames(out) <- as.character(sim$effects$group)
+  out
+}
+
 current_group_coefficients <- function(dat, fit, model) {
   groups <- levels(dat$group)
 
@@ -189,6 +205,34 @@ estimated_sds <- function(fit) {
   out
 }
 
+coefficient_recovery <- function(sim, fit, model) {
+  truth <- truth_coefficients(sim)
+  no_pool <- no_pool_coefficients(sim$data)
+  selected <- current_group_coefficients(sim$data, fit, model)
+  groups <- rownames(truth)
+
+  list(
+    table = data.frame(
+      Group = groups,
+      `True int.` = truth[groups, "intercept"],
+      `No pool int.` = no_pool[groups, "intercept"],
+      `Model int.` = selected[groups, "intercept"],
+      `True slope` = truth[groups, "slope"],
+      `No pool slope` = no_pool[groups, "slope"],
+      `Model slope` = selected[groups, "slope"],
+      check.names = FALSE
+    ),
+    no_pool_rmse = c(
+      intercept = sqrt(mean((no_pool[, "intercept"] - truth[, "intercept"])^2)),
+      slope = sqrt(mean((no_pool[, "slope"] - truth[, "slope"])^2))
+    ),
+    selected_rmse = c(
+      intercept = sqrt(mean((selected[, "intercept"] - truth[, "intercept"])^2)),
+      slope = sqrt(mean((selected[, "slope"] - truth[, "slope"])^2))
+    )
+  )
+}
+
 model_labels <- c(
   pooled = "Complete pooling",
   separate = "No pooling",
@@ -239,10 +283,13 @@ ui <- page_fillable(
       .mixed-qq { grid-column: 3; grid-row: 2; }
       .mixed-random { grid-column: 1; grid-row: 3; }
       .mixed-shrinkage { grid-column: 2; grid-row: 3; }
-      .mixed-variance { grid-column: 3; grid-row: 3; }
+      .mixed-recovery { grid-column: 3; grid-row: 3; }
 
       .mixed-grid .card { min-width: 0; min-height: 0; }
       .mixed-grid .shiny-plot-output { height: 100% !important; min-height: 0; }
+      .recovery-wrap { overflow: auto; padding: 0.25rem 0.4rem; font-size: 0.72rem; }
+      .recovery-wrap table { white-space: nowrap; margin-bottom: 0; }
+      .recovery-summary { padding: 0.25rem 0.4rem 0; font-size: 0.72rem; }
 
       @media (max-width: 1100px) {
         .mixed-grid {
@@ -255,7 +302,7 @@ ui <- page_fillable(
         }
 
         .mixed-main { grid-column: 1 / -1; grid-row: auto; min-height: 65vh; }
-        .mixed-residuals, .mixed-qq, .mixed-random, .mixed-shrinkage, .mixed-variance {
+        .mixed-residuals, .mixed-qq, .mixed-random, .mixed-shrinkage, .mixed-recovery {
           grid-column: auto;
           grid-row: auto;
           min-height: 260px;
@@ -352,9 +399,10 @@ ui <- page_fillable(
         plotOutput("shrinkage_plot", width = "100%", height = "100%")
       ),
       card(
-        class = "mixed-variance",
-        card_header("Variance components"),
-        plotOutput("variance_plot", width = "100%", height = "100%")
+        class = "mixed-recovery",
+        card_header("Coefficient recovery"),
+        tags$div(class = "recovery-summary", uiOutput("recovery_summary")),
+        tags$div(class = "recovery-wrap", tableOutput("recovery_table"))
       )
     )
   )
@@ -381,6 +429,10 @@ server <- function(input, output, session) {
     fit_selected_model(sim()$data, input$model)
   })
 
+  recovery <- reactive({
+    coefficient_recovery(sim(), fit(), input$model)
+  })
+
   output$truth_formula <- renderUI({
     truth <- sim()$truth
 
@@ -510,12 +562,12 @@ server <- function(input, output, session) {
     truth <- truth_lines(simulation)
     fitted <- fitted_lines(dat, mod)
     groups <- levels(dat$group)
-    cols <- setNames(hcl.colors(length(groups), "Dark 3"), groups)
+    cols <- group_colors(groups)
 
     plot(
       dat$x,
       dat$y,
-      col = adjustcolor(cols[as.character(dat$group)], alpha.f = 0.65),
+      col = adjustcolor(cols[as.character(dat$group)], alpha.f = 0.7),
       pch = 16,
       xlab = "x",
       ylab = "y",
@@ -524,7 +576,7 @@ server <- function(input, output, session) {
 
     for (g in groups) {
       dtruth <- truth[truth$group == g, ]
-      lines(dtruth$x, dtruth$y, col = adjustcolor(cols[g], alpha.f = 0.45), lty = 2, lwd = 2)
+      lines(dtruth$x, dtruth$y, col = adjustcolor(cols[g], alpha.f = 0.5), lty = 2, lwd = 2)
     }
 
     if (input$model == "pooled") {
@@ -547,7 +599,10 @@ server <- function(input, output, session) {
   })
 
   output$residual_plot <- renderPlot({
+    dat <- sim()$data
     mod <- fit()
+    groups <- levels(dat$group)
+    cols <- group_colors(groups)
     x <- fitted(mod)
     y <- residuals(mod)
 
@@ -555,7 +610,7 @@ server <- function(input, output, session) {
       x,
       y,
       pch = 16,
-      col = adjustcolor(primary_color, alpha.f = 0.55),
+      col = adjustcolor(cols[as.character(dat$group)], alpha.f = 0.72),
       xlab = "Fitted",
       ylab = "Residual",
       main = ""
@@ -565,8 +620,22 @@ server <- function(input, output, session) {
   })
 
   output$qq_plot <- renderPlot({
+    dat <- sim()$data
     r <- residuals(fit())
-    qqnorm(r, pch = 16, col = adjustcolor(primary_color, alpha.f = 0.55), main = "")
+    groups <- levels(dat$group)
+    cols <- group_colors(groups)
+    order_r <- order(r)
+    qq <- qqnorm(r, plot.it = FALSE)
+
+    plot(
+      qq$x,
+      qq$y,
+      pch = 16,
+      col = adjustcolor(cols[as.character(dat$group)[order_r]], alpha.f = 0.72),
+      xlab = "Theoretical Quantiles",
+      ylab = "Sample Quantiles",
+      main = ""
+    )
     qqline(r, lwd = 2)
   })
 
@@ -582,6 +651,7 @@ server <- function(input, output, session) {
 
     re <- ranef(mod)$group
     groups <- rownames(re)
+    cols <- group_colors(groups)
     values <- as.matrix(re)
     xr <- range(c(0, values))
     pad <- max(diff(xr) * 0.08, 0.1)
@@ -600,7 +670,14 @@ server <- function(input, output, session) {
 
     pchs <- c(16, 1)
     for (j in seq_len(ncol(values))) {
-      points(values[, j], seq_along(groups), pch = pchs[j], cex = 1.15)
+      points(
+        values[, j],
+        seq_along(groups),
+        pch = pchs[j],
+        col = cols[groups],
+        cex = 1.15,
+        lwd = 1.5
+      )
     }
 
     legend(
@@ -617,6 +694,7 @@ server <- function(input, output, session) {
     dat <- sim()$data
     mod <- fit()
     groups <- levels(dat$group)
+    cols <- group_colors(groups)
     no_pool <- no_pool_coefficients(dat)
     current <- current_group_coefficients(dat, mod, input$model)
     population <- population_coefficients(dat, mod, input$model)
@@ -640,21 +718,34 @@ server <- function(input, output, session) {
       no_pool[, "slope"],
       current[, "intercept"],
       current[, "slope"],
-      col = "grey70"
+      col = adjustcolor(cols[groups], alpha.f = 0.55),
+      lwd = 1.5
+    )
+    points(
+      no_pool[, "intercept"],
+      no_pool[, "slope"],
+      pch = 1,
+      col = cols[groups],
+      cex = 1.15,
+      lwd = 1.5
+    )
+    points(
+      current[, "intercept"],
+      current[, "slope"],
+      pch = 16,
+      col = cols[groups],
+      cex = 1.15
     )
-    points(no_pool[, "intercept"], no_pool[, "slope"], pch = 1, cex = 1.1)
-    points(current[, "intercept"], current[, "slope"], pch = 16, cex = 1.1)
     points(population[1], population[2], pch = 8, cex = 1.5, lwd = 2)
 
-    if (input$model %in% c("ri", "rs", "ris")) {
-      text(
-        current[, "intercept"],
-        current[, "slope"],
-        labels = groups,
-        pos = 3,
-        cex = 0.7
-      )
-    }
+    text(
+      current[, "intercept"],
+      current[, "slope"],
+      labels = groups,
+      col = cols[groups],
+      pos = 3,
+      cex = 0.7
+    )
 
     legend(
       "topright",
@@ -665,50 +756,32 @@ server <- function(input, output, session) {
     )
   })
 
-  output$variance_plot <- renderPlot({
-    truth <- sim()$truth
-    mod <- fit()
-
-    true_sd <- c(
-      intercept = truth$sigma_b0,
-      slope = truth$sigma_b1,
-      residual = truth$sigma_e
-    )
-    estimated_sd <- estimated_sds(mod)
-
-    values <- c(true_sd, estimated_sd)
-    ymax <- max(values, na.rm = TRUE) * 1.2
-    if (!is.finite(ymax) || ymax == 0) ymax <- 1
+  output$recovery_summary <- renderUI({
+    x <- recovery()
 
-    plot(
-      c(0.7, 3.3),
-      c(0, ymax),
-      type = "n",
-      xaxt = "n",
-      xlab = "",
-      ylab = "Standard deviation",
-      main = ""
-    )
-    axis(1, at = 1:3, labels = c("Random intercept", "Random slope", "Residual"))
-
-    for (i in seq_along(true_sd)) {
-      if (!is.na(estimated_sd[i])) {
-        segments(i, true_sd[i], i, estimated_sd[i], col = "grey70")
-      }
-    }
-
-    points(1:3, true_sd, pch = 1, cex = 1.35, lwd = 2)
-    keep <- !is.na(estimated_sd)
-    points((1:3)[keep], estimated_sd[keep], pch = 16, cex = 1.15)
-
-    legend(
-      "topright",
-      legend = c("True", "Estimated"),
-      pch = c(1, 16),
-      bty = "n",
-      cex = 0.82
+    tags$div(
+      tags$div(
+        tags$strong("No pooling RMSE: "),
+        sprintf("intercept %.2f · slope %.2f", x$no_pool_rmse["intercept"], x$no_pool_rmse["slope"])
+      ),
+      tags$div(
+        tags$strong(paste0(model_labels[[input$model]], " RMSE: ")),
+        sprintf("intercept %.2f · slope %.2f", x$selected_rmse["intercept"], x$selected_rmse["slope"])
+      )
     )
   })
+
+  output$recovery_table <- renderTable(
+    {
+      recovery()$table
+    },
+    digits = 2,
+    striped = TRUE,
+    bordered = FALSE,
+    hover = TRUE,
+    spacing = "xs",
+    rownames = FALSE
+  )
 }
 
 shinyApp(ui, server)

From 80a0d4f5a193131700736cd17f94e274404d83f1 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 02:17:29 -0400
Subject: [PATCH 15/19] Document coefficient recovery view

---
 mixed-models/readme.md | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/mixed-models/readme.md b/mixed-models/readme.md
index 8865ab2..39895af 100644
--- a/mixed-models/readme.md
+++ b/mixed-models/readme.md
@@ -39,11 +39,13 @@ Groups with less information generally shrink more. Groups with more observation
 
 ### What to inspect
 
+Group colors are kept consistent across the main view, residual diagnostics, random effects, and shrinkage plot so the same group can be followed through the model.
+
 - **Main view:** simulated truth and fitted group relationships.
-- **Residuals vs fitted:** remaining structure or changing residual spread.
-- **Normal Q-Q:** whether residuals look compatible with a normal-error assumption.
+- **Residuals vs fitted:** remaining structure or changing residual spread, colored by group.
+- **Normal Q-Q:** whether residuals look compatible with a normal-error assumption and whether departures concentrate in particular groups.
 - **Random effects:** estimated group deviations around zero.
 - **Shrinkage:** no-pooling estimates compared with the estimates from the selected model.
-- **Variance components:** true simulation standard deviations compared with those estimated by the fitted mixed model.
+- **Coefficient recovery:** true simulated intercepts and slopes compared with no-pooling and selected-model estimates. RMSE summarizes how closely each approach recovers the truth.
 
 A singular mixed-model fit is informative here: it often means the fitted random-effect covariance has reached a boundary, commonly because one random-effect variance is estimated close to zero.

From c8cd734b7244aba77b74f9d6b85ab890238e0cef Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 02:32:46 -0400
Subject: [PATCH 16/19] Add no-model view and train-test comparison

---
 mixed-models/app.R | 202 +++++++++++++++++++++++++++++----------------
 1 file changed, 133 insertions(+), 69 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index b7ca0b4..bce880e 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -26,6 +26,11 @@ group_colors <- function(groups) {
   setNames(palette[index], groups)
 }
 
+empty_plot <- function(text) {
+  plot.new()
+  text(0.5, 0.5, text, cex = 0.95, col = "grey40")
+}
+
 simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   set.seed(seed)
 
@@ -91,6 +96,7 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
 fit_selected_model <- function(dat, model) {
   switch(
     model,
+    none = NULL,
     pooled = lm(y ~ x, data = dat),
     separate = lm(y ~ x * group, data = dat),
     ri = lmer(y ~ x + (1 | group), data = dat, REML = TRUE),
@@ -99,6 +105,51 @@ fit_selected_model <- function(dat, model) {
   )
 }
 
+predict_selected_model <- function(fit, newdata) {
+  if (inherits(fit, "merMod")) {
+    predict(fit, newdata = newdata, allow.new.levels = TRUE)
+  } else {
+    predict(fit, newdata = newdata)
+  }
+}
+
+split_grouped_data <- function(dat, seed, train_share = 0.7) {
+  set.seed(seed + 10000L)
+
+  train_rows <- unlist(
+    lapply(split(seq_len(nrow(dat)), dat$group), function(index) {
+      n_train <- floor(length(index) * train_share)
+      n_train <- max(2L, min(length(index) - 1L, n_train))
+      sample(index, n_train)
+    }),
+    use.names = FALSE
+  )
+
+  list(
+    train = dat[train_rows, , drop = FALSE],
+    test = dat[-train_rows, , drop = FALSE]
+  )
+}
+
+compare_models <- function(dat, seed) {
+  split <- split_grouped_data(dat, seed)
+  models <- c("pooled", "separate", "ri", "rs", "ris")
+
+  rows <- lapply(models, function(model) {
+    mod <- suppressWarnings(fit_selected_model(split$train, model))
+    pred_train <- predict_selected_model(mod, split$train)
+    pred_test <- predict_selected_model(mod, split$test)
+
+    data.frame(
+      model = model,
+      train_rmse = sqrt(mean((split$train$y - pred_train)^2)),
+      test_rmse = sqrt(mean((split$test$y - pred_test)^2))
+    )
+  })
+
+  do.call(rbind, rows)
+}
+
 truth_lines <- function(sim) {
   groups <- levels(sim$data$group)
 
@@ -130,7 +181,7 @@ fitted_lines <- function(dat, fit) {
         x = seq(min(observed$x), max(observed$x), length.out = 100),
         group = factor(g, levels = groups)
       )
-      grid$y <- predict(fit, newdata = grid)
+      grid$y <- predict_selected_model(fit, grid)
       grid
     })
   )
@@ -143,15 +194,6 @@ no_pool_coefficients <- function(dat) {
   out
 }
 
-truth_coefficients <- function(sim) {
-  out <- cbind(
-    intercept = sim$truth$beta0 + sim$effects$b0,
-    slope = sim$truth$beta1 + sim$effects$b1
-  )
-  rownames(out) <- as.character(sim$effects$group)
-  out
-}
-
 current_group_coefficients <- function(dat, fit, model) {
   groups <- levels(dat$group)
 
@@ -205,35 +247,8 @@ estimated_sds <- function(fit) {
   out
 }
 
-coefficient_recovery <- function(sim, fit, model) {
-  truth <- truth_coefficients(sim)
-  no_pool <- no_pool_coefficients(sim$data)
-  selected <- current_group_coefficients(sim$data, fit, model)
-  groups <- rownames(truth)
-
-  list(
-    table = data.frame(
-      Group = groups,
-      `True int.` = truth[groups, "intercept"],
-      `No pool int.` = no_pool[groups, "intercept"],
-      `Model int.` = selected[groups, "intercept"],
-      `True slope` = truth[groups, "slope"],
-      `No pool slope` = no_pool[groups, "slope"],
-      `Model slope` = selected[groups, "slope"],
-      check.names = FALSE
-    ),
-    no_pool_rmse = c(
-      intercept = sqrt(mean((no_pool[, "intercept"] - truth[, "intercept"])^2)),
-      slope = sqrt(mean((no_pool[, "slope"] - truth[, "slope"])^2))
-    ),
-    selected_rmse = c(
-      intercept = sqrt(mean((selected[, "intercept"] - truth[, "intercept"])^2)),
-      slope = sqrt(mean((selected[, "slope"] - truth[, "slope"])^2))
-    )
-  )
-}
-
 model_labels <- c(
+  none = "No model",
   pooled = "Complete pooling",
   separate = "No pooling",
   ri = "Random intercept",
@@ -283,13 +298,13 @@ ui <- page_fillable(
       .mixed-qq { grid-column: 3; grid-row: 2; }
       .mixed-random { grid-column: 1; grid-row: 3; }
       .mixed-shrinkage { grid-column: 2; grid-row: 3; }
-      .mixed-recovery { grid-column: 3; grid-row: 3; }
+      .mixed-comparison { grid-column: 3; grid-row: 3; }
 
       .mixed-grid .card { min-width: 0; min-height: 0; }
       .mixed-grid .shiny-plot-output { height: 100% !important; min-height: 0; }
-      .recovery-wrap { overflow: auto; padding: 0.25rem 0.4rem; font-size: 0.72rem; }
-      .recovery-wrap table { white-space: nowrap; margin-bottom: 0; }
-      .recovery-summary { padding: 0.25rem 0.4rem 0; font-size: 0.72rem; }
+      .comparison-wrap { overflow: auto; padding: 0.25rem 0.4rem; font-size: 0.74rem; }
+      .comparison-wrap table { white-space: nowrap; margin-bottom: 0; }
+      .comparison-summary { padding: 0.3rem 0.5rem 0; font-size: 0.74rem; }
 
       @media (max-width: 1100px) {
         .mixed-grid {
@@ -302,7 +317,7 @@ ui <- page_fillable(
         }
 
         .mixed-main { grid-column: 1 / -1; grid-row: auto; min-height: 65vh; }
-        .mixed-residuals, .mixed-qq, .mixed-random, .mixed-shrinkage, .mixed-recovery {
+        .mixed-residuals, .mixed-qq, .mixed-random, .mixed-shrinkage, .mixed-comparison {
           grid-column: auto;
           grid-row: auto;
           min-height: 260px;
@@ -336,7 +351,7 @@ ui <- page_fillable(
           "Keeps the data fixed and changes how group structure is modeled."
         ),
         choices = setNames(names(model_labels), model_labels),
-        selected = "ri"
+        selected = "none"
       ),
       tags$div(
         class = "formula-block",
@@ -399,10 +414,10 @@ ui <- page_fillable(
         plotOutput("shrinkage_plot", width = "100%", height = "100%")
       ),
       card(
-        class = "mixed-recovery",
-        card_header("Coefficient recovery"),
-        tags$div(class = "recovery-summary", uiOutput("recovery_summary")),
-        tags$div(class = "recovery-wrap", tableOutput("recovery_table"))
+        class = "mixed-comparison",
+        card_header("Train / test RMSE"),
+        tags$div(class = "comparison-summary", uiOutput("comparison_summary")),
+        tags$div(class = "comparison-wrap", tableOutput("comparison_table"))
       )
     )
   )
@@ -429,8 +444,9 @@ server <- function(input, output, session) {
     fit_selected_model(sim()$data, input$model)
   })
 
-  recovery <- reactive({
-    coefficient_recovery(sim(), fit(), input$model)
+  comparison <- reactive({
+    req(input$model != "none")
+    compare_models(sim()$data, seed())
   })
 
   output$truth_formula <- renderUI({
@@ -470,6 +486,10 @@ server <- function(input, output, session) {
   })
 
   output$fit_formula <- renderUI({
+    if (input$model == "none") {
+      return(tags$div(class = "text-muted", "No fitted model"))
+    }
+
     math <- switch(
       input$model,
       pooled = "y_{ij} = \\beta_0 + \\beta_1x_{ij} + \\varepsilon_{ij}",
@@ -497,6 +517,7 @@ server <- function(input, output, session) {
   output$pooling_note <- renderUI({
     text <- switch(
       input$model,
+      none = "Explore the grouped data before fitting a model.",
       pooled = "Complete pooling · one relationship is shared by every group.",
       separate = "No pooling · each group gets its own OLS relationship.",
       "Partial pooling · group deviations are estimated jointly and shrink toward zero."
@@ -506,6 +527,15 @@ server <- function(input, output, session) {
   })
 
   output$model_check <- renderUI({
+    if (input$model == "none") {
+      return(
+        tags$div(
+          class = "small mt-2 mb-2 text-muted",
+          "No model selected"
+        )
+      )
+    }
+
     mod <- fit()
     residual_mean <- mean(residuals(mod))
 
@@ -558,9 +588,6 @@ server <- function(input, output, session) {
   output$main_plot <- renderPlot({
     simulation <- sim()
     dat <- simulation$data
-    mod <- fit()
-    truth <- truth_lines(simulation)
-    fitted <- fitted_lines(dat, mod)
     groups <- levels(dat$group)
     cols <- group_colors(groups)
 
@@ -574,6 +601,14 @@ server <- function(input, output, session) {
       main = ""
     )
 
+    if (input$model == "none") {
+      return(invisible())
+    }
+
+    mod <- fit()
+    truth <- truth_lines(simulation)
+    fitted <- fitted_lines(dat, mod)
+
     for (g in groups) {
       dtruth <- truth[truth$group == g, ]
       lines(dtruth$x, dtruth$y, col = adjustcolor(cols[g], alpha.f = 0.5), lty = 2, lwd = 2)
@@ -599,6 +634,11 @@ server <- function(input, output, session) {
   })
 
   output$residual_plot <- renderPlot({
+    if (input$model == "none") {
+      empty_plot("Select a fitted model")
+      return(invisible())
+    }
+
     dat <- sim()$data
     mod <- fit()
     groups <- levels(dat$group)
@@ -620,6 +660,11 @@ server <- function(input, output, session) {
   })
 
   output$qq_plot <- renderPlot({
+    if (input$model == "none") {
+      empty_plot("Select a fitted model")
+      return(invisible())
+    }
+
     dat <- sim()$data
     r <- residuals(fit())
     groups <- levels(dat$group)
@@ -640,12 +685,15 @@ server <- function(input, output, session) {
   })
 
   output$random_effects_plot <- renderPlot({
+    if (input$model == "none") {
+      empty_plot("Select a fitted model")
+      return(invisible())
+    }
+
     mod <- fit()
 
     if (!inherits(mod, "merMod")) {
-      plot.new()
-      text(0.5, 0.55, "No random effects in this model", cex = 1.05)
-      text(0.5, 0.43, "Complete/no pooling do not estimate b_j", cex = 0.85)
+      empty_plot("This model has no random effects")
       return(invisible())
     }
 
@@ -691,6 +739,11 @@ server <- function(input, output, session) {
   })
 
   output$shrinkage_plot <- renderPlot({
+    if (input$model == "none") {
+      empty_plot("Select a fitted model")
+      return(invisible())
+    }
+
     dat <- sim()$data
     mod <- fit()
     groups <- levels(dat$group)
@@ -756,24 +809,35 @@ server <- function(input, output, session) {
     )
   })
 
-  output$recovery_summary <- renderUI({
-    x <- recovery()
+  output$comparison_summary <- renderUI({
+    if (input$model == "none") {
+      return(tags$span(class = "text-muted", "Select a fitted model to compare generalization."))
+    }
 
-    tags$div(
-      tags$div(
-        tags$strong("No pooling RMSE: "),
-        sprintf("intercept %.2f · slope %.2f", x$no_pool_rmse["intercept"], x$no_pool_rmse["slope"])
-      ),
-      tags$div(
-        tags$strong(paste0(model_labels[[input$model]], " RMSE: ")),
-        sprintf("intercept %.2f · slope %.2f", x$selected_rmse["intercept"], x$selected_rmse["slope"])
-      )
+    tags$span(
+      "70% train / 30% test within each group. Test contains unseen observations from the same groups."
     )
   })
 
-  output$recovery_table <- renderTable(
+  output$comparison_table <- renderTable(
     {
-      recovery()$table
+      if (input$model == "none") {
+        return(NULL)
+      }
+
+      x <- comparison()
+      best <- which.min(x$test_rmse)
+
+      data.frame(
+        Model = paste0(
+          unname(model_labels[x$model]),
+          ifelse(x$model == input$model, "  ← selected", "")
+        ),
+        Train = x$train_rmse,
+        Test = x$test_rmse,
+        Best = ifelse(seq_len(nrow(x)) == best, "✓", ""),
+        check.names = FALSE
+      )
     },
     digits = 2,
     striped = TRUE,

From 4c455dbe6733296e2d85091d9803d1e14ffa41a5 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 02:33:03 -0400
Subject: [PATCH 17/19] Explain no-model view and train-test RMSE

---
 mixed-models/readme.md | 22 +++++++++++++++-------
 1 file changed, 15 insertions(+), 7 deletions(-)

diff --git a/mixed-models/readme.md b/mixed-models/readme.md
index 39895af..951db13 100644
--- a/mixed-models/readme.md
+++ b/mixed-models/readme.md
@@ -5,6 +5,8 @@ This experiment separates two ideas that are easy to mix up:
 1. **Data structure** controls how the simulated data are generated.
 2. **Fitted model** controls how we choose to explain those same data.
 
+Start with **No model** to inspect the grouped observations without fitted or true relationship lines. Then add models and compare what each one assumes.
+
 The simulated relationship is
 
 \\[
@@ -27,7 +29,7 @@ For scenarios with intercept differences, groups are observed over partly differ
 
 **Complete pooling** ignores group differences and estimates one relationship for everyone.
 
-**No pooling** estimates a separate relationship for each group.
+**No pooling** estimates a separate relationship for each group using only that group's observations.
 
 **Partial pooling** is the mixed-model middle ground: groups have their own effects, but those effects are estimated together and share information.
 
@@ -37,15 +39,21 @@ Shrinkage is the visible result of partial pooling. A group estimate is pulled t
 
 Groups with less information generally shrink more. Groups with more observations, less noise, or stronger evidence of genuine between-group differences generally shrink less.
 
-### What to inspect
+### Train and test RMSE
+
+The same simulated dataset is split within every group into approximately 70% training observations and 30% test observations. Each candidate model is fitted on the training observations and evaluated on both sets.
 
-Group colors are kept consistent across the main view, residual diagnostics, random effects, and shrinkage plot so the same group can be followed through the model.
+A low training RMSE only says that a model describes the observations it already saw well. The test RMSE asks whether that fitted relationship also predicts unseen observations from the **same groups**. This helps show why no pooling can fit small groups very closely yet produce less stable group-specific relationships than a partially pooled model.
+
+This is not yet a test on completely new groups. Predicting a group that was absent during fitting is a separate hierarchical-model question.
+
+### What to inspect
 
-- **Main view:** simulated truth and fitted group relationships.
+- **Main view:** raw grouped observations, then simulated truth and fitted group relationships after selecting a model.
 - **Residuals vs fitted:** remaining structure or changing residual spread, colored by group.
-- **Normal Q-Q:** whether residuals look compatible with a normal-error assumption and whether departures concentrate in particular groups.
+- **Normal Q-Q:** whether residuals look compatible with a normal-error assumption, with group colors retained.
 - **Random effects:** estimated group deviations around zero.
-- **Shrinkage:** no-pooling estimates compared with the estimates from the selected model.
-- **Coefficient recovery:** true simulated intercepts and slopes compared with no-pooling and selected-model estimates. RMSE summarizes how closely each approach recovers the truth.
+- **Shrinkage:** no-pooling estimates compared with estimates from the selected model.
+- **Train / test RMSE:** predictive performance of all candidate models on the same split.
 
 A singular mixed-model fit is informative here: it often means the fitted random-effect covariance has reached a boundary, commonly because one random-effect variance is estimated close to zero.

From efc5784f6a5c8bbbf5be37e532e80f14e943001c Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 13:30:49 -0400
Subject: [PATCH 18/19] Add balanced and unbalanced group size scenarios

---
 mixed-models/app.R | 125 ++++++++++++++++++++++++++++++++++++---------
 1 file changed, 101 insertions(+), 24 deletions(-)

diff --git a/mixed-models/app.R b/mixed-models/app.R
index bce880e..ecc801d 100644
--- a/mixed-models/app.R
+++ b/mixed-models/app.R
@@ -31,7 +31,21 @@ empty_plot <- function(text) {
   text(0.5, 0.5, text, cex = 0.95, col = "grey40")
 }
 
-simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
+make_group_sizes <- function(pattern, n_groups, n_per_group) {
+  if (pattern == "balanced") {
+    return(rep(as.integer(n_per_group), n_groups))
+  }
+
+  as.integer(round(exp(seq(log(4), log(60), length.out = n_groups))))
+}
+
+simulate_grouped_data <- function(
+  structure,
+  n_groups,
+  n_per_group,
+  group_size_pattern,
+  seed
+) {
   set.seed(seed)
 
   beta0 <- 2
@@ -42,6 +56,7 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   rho <- if (structure == "both") 0.35 else 0
 
   groups <- paste0("G", seq_len(n_groups))
+  sizes <- make_group_sizes(group_size_pattern, n_groups, n_per_group)
   z0 <- rnorm(n_groups)
   z1 <- rnorm(n_groups)
   b0 <- sigma_b0 * z0
@@ -62,13 +77,15 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   dat <- do.call(
     rbind,
     lapply(seq_len(n_groups), function(i) {
+      n_i <- sizes[i]
+
       x <- if (structure %in% c("intercept", "both")) {
-        sort(x_centers[i] + runif(n_per_group, -0.85, 0.85))
+        sort(x_centers[i] + runif(n_i, -0.85, 0.85))
       } else {
-        sort(runif(n_per_group, -2, 2))
+        sort(runif(n_i, -2, 2))
       }
 
-      epsilon <- rnorm(n_per_group, 0, sigma_e)
+      epsilon <- rnorm(n_i, 0, sigma_e)
       y <- beta0 + b0[i] + (beta1 + b1[i]) * x + epsilon
 
       data.frame(
@@ -82,6 +99,7 @@ simulate_grouped_data <- function(structure, n_groups, n_per_group, seed) {
   list(
     data = dat,
     effects = effects,
+    group_sizes = setNames(sizes, groups),
     truth = list(
       beta0 = beta0,
       beta1 = beta1,
@@ -136,14 +154,40 @@ compare_models <- function(dat, seed) {
   models <- c("pooled", "separate", "ri", "rs", "ris")
 
   rows <- lapply(models, function(model) {
-    mod <- suppressWarnings(fit_selected_model(split$train, model))
-    pred_train <- predict_selected_model(mod, split$train)
-    pred_test <- predict_selected_model(mod, split$test)
+    mod <- tryCatch(
+      suppressWarnings(fit_selected_model(split$train, model)),
+      error = function(e) NULL
+    )
+
+    if (is.null(mod)) {
+      return(
+        data.frame(
+          model = model,
+          train_rmse = NA_real_,
+          test_rmse = NA_real_,
+          equal_group_test_rmse = NA_real_
+        )
+      )
+    }
+
+    pred_train <- tryCatch(
+      predict_selected_model(mod, split$train),
+      error = function(e) rep(NA_real_, nrow(split$train))
+    )
+    pred_test <- tryCatch(
+      predict_selected_model(mod, split$test),
+      error = function(e) rep(NA_real_, nrow(split$test))
+    )
+
+    train_error2 <- (split$train$y - pred_train)^2
+    test_error2 <- (split$test$y - pred_test)^2
+    group_mse <- tapply(test_error2, split$test$group, mean, na.rm = TRUE)
 
     data.frame(
       model = model,
-      train_rmse = sqrt(mean((split$train$y - pred_train)^2)),
-      test_rmse = sqrt(mean((split$test$y - pred_test)^2))
+      train_rmse = sqrt(mean(train_error2, na.rm = TRUE)),
+      test_rmse = sqrt(mean(test_error2, na.rm = TRUE)),
+      equal_group_test_rmse = sqrt(mean(group_mse, na.rm = TRUE))
     )
   })
 
@@ -249,8 +293,8 @@ estimated_sds <- function(fit) {
 
 model_labels <- c(
   none = "No model",
-  pooled = "Complete pooling",
-  separate = "No pooling",
+  pooled = "Global model",
+  separate = "Group-specific models",
   ri = "Random intercept",
   rs = "Random slope",
   ris = "Random intercept + slope"
@@ -302,9 +346,9 @@ ui <- page_fillable(
 
       .mixed-grid .card { min-width: 0; min-height: 0; }
       .mixed-grid .shiny-plot-output { height: 100% !important; min-height: 0; }
-      .comparison-wrap { overflow: auto; padding: 0.25rem 0.4rem; font-size: 0.74rem; }
+      .comparison-wrap { overflow: auto; padding: 0.25rem 0.4rem; font-size: 0.72rem; }
       .comparison-wrap table { white-space: nowrap; margin-bottom: 0; }
-      .comparison-summary { padding: 0.3rem 0.5rem 0; font-size: 0.74rem; }
+      .comparison-summary { padding: 0.3rem 0.5rem 0; font-size: 0.72rem; }
 
       @media (max-width: 1100px) {
         .mixed-grid {
@@ -367,14 +411,27 @@ ui <- page_fillable(
         value = 6,
         step = 1
       ),
-      sliderInput(
-        "n_per_group",
-        tags$small("Observations per group"),
-        min = 5,
-        max = 40,
-        value = 18,
-        step = 1
+      radioButtons(
+        "group_size_pattern",
+        input_label_vdl(
+          "Group sizes",
+          "Balanced gives every group the same amount of data. Unbalanced mixes small and large groups."
+        ),
+        choices = c("Balanced" = "balanced", "Unbalanced" = "unbalanced"),
+        selected = "unbalanced"
+      ),
+      conditionalPanel(
+        condition = "input.group_size_pattern == 'balanced'",
+        sliderInput(
+          "n_per_group",
+          tags$small("Observations per group"),
+          min = 5,
+          max = 40,
+          value = 18,
+          step = 1
+        )
       ),
+      uiOutput("group_size_note"),
       actionButton("resimulate", "Resimulate data", width = "100%"),
       uiOutput("model_check"),
       accordion(
@@ -436,6 +493,7 @@ server <- function(input, output, session) {
       structure = input$data_structure,
       n_groups = input$n_groups,
       n_per_group = input$n_per_group,
+      group_size_pattern = input$group_size_pattern,
       seed = seed()
     )
   })
@@ -449,6 +507,23 @@ server <- function(input, output, session) {
     compare_models(sim()$data, seed())
   })
 
+  output$group_size_note <- renderUI({
+    sizes <- make_group_sizes(
+      input$group_size_pattern,
+      input$n_groups,
+      input$n_per_group
+    )
+    groups <- paste0("G", seq_len(input$n_groups))
+
+    tags$div(
+      class = "small text-muted mb-2",
+      paste0(
+        "n by group: ",
+        paste(paste0(groups, "=", sizes), collapse = " · ")
+      )
+    )
+  })
+
   output$truth_formula <- renderUI({
     truth <- sim()$truth
 
@@ -518,8 +593,8 @@ server <- function(input, output, session) {
     text <- switch(
       input$model,
       none = "Explore the grouped data before fitting a model.",
-      pooled = "Complete pooling · one relationship is shared by every group.",
-      separate = "No pooling · each group gets its own OLS relationship.",
+      pooled = "Complete pooling · one global relationship is shared by every group.",
+      separate = "No pooling · each group gets an independent OLS relationship.",
       "Partial pooling · group deviations are estimated jointly and shrink toward zero."
     )
 
@@ -815,7 +890,7 @@ server <- function(input, output, session) {
     }
 
     tags$span(
-      "70% train / 30% test within each group. Test contains unseen observations from the same groups."
+      "70% train / 30% test within each group. Test weights rows; Equal-group gives every group the same weight."
     )
   })
 
@@ -826,7 +901,8 @@ server <- function(input, output, session) {
       }
 
       x <- comparison()
-      best <- which.min(x$test_rmse)
+      score <- replace(x$equal_group_test_rmse, is.na(x$equal_group_test_rmse), Inf)
+      best <- which.min(score)
 
       data.frame(
         Model = paste0(
@@ -835,6 +911,7 @@ server <- function(input, output, session) {
         ),
         Train = x$train_rmse,
         Test = x$test_rmse,
+        `Equal-group test` = x$equal_group_test_rmse,
         Best = ifelse(seq_len(nrow(x)) == best, "✓", ""),
         check.names = FALSE
       )

From 9a113e4a1a697c572c7bf842a4ea5de913f208a9 Mon Sep 17 00:00:00 2001
From: Joshua Kunst <jbkunst@gmail.com>
Date: Tue, 18 Aug 2026 13:31:20 -0400
Subject: [PATCH 19/19] Explain unbalanced groups and equal-group RMSE

---
 mixed-models/readme.md | 26 +++++++++++++++++++-------
 1 file changed, 19 insertions(+), 7 deletions(-)

diff --git a/mixed-models/readme.md b/mixed-models/readme.md
index 951db13..592f093 100644
--- a/mixed-models/readme.md
+++ b/mixed-models/readme.md
@@ -23,15 +23,23 @@ where the group deviations are centered at zero and the observation noise is nor
 
 Depending on the selected data structure, the random-intercept or random-slope variance can be exactly zero.
 
-For scenarios with intercept differences, groups are observed over partly different ranges of \\(x\\). Those ranges are generated independently of the random effects; they simply make the contrast between a pooled relationship and within-group relationships easier to see.
+For scenarios with intercept differences, groups are observed over partly different ranges of \\(x\\). Those ranges are generated independently of the random effects; they simply make the contrast between a global relationship and within-group relationships easier to see.
+
+### Group sizes
+
+**Balanced** gives every group the same number of observations.
+
+**Unbalanced** deliberately mixes very small and large groups. With six groups the sizes are approximately 4, 7, 12, 20, 35, and 60 observations. This is useful for seeing when estimating every group independently becomes unstable and when partial pooling can help the smaller groups borrow information from the population.
 
 ### Pooling
 
-**Complete pooling** ignores group differences and estimates one relationship for everyone.
+The selector uses descriptive model names while the text below it keeps the standard pooling terminology.
 
-**No pooling** estimates a separate relationship for each group using only that group's observations.
+**Global model** is complete pooling: it ignores group differences and estimates one relationship for everyone.
 
-**Partial pooling** is the mixed-model middle ground: groups have their own effects, but those effects are estimated together and share information.
+**Group-specific models** are no pooling: each group gets a separate relationship estimated only from that group's observations.
+
+**Random intercept**, **random slope**, and **random intercept + slope** are partial-pooling models: groups may differ, but their deviations are estimated jointly.
 
 ### Shrinkage
 
@@ -43,7 +51,11 @@ Groups with less information generally shrink more. Groups with more observation
 
 The same simulated dataset is split within every group into approximately 70% training observations and 30% test observations. Each candidate model is fitted on the training observations and evaluated on both sets.
 
-A low training RMSE only says that a model describes the observations it already saw well. The test RMSE asks whether that fitted relationship also predicts unseen observations from the **same groups**. This helps show why no pooling can fit small groups very closely yet produce less stable group-specific relationships than a partially pooled model.
+A low training RMSE only says that a model describes observations it already saw well. Test RMSE asks whether that fitted relationship predicts unseen observations from the **same groups**.
+
+With unbalanced groups, ordinary test RMSE gives more influence to large groups because they contribute more rows. The app therefore also reports **Equal-group test RMSE**: it computes test error within each group and then gives every group the same weight. This makes performance on small groups visible instead of letting the largest groups dominate the summary.
+
+Neither metric is expected to make a mixed model win every simulation. The point is to see **when** sharing information helps and when a group-specific model already has enough data to work well on its own.
 
 This is not yet a test on completely new groups. Predicting a group that was absent during fitting is a separate hierarchical-model question.
 
@@ -53,7 +65,7 @@ This is not yet a test on completely new groups. Predicting a group that was abs
 - **Residuals vs fitted:** remaining structure or changing residual spread, colored by group.
 - **Normal Q-Q:** whether residuals look compatible with a normal-error assumption, with group colors retained.
 - **Random effects:** estimated group deviations around zero.
-- **Shrinkage:** no-pooling estimates compared with estimates from the selected model.
-- **Train / test RMSE:** predictive performance of all candidate models on the same split.
+- **Shrinkage:** group-specific estimates compared with estimates from the selected model.
+- **Train / test RMSE:** predictive performance of all candidate models on the same split, including an equal-group test score.
 
 A singular mixed-model fit is informative here: it often means the fitted random-effect covariance has reached a boundary, commonly because one random-effect variance is estimated close to zero.