Figure gallery • drmTMB

This gallery shows the kinds of figures drmTMB users can build for model interpretation, uncertainty display, fitted interactions, fitted correlation summaries, and later simulation diagnostics. It is different from the internal count-pilot diagnostics used while developing the first Poisson/NB2 simulation plumbing.

The examples here use small simulated data so the page can build quickly. They are meant to show figure patterns, not biological conclusions. The interaction figures deliberately keep three pieces visible: the observed data, the fitted surface or cell means, and the conditioning values used to make the display.

library(drmTMB)
#> 
#> Attaching package: 'drmTMB'
#> The following object is masked from 'package:base':
#> 
#>     beta
library(ggplot2)

theme_drmtmb_gallery <- function() {
  theme_minimal(base_size = 11) +
    theme(
      panel.grid.minor = element_blank(),
      panel.grid.major.x = element_blank(),
      axis.title = element_text(colour = "grey15"),
      plot.title.position = "plot",
      plot.title = element_text(face = "bold"),
      plot.subtitle = element_text(colour = "grey30"),
      strip.text = element_text(face = "bold"),
      legend.position = "bottom",
      legend.title = element_text(colour = "grey20"),
      legend.key.width = grid::unit(1.1, "lines")
    )
}

gallery_colours <- c(
  blue = "#0072B2",
  green = "#009E73",
  orange = "#D55E00",
  yellow = "#E69F00",
  pink = "#CC79A7",
  sky = "#56B4E9",
  charcoal = "#4D4D4D"
)
gallery_colour <- function(name) unname(gallery_colours[name])

habitat_palette <- c(
  "reef" = gallery_colour("blue"),
  "kelp" = gallery_colour("orange")
)
season_palette <- c(
  "dry" = gallery_colour("blue"),
  "wet" = gallery_colour("yellow")
)
surface_palette <- c(
  "Gaussian location-scale" = gallery_colour("blue"),
  "Proportion" = gallery_colour("green"),
  "Count" = gallery_colour("pink")
)
make_confidence_distribution <- function(
  parameter,
  estimate,
  std_error,
  transform = c("identity", "exp"),
  null_value = 0,
  n = 500
) {
  transform <- match.arg(transform)
  theta <- seq(
    estimate - 4 * std_error,
    estimate + 4 * std_error,
    length.out = n
  )
  value <- switch(
    transform,
    identity = theta,
    exp = exp(theta)
  )
  density <- switch(
    transform,
    identity = dnorm(theta, mean = estimate, sd = std_error),
    exp = dnorm(theta, mean = estimate, sd = std_error) / value
  )
  data.frame(
    parameter = parameter,
    theta = theta,
    value = value,
    cdf = pnorm(theta, mean = estimate, sd = std_error),
    density = density,
    null_value = switch(transform, identity = null_value, exp = exp(null_value)),
    estimate_value = switch(transform, identity = estimate, exp = exp(estimate))
  )
}

make_confidence_intervals <- function(
  parameter,
  estimate,
  std_error,
  transform = c("identity", "exp"),
  levels = c(0.66, 0.95)
) {
  transform <- match.arg(transform)
  z_low <- qnorm((1 - levels) / 2)
  z_high <- qnorm(1 - (1 - levels) / 2)
  low <- estimate + z_low * std_error
  high <- estimate + z_high * std_error
  data.frame(
    parameter = parameter,
    level = sprintf("%d%%", round(100 * levels)),
    y = seq(0.10, 0.22, length.out = length(levels)),
    xmin = switch(transform, identity = low, exp = exp(low)),
    xmax = switch(transform, identity = high, exp = exp(high))
  )
}
scale_drmtmb_colour <- scale_colour_manual(
  values = habitat_palette,
  limits = names(habitat_palette)
)
scale_drmtmb_fill <- scale_fill_manual(
  values = habitat_palette,
  limits = names(habitat_palette)
)

Fitted mean with a confidence band

Start with the figure most readers expect: raw observations, fitted mean, and a 95% confidence band on the response scale. This example is also a categorical-by-continuous interaction: the temperature slope differs by habitat.

set.seed(2591)
n <- 160
habitat_levels <- c("reef", "kelp")
fish <- data.frame(
  temperature = runif(n, -1.6, 1.6),
  habitat = factor(
    sample(habitat_levels, n, replace = TRUE),
    levels = habitat_levels
  )
)
eta_mu <- 0.4 +
  0.7 * fish$temperature +
  0.35 * (fish$habitat == "kelp") +
  0.18 * fish$temperature * (fish$habitat == "kelp")
eta_sigma <- -0.55 + 0.35 * fish$temperature
fish$growth <- rnorm(n, mean = eta_mu, sd = exp(eta_sigma))

fit_growth <- drmTMB(
  bf(growth ~ temperature * habitat, sigma ~ temperature),
  family = gaussian(),
  data = fish
)

grid_growth <- expand.grid(
  temperature = seq(-1.6, 1.6, length.out = 80),
  habitat = levels(fish$habitat)
)
pred_growth <- predict_parameters(
  fit_growth,
  newdata = grid_growth,
  dpar = c("mu", "sigma"),
  conf.int = TRUE
)
pred_growth$panel <- factor(
  ifelse(
    pred_growth$dpar == "mu",
    "Expected growth (mu)",
    "Residual SD (sigma)"
  ),
  levels = c("Expected growth (mu)", "Residual SD (sigma)")
)

pred_mu <- pred_growth[pred_growth$dpar == "mu", , drop = FALSE]

ggplot(fish, aes(temperature, growth, colour = habitat)) +
  geom_point(alpha = 0.35, size = 1.4) +
  geom_ribbon(
    data = pred_mu,
    aes(
      x = temperature,
      y = estimate,
      ymin = conf.low,
      ymax = conf.high,
      fill = habitat
    ),
    inherit.aes = FALSE,
    alpha = 0.16,
    colour = NA
  ) +
  geom_line(data = pred_mu, aes(y = estimate), linewidth = 0.8) +
  scale_drmtmb_colour +
  scale_drmtmb_fill +
  labs(
    title = "Expected growth along a temperature gradient",
    subtitle = "Raw observations, fitted mean, and 95% confidence band",
    x = "Temperature anomaly",
    y = "Growth",
    colour = "Habitat",
    fill = "Habitat"
  ) +
  theme_drmtmb_gallery()

Scatterplot of simulated growth against temperature for two habitats, with fitted mean lines and 95 percent confidence bands for each habitat.

Distributional-parameter surfaces

The same fit can show more than the mean. predict_parameters() returns a long table with one row per distributional parameter, and plot_parameter_surface() draws estimate lines plus finite interval bands when they are present.

plot_parameter_surface(
  pred_growth,
  x = "temperature",
  colour = "habitat",
  dpar = c("mu", "sigma"),
  facet = "panel"
) +
  scale_drmtmb_colour +
  scale_drmtmb_fill +
  labs(
    title = "Location and residual-scale surfaces",
    subtitle = "Both panels use the same prediction grid",
    x = "Temperature anomaly",
    y = "Estimate on response scale",
    colour = "Habitat"
  ) +
  theme_drmtmb_gallery()

Two-panel plot showing fitted expected growth and residual standard deviation surfaces along temperature for reef and kelp habitats.

The mu panel is the expected response. The sigma panel is the fitted residual standard deviation, so an upward line means lower predictability or larger residual heterogeneity at warmer values.

For scale models, the observed-data analogue is not the raw response itself. Raw growth lives on the location scale. A useful scale check is the residual magnitude after subtracting the fitted mu: under a Gaussian model, E(|residual|) = sigma * sqrt(2 / pi). This keeps the observed scatter visible without putting raw response points on a sigma axis.

pred_row_growth <- predict_parameters(
  fit_growth,
  newdata = fish,
  dpar = c("mu", "sigma")
)
row_mu <- pred_row_growth[pred_row_growth$dpar == "mu", c("row", "estimate")]
row_sigma <- pred_row_growth[
  pred_row_growth$dpar == "sigma",
  c("row", "estimate")
]
names(row_mu)[2] <- "fitted_mu"
names(row_sigma)[2] <- "fitted_sigma"
scale_check <- merge(row_mu, row_sigma, by = "row")
scale_check <- cbind(fish[scale_check$row, ], scale_check)
scale_check$abs_residual <- abs(scale_check$growth - scale_check$fitted_mu)

pred_sigma <- pred_growth[pred_growth$dpar == "sigma", , drop = FALSE]
pred_sigma$expected_abs_residual <- pred_sigma$estimate * sqrt(2 / pi)
pred_sigma$quantity <- "Expected absolute residual"

ggplot(scale_check, aes(temperature, abs_residual, colour = habitat)) +
  geom_point(alpha = 0.22, size = 1.2) +
  geom_line(
    data = pred_sigma,
    aes(
      y = expected_abs_residual,
      colour = quantity,
      group = quantity
    ),
    linewidth = 0.9
  ) +
  scale_colour_manual(values = c(
    habitat_palette,
    "Expected absolute residual" = gallery_colour("green")
  )) +
  labs(
    title = "Residual magnitudes carry the scale information",
    subtitle = "The green curve is sigma converted to expected absolute residual size",
    x = "Temperature anomaly",
    y = "|Observed growth - fitted mu|",
    colour = NULL
  ) +
  theme_drmtmb_gallery()

Scatterplot of absolute residual magnitudes against temperature, with a fitted expected absolute residual curve derived from the residual sigma surface.

Confidence intervals are only one slice through the uncertainty surface. A confidence distribution is a frequentist, posterior-like curve over parameter values whose quantiles are confidence limits. The code below keeps that estimate-plus-standard-error contract dependency-free, then draws a compact cloud so the no-effect line, estimate, and interval width are visible together. The scale coefficient is shown as a residual-SD ratio: values above one mean that one unit of temperature increases residual scale.

coef_growth <- summary(fit_growth)$coefficients
confidence_curves <- rbind(
  make_confidence_distribution(
    parameter = "mu temperature slope",
    estimate = coef_growth["mu:temperature", "estimate"],
    std_error = coef_growth["mu:temperature", "std_error"],
    transform = "identity",
    null_value = 0
  ),
  make_confidence_distribution(
    parameter = "sigma temperature SD ratio",
    estimate = coef_growth["sigma:temperature", "estimate"],
    std_error = coef_growth["sigma:temperature", "std_error"],
    transform = "exp",
    null_value = 0
  )
)
confidence_curves$parameter <- factor(
  confidence_curves$parameter,
  levels = c("mu temperature slope", "sigma temperature SD ratio")
)
confidence_curves$density_scaled <- ave(
  confidence_curves$density,
  confidence_curves$parameter,
  FUN = function(x) x / max(x, na.rm = TRUE)
)
confidence_intervals <- rbind(
  make_confidence_intervals(
    parameter = "mu temperature slope",
    estimate = coef_growth["mu:temperature", "estimate"],
    std_error = coef_growth["mu:temperature", "std_error"],
    transform = "identity"
  ),
  make_confidence_intervals(
    parameter = "sigma temperature SD ratio",
    estimate = coef_growth["sigma:temperature", "estimate"],
    std_error = coef_growth["sigma:temperature", "std_error"],
    transform = "exp"
  )
)
confidence_intervals$parameter <- factor(
  confidence_intervals$parameter,
  levels = levels(confidence_curves$parameter)
)
confidence_intervals$level <- factor(
  confidence_intervals$level,
  levels = c("95%", "66%")
)
confidence_intervals <- confidence_intervals[
  order(confidence_intervals$parameter, confidence_intervals$level),
]
confidence_reference <- unique(
  confidence_curves[c("parameter", "null_value", "estimate_value")]
)
confidence_reference$null_label <- c("no slope", "no scale change")
confidence_reference$label_x <- confidence_reference$null_value + c(0.03, 0.02)
confidence_palette <- c(
  "mu temperature slope" = gallery_colour("blue"),
  "sigma temperature SD ratio" = gallery_colour("green")
)

ggplot(confidence_curves, aes(value, density_scaled, fill = parameter)) +
  geom_area(alpha = 0.24, colour = NA) +
  geom_line(aes(colour = parameter), linewidth = 0.65, alpha = 0.85) +
  geom_vline(
    data = confidence_reference,
    aes(xintercept = null_value),
    inherit.aes = FALSE,
    linetype = "dashed",
    colour = "grey45",
    linewidth = 0.45
  ) +
  geom_text(
    data = confidence_reference,
    aes(x = label_x, y = 0.98, label = null_label),
    inherit.aes = FALSE,
    colour = "grey35",
    size = 3,
    hjust = 0,
    vjust = 0
  ) +
  geom_segment(
    data = confidence_intervals,
    aes(
      x = xmin,
      xend = xmax,
      y = 0,
      yend = 0,
      colour = parameter,
      linewidth = level,
      alpha = level
    ),
    inherit.aes = FALSE,
    lineend = "round"
  ) +
  geom_point(
    data = confidence_reference,
    aes(x = estimate_value, y = 0, colour = parameter),
    inherit.aes = FALSE,
    size = 2.4
  ) +
  facet_wrap(~parameter, scales = "free_x") +
  scale_colour_manual(values = confidence_palette) +
  scale_fill_manual(values = confidence_palette) +
  scale_linewidth_manual(values = c("95%" = 1.2, "66%" = 4.8)) +
  scale_alpha_manual(values = c("95%" = 0.50, "66%" = 0.85)) +
  scale_y_continuous(
    breaks = NULL,
    limits = c(-0.03, 1.10),
    expand = expansion(mult = c(0, 0.02))
  ) +
  labs(
    title = "Confidence clouds show the effect away from no change",
    subtitle = "Dense colour marks more compatible values; bars show central 66% and 95% confidence intervals",
    x = "Effect value",
    y = NULL,
    colour = "Parameter"
  ) +
  theme_drmtmb_gallery() +
  theme(
    legend.position = "none",
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )

Two-panel confidence cloud plot for the fitted mu temperature slope and the residual SD ratio temperature effect, with density shapes, no-effect reference lines, estimate dots, and central interval bars.

Distributional parameters are not all residual standard deviations. Shape, inflation, hurdle, and residual-correlation parameters answer different questions and often live on bounded response scales. The next panel set uses three fitted models and the same plotting helper, but each panel names the estimand explicitly.

set.seed(2594)
n_tail <- 2500
tail_dat <- data.frame(tail_pressure = runif(n_tail, -1.5, 1.5))
tail_nu <- 2 + exp(log(4) - tail_dat$tail_pressure)
tail_dat$robust_trait <- 0.35 +
  0.45 * tail_dat$tail_pressure +
  0.55 * stats::rt(n_tail, df = tail_nu)
fit_tail <- drmTMB(
  bf(robust_trait ~ tail_pressure, sigma ~ 1, nu ~ tail_pressure),
  family = student(),
  data = tail_dat
)
tail_grid <- data.frame(tail_pressure = seq(-1.5, 1.5, length.out = 80))
pred_tail <- predict_parameters(
  fit_tail,
  newdata = tail_grid,
  dpar = "nu",
  conf.int = TRUE
)
pred_tail$gradient <- pred_tail$tail_pressure
pred_tail$panel <- "Tail weight (nu)"
pred_tail$model <- "Tail weight"

set.seed(2595)
n_count <- 180
count_dat <- data.frame(detection_difficulty = runif(n_count, -1.2, 1.2))
mu_count <- exp(0.9 + 0.35 * count_dat$detection_difficulty)
zi_count <- plogis(-1.1 + 1.25 * count_dat$detection_difficulty)
count_dat$count <- ifelse(
  runif(n_count) < zi_count,
  0,
  stats::rpois(n_count, mu_count)
)
fit_zi <- drmTMB(
  bf(count ~ detection_difficulty, zi ~ detection_difficulty),
  family = poisson(link = "log"),
  data = count_dat
)
zi_grid <- data.frame(detection_difficulty = seq(-1.2, 1.2, length.out = 80))
pred_zi <- predict_parameters(
  fit_zi,
  newdata = zi_grid,
  dpar = "zi",
  conf.int = TRUE
)
pred_zi$gradient <- pred_zi$detection_difficulty
pred_zi$panel <- "Structural zero (zi)"
pred_zi$model <- "Structural-zero probability"

set.seed(2596)
n_pair <- 160
pair_dat <- data.frame(disturbance = runif(n_pair, -1.4, 1.4))
rho_true <- tanh(-0.15 + 0.75 * pair_dat$disturbance)
z1 <- stats::rnorm(n_pair)
z2 <- stats::rnorm(n_pair)
pair_dat$activity <- 0.25 + 0.35 * pair_dat$disturbance + z1
pair_dat$boldness <- -0.1 +
  0.2 * pair_dat$disturbance +
  rho_true * z1 +
  sqrt(1 - rho_true^2) * z2
fit_pair <- drmTMB(
  bf(
    mu1 = activity ~ disturbance,
    mu2 = boldness ~ disturbance,
    sigma1 ~ 1,
    sigma2 ~ 1,
    rho12 ~ disturbance
  ),
  family = c(gaussian(), gaussian()),
  data = pair_dat
)
rho_grid <- data.frame(disturbance = seq(-1.4, 1.4, length.out = 80))
pred_rho <- predict_parameters(
  fit_pair,
  newdata = rho_grid,
  dpar = "rho12",
  conf.int = TRUE
)
pred_rho$gradient <- pred_rho$disturbance
pred_rho$panel <- "Residual rho12"
pred_rho$model <- "Residual correlation"

pred_extra <- rbind(
  pred_tail[intersect(names(pred_tail), names(pred_zi))],
  pred_zi[intersect(names(pred_tail), names(pred_zi))],
  pred_rho[intersect(names(pred_tail), names(pred_zi))]
)
pred_extra$panel <- factor(
  pred_extra$panel,
  levels = c(
    "Tail weight (nu)",
    "Structural zero (zi)",
    "Residual rho12"
  )
)

plot_parameter_surface(
  pred_extra,
  x = "gradient",
  colour = "model",
  facet = "panel",
  dpar = c("nu", "zi", "rho12"),
  point = FALSE
) +
  geom_hline(
    data = data.frame(
      .drmTMB_plot_facet = factor(
        "Residual rho12",
        levels = levels(pred_extra$panel)
      )
    ),
    aes(yintercept = 0),
    inherit.aes = FALSE,
    linetype = "dashed",
    colour = "grey55"
  ) +
  scale_colour_manual(values = c(
    "Tail weight" = "#0072B2",
    "Structural-zero probability" = "#009E73",
    "Residual correlation" = "#D55E00"
  )) +
  scale_fill_manual(values = c(
    "Tail weight" = "#0072B2",
    "Structural-zero probability" = "#009E73",
    "Residual correlation" = "#D55E00"
  )) +
  labs(
    title = "Shape, inflation, and residual-correlation panels",
    subtitle = "Each panel is a different distributional parameter on its response scale",
    x = "Predictor gradient",
    y = "Estimate on response scale",
    colour = "Parameter",
    fill = "Parameter"
  ) +
  theme_drmtmb_gallery() +
  theme(legend.position = "none")

Three-panel plot showing a Student-t degrees-of-freedom shape parameter, a zero-inflation probability, and a residual rho12 correlation along a predictor gradient with 95 percent confidence bands.

The nu panel is a Student-t degrees-of-freedom parameter: lower values mean heavier tails, and here the fitted model estimates a changing tail shape through nu ~ tail_pressure. The zi panel is a structural-zero probability, not a mean count. The rho12 panel is residual coupling between two responses after their own means and scales have been modelled.

Random-effect and variance-component figures

Random-effect figures need a different visual grammar from residual-scale figures. Residual sigma is the observation-level standard deviation. A group-level SD, such as the SD of site intercepts or temperature slopes, is among-group variation in a fitted random effect.

set.seed(2621)
n_site <- 24
n_site_obs <- 6
site_levels <- seq_len(n_site)
site_random <- data.frame(
  site = factor(site_levels),
  intercept = rnorm(n_site, mean = 0, sd = 0.42),
  slope = rnorm(n_site, mean = 0, sd = 0.24)
)
reef_re <- data.frame(
  site = factor(rep(site_levels, each = n_site_obs), levels = site_levels)
)
reef_re$temperature <- rep(seq(-1, 1, length.out = n_site_obs), n_site)
site_lookup <- as.integer(reef_re$site)
reef_re$growth <- 0.25 +
  0.62 * reef_re$temperature +
  site_random$intercept[site_lookup] +
  site_random$slope[site_lookup] * reef_re$temperature +
  rnorm(nrow(reef_re), mean = 0, sd = 0.34)

fit_re_slopes <- drmTMB(
  bf(growth ~ temperature + (1 + temperature | site), sigma ~ 1),
  family = gaussian(),
  data = reef_re,
  control = drm_control(optimizer = list(eval.max = 160L, iter.max = 160L))
)

component_rows <- summary(fit_re_slopes)$parameters
variance_components <- data.frame(
  component = c(
    "Residual sigma",
    "Site intercept SD",
    "Site temperature-slope SD"
  ),
  component_type = c(
    "Residual scale",
    "Group intercept SD",
    "Group slope SD"
  ),
  estimate = c(
    component_rows["sigma", "estimate"],
    component_rows[
      "sd:mu:(1 + temperature | site):(Intercept)",
      "estimate"
    ],
    component_rows[
      "sd:mu:(1 + temperature | site):temperature",
      "estimate"
    ]
  )
)
variance_components$component <- factor(
  variance_components$component,
  levels = rev(variance_components$component)
)
variance_components$component_type <- factor(
  variance_components$component_type,
  levels = c("Residual scale", "Group intercept SD", "Group slope SD")
)

component_palette <- c(
  "Residual scale" = gallery_colour("blue"),
  "Group intercept SD" = gallery_colour("green"),
  "Group slope SD" = gallery_colour("orange")
)

ggplot(variance_components, aes(estimate, component, colour = component_type)) +
  geom_segment(
    aes(x = 0, xend = estimate, yend = component),
    linewidth = 0.7,
    alpha = 0.55
  ) +
  geom_point(size = 2.8) +
  scale_colour_manual(values = component_palette) +
  labs(
    title = "Residual scale and group-level SDs are different quantities",
    subtitle = "The slope SD measures among-site variation in temperature responses",
    x = "Fitted standard deviation",
    y = NULL,
    colour = "Component"
  ) +
  theme_drmtmb_gallery() +
  guides(colour = "none")

Dot plot comparing residual sigma, among-site intercept SD, and among-site temperature-slope SD on their fitted response-scale units.

The dot plot compares three standard deviations, but only one is residual sigma. The other two are random-effect SDs: how much sites differ in baseline growth, and how much sites differ in the temperature slope. For intervals on these SDs, use profile_targets(fit_re_slopes) to choose direct profile-ready targets rather than treating the quick point plot as an uncertainty analysis.

Readers used to Bayesian model summaries often expect one display that keeps fixed effects, standard deviations, and correlations together. For this frequentist fit, a direct interval plot is clearer than a cloud: the point is the estimate, the thick bar is the central 66% Wald interval, and the thin bar is the 95% Wald interval.

fixed_rows <- summary(fit_re_slopes)$coefficients
parameter_rows <- summary(fit_re_slopes)$parameters
coefficient_specs <- data.frame(
  group = c(
    "Fixed effects",
    "Fixed effects",
    "SD parameters",
    "SD parameters",
    "SD parameters",
    "Correlation"
  ),
  parameter = c(
    "mu intercept",
    "mu temperature slope",
    "residual sigma",
    "site intercept SD",
    "site temperature-slope SD",
    "intercept-slope correlation"
  ),
  estimate = c(
    fixed_rows["mu:(Intercept)", "estimate"],
    fixed_rows["mu:temperature", "estimate"],
    parameter_rows["sigma", "estimate"],
    parameter_rows[
      "sd:mu:(1 + temperature | site):(Intercept)",
      "estimate"
    ],
    parameter_rows[
      "sd:mu:(1 + temperature | site):temperature",
      "estimate"
    ],
    parameter_rows[
      "cor:mu:cor((Intercept),temperature | site)",
      "estimate"
    ]
  ),
  std_error = c(
    fixed_rows["mu:(Intercept)", "std_error"],
    fixed_rows["mu:temperature", "std_error"],
    parameter_rows["sigma", "std_error"],
    parameter_rows[
      "sd:mu:(1 + temperature | site):(Intercept)",
      "std_error"
    ],
    parameter_rows[
      "sd:mu:(1 + temperature | site):temperature",
      "std_error"
    ],
    parameter_rows[
      "cor:mu:cor((Intercept),temperature | site)",
      "std_error"
    ]
  ),
  null_value = c(0, 0, 0, 0, 0, 0)
)
coefficient_levels <- rev(coefficient_specs$parameter)
coefficient_specs$parameter <- factor(
  coefficient_specs$parameter,
  levels = coefficient_levels
)
coefficient_intervals <- do.call(
  rbind,
  lapply(seq_len(nrow(coefficient_specs)), function(i) {
    row <- coefficient_specs[i, ]
    make_confidence_intervals(
      parameter = row$parameter,
      estimate = row$estimate,
      std_error = row$std_error,
      transform = "identity"
    )
  })
)
coefficient_intervals <- merge(
  coefficient_intervals,
  coefficient_specs[c("parameter", "group")],
  by = "parameter"
)
coefficient_intervals$parameter <- factor(
  coefficient_intervals$parameter,
  levels = coefficient_levels
)
coefficient_intervals$level <- factor(
  coefficient_intervals$level,
  levels = c("95%", "66%")
)
coefficient_intervals <- coefficient_intervals[
  order(coefficient_intervals$parameter, coefficient_intervals$level),
]

interval_group_palette <- c(
  "Fixed effects" = gallery_colour("blue"),
  "SD parameters" = gallery_colour("green"),
  "Correlation" = gallery_colour("orange")
)

ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey55") +
  geom_segment(
    data = coefficient_intervals[coefficient_intervals$level == "95%", ],
    aes(
      x = xmin,
      xend = xmax,
      y = parameter,
      yend = parameter,
      colour = group
    ),
    linewidth = 1.0,
    alpha = 0.45,
    lineend = "round"
  ) +
  geom_segment(
    data = coefficient_intervals[coefficient_intervals$level == "66%", ],
    aes(
      x = xmin,
      xend = xmax,
      y = parameter,
      yend = parameter,
      colour = group
    ),
    linewidth = 3.4,
    alpha = 0.72,
    lineend = "round"
  ) +
  geom_point(
    data = coefficient_specs,
    aes(x = estimate, y = parameter, colour = group),
    size = 2.3
  ) +
  scale_colour_manual(values = interval_group_palette) +
  scale_y_discrete(drop = FALSE) +
  labs(
    title = "Coefficient intervals separate fixed effects, SDs, and correlations",
    subtitle = "Points are estimates; thick and thin bars show 66% and 95% Wald intervals",
    x = "Estimate on displayed scale",
    y = NULL,
    colour = "Parameter class"
  ) +
  theme_drmtmb_gallery() +
  theme(
    legend.position = "bottom",
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )

Coefficient interval panel showing fixed effects, residual and group-level standard deviations, and the random intercept-slope correlation, with zero reference line, estimate dots, and 66 and 95 percent interval bars.

A random-slope figure should show repeated site behaviour. The translucent lines below are conditional site-level trajectories: each site borrows strength from the population fit, but sites can still differ in both intercept and temperature slope.

site_intercept <- ranef(fit_re_slopes, "mu")$terms[[
  "(1 + temperature | site):(Intercept)"
]]
site_slope_deviation <- ranef(fit_re_slopes, "mu")$terms[[
  "(1 + temperature | site):temperature"
]]
site_effects <- data.frame(
  site = names(site_slope_deviation),
  intercept_deviation = unname(site_intercept),
  slope_deviation = unname(site_slope_deviation)
)
site_effects$slope_direction <- factor(
  ifelse(
    site_effects$slope_deviation < 0,
    "Lower than average",
    "Higher than average"
  ),
  levels = c("Lower than average", "Higher than average")
)
fixed_re <- summary(fit_re_slopes)$coefficients
site_lines <- expand.grid(
  site = names(site_slope_deviation),
  temperature = seq(-1, 1, length.out = 50),
  KEEP.OUT.ATTRS = FALSE
)
site_lines <- merge(site_lines, site_effects, by = "site")
site_lines$fitted_growth <- (
  fixed_re["mu:(Intercept)", "estimate"] +
    site_lines$intercept_deviation
) + (
  fixed_re["mu:temperature", "estimate"] +
    site_lines$slope_deviation
) * site_lines$temperature
population_line <- data.frame(
  temperature = seq(-1, 1, length.out = 80)
)
population_line$fitted_growth <- fixed_re["mu:(Intercept)", "estimate"] +
  fixed_re["mu:temperature", "estimate"] * population_line$temperature

ggplot() +
  geom_point(
    data = reef_re,
    aes(temperature, growth),
    colour = "grey55",
    alpha = 0.38,
    size = 1.15
  ) +
  geom_line(
    data = site_lines,
    aes(
      temperature,
      fitted_growth,
      group = site,
      colour = slope_direction
    ),
    alpha = 0.48,
    linewidth = 0.55
  ) +
  geom_line(
    data = population_line,
    aes(temperature, fitted_growth),
    colour = gallery_colour("charcoal"),
    linewidth = 1.05
  ) +
  scale_colour_manual(values = c(
    "Lower than average" = gallery_colour("blue"),
    "Higher than average" = gallery_colour("orange")
  )) +
  labs(
    title = "Translucent site trajectories show random slopes",
    subtitle = "Grey points are repeated site observations; the dark line is the population slope",
    x = "Temperature anomaly",
    y = "Growth",
    colour = "Site slope"
  ) +
  theme_drmtmb_gallery()

Repeated-measures random-slope plot with grey observed site data, translucent site-level fitted trajectories, and a darker population-average fitted line.

A direct sd(site) ~ reef_cover formula asks a third question: does a site-level predictor change the among-site SD in expected growth? The prediction table keeps this target as sd(site), not residual sigma.

set.seed(2622)
n_sd_site <- 28
n_sd_each <- 5
site_sd_data <- data.frame(
  site = factor(seq_len(n_sd_site)),
  reef_cover = runif(n_sd_site, -1.2, 1.2)
)
site_sd <- exp(-0.6 + 0.45 * site_sd_data$reef_cover)
site_deviation <- rnorm(n_sd_site, mean = 0, sd = site_sd)
reef_sd <- site_sd_data[rep(seq_len(n_sd_site), each = n_sd_each), ]
reef_sd$temperature <- rep(seq(-1.2, 1.2, length.out = n_sd_each), n_sd_site)
reef_sd$growth <- 0.4 +
  0.5 * reef_sd$temperature +
  site_deviation[as.integer(reef_sd$site)] +
  rnorm(nrow(reef_sd), mean = 0, sd = 0.35)

fit_site_sd <- drmTMB(
  bf(growth ~ temperature + (1 | site), sigma ~ 1, sd(site) ~ reef_cover),
  family = gaussian(),
  data = reef_sd,
  control = drm_control(optimizer = list(eval.max = 180L, iter.max = 180L))
)

grid_site_sd <- prediction_grid(
  fit_site_sd,
  focal = "reef_cover",
  at = list(reef_cover = seq(-1.2, 1.2, length.out = 80))
)
pred_site_sd <- predict_parameters(
  fit_site_sd,
  newdata = grid_site_sd,
  dpar = "sd(site)",
  conf.int = TRUE
)
pred_site_sd$quantity <- "Among-site SD"

plot_parameter_surface(
  pred_site_sd,
  x = "reef_cover",
  colour = "quantity",
  facet = NULL,
  point = FALSE
) +
  scale_colour_manual(values = c(
    "Among-site SD" = gallery_colour("green")
  )) +
  labs(
    title = "Among-site SD surface",
    subtitle = "Requested Wald intervals remain unavailable for direct random-effect SD surfaces",
    x = "Reef cover",
    y = "Among-site SD in expected growth",
    colour = NULL
  ) +
  theme_drmtmb_gallery() +
  guides(colour = "none")

Line plot showing a fitted among-site random-intercept SD surface along reef cover, with no confidence band because direct random-effect SD surfaces do not yet supply Wald intervals.

This line is a fitted random-effect SD surface. It should not be read as a raw growth curve, and it should not be overlaid with raw growth points. The table behind the plot reports conf.status = "wald_unavailable", which is the honest interval boundary until a profile or bootstrap route is attached to these derived surfaces.

Discrete comparisons

For a factor predictor, the same plotting helper uses interval bars rather than ribbons. This is useful for compact tutorial figures and estimated marginal comparisons.

habitat_grid <- data.frame(
  temperature = 0,
  habitat = factor(c("reef", "kelp"), levels = levels(fish$habitat))
)
pred_habitat <- predict_parameters(
  fit_growth,
  newdata = habitat_grid,
  dpar = "mu",
  conf.int = TRUE
)

plot_parameter_surface(
  pred_habitat,
  x = "habitat",
  colour = "habitat",
  line = FALSE,
  point = TRUE
) +
  scale_drmtmb_colour +
  labs(
    title = "Estimated habitat contrast at average temperature",
    subtitle = "Intervals are Wald intervals from predict_parameters()",
    x = NULL,
    y = "Expected growth",
    colour = "Habitat"
  ) +
  theme_drmtmb_gallery()

Point interval plot comparing expected growth between reef and kelp habitats at average temperature.

Categorical by categorical interaction

For two categorical predictors, show every fitted cell and its interval. Use position dodging so the uncertainty intervals do not sit on top of each other. Raw observations can sit behind the intervals when they help readers see the data support for each cell.

set.seed(2592)
n_cell <- 30
reef <- expand.grid(
  habitat = factor(habitat_levels, levels = habitat_levels),
  season = factor(c("dry", "wet"), levels = c("dry", "wet")),
  replicate = seq_len(n_cell)
)
reef$mu <- with(
  reef,
  0.6 +
    0.35 * (habitat == "kelp") -
    0.15 * (season == "wet") +
    0.45 * (habitat == "kelp") * (season == "wet")
)
reef$growth <- rnorm(nrow(reef), mean = reef$mu, sd = 0.45)

fit_cat_cat <- drmTMB(
  bf(growth ~ habitat * season, sigma ~ 1),
  family = gaussian(),
  data = reef
)

grid_cat_cat <- expand.grid(
  habitat = levels(reef$habitat),
  season = levels(reef$season)
)
pred_cat_cat <- predict_parameters(
  fit_cat_cat,
  newdata = grid_cat_cat,
  dpar = "mu",
  conf.int = TRUE
)

raw_cat_cat <- reef

ggplot(pred_cat_cat, aes(habitat, estimate, colour = season)) +
  geom_point(
    data = raw_cat_cat,
    aes(y = growth, colour = season),
    position = position_jitterdodge(
      jitter.width = 0.08,
      jitter.height = 0,
      dodge.width = 0.45
    ),
    alpha = 0.22,
    size = 1.1,
    show.legend = FALSE
  ) +
  geom_pointrange(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(width = 0.45),
    linewidth = 0.65,
    fatten = 2.3
  ) +
  scale_colour_manual(values = season_palette) +
  labs(
    title = "Categorical by categorical interaction",
    subtitle = "Small points are observations; intervals are fitted cell means",
    x = "Habitat",
    y = "Growth",
    colour = "Season"
  ) +
  theme_drmtmb_gallery()
#> Warning: The `fatten` argument of `geom_pointrange()` is deprecated as of ggplot2 4.0.0.
#> ℹ Please use the `size` aesthetic instead.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

Categorical by categorical interaction plot showing raw growth observations and fitted growth intervals for each habitat and season combination.

Estimated marginal means and marginal summaries

For the first supported emmeans path, a fixed-effect univariate mu model can produce estimated marginal means that are easy to turn into publication-style figures. The same rule applies: keep the target and conditioning values visible. These displays are mu summaries, not sigma, fitted-response, random-effect, or bivariate summaries.

emm_habitat <- suppressMessages(
  as.data.frame(
    emmeans::emmeans(
      fit_growth,
      ~ habitat,
      at = list(temperature = 0),
      type = "response"
    )
  )
)

ggplot(emm_habitat, aes(habitat, emmean, colour = habitat)) +
  geom_pointrange(
    aes(ymin = asymp.LCL, ymax = asymp.UCL),
    linewidth = 0.7,
    fatten = 2.4
  ) +
  scale_drmtmb_colour +
  labs(
    title = "Estimated marginal mean growth",
    subtitle = "Supported fixed-effect mu target, conditioned at temperature = 0",
    x = NULL,
    y = "Expected growth",
    colour = "Habitat"
  ) +
  theme_drmtmb_gallery()

Point interval plot of estimated marginal mean growth for reef and kelp habitats conditioned at temperature zero.

emmeans reference grids can also condition on a factor. In the next display, the mu cell means are grouped by season, so the reader sees the fitted habitat contrast inside each season rather than a single averaged contrast.

emm_season <- suppressMessages(
  as.data.frame(
    emmeans::emmeans(
      fit_cat_cat,
      ~ habitat | season,
      type = "response"
    )
  )
)

ggplot(emm_season, aes(habitat, emmean, colour = season)) +
  geom_pointrange(
    aes(ymin = asymp.LCL, ymax = asymp.UCL),
    position = position_dodge(width = 0.35),
    linewidth = 0.65,
    fatten = 2.3
  ) +
  scale_colour_manual(values = season_palette) +
  labs(
    title = "Factor-conditioned estimated marginal means",
    subtitle = "Habitat mu summaries are shown separately for each season",
    x = "Habitat",
    y = "Expected growth",
    colour = "Season"
  ) +
  theme_drmtmb_gallery()

Point interval plot of estimated marginal mean growth for reef and kelp habitats shown separately for dry and wet seasons.

For a numeric interaction, use an explicit reference grid. Here temperature is held at three named values so the habitat contrast does not hide the conditioning rule.

emm_temperature <- suppressMessages(
  as.data.frame(
    emmeans::emmeans(
      fit_growth,
      ~ habitat | temperature,
      at = list(temperature = c(-1, 0, 1)),
      type = "response"
    )
  )
)
emm_temperature$temperature_slice <- factor(
  emm_temperature$temperature,
  levels = c(-1, 0, 1),
  labels = c("cooler", "average", "warmer")
)

ggplot(emm_temperature, aes(temperature_slice, emmean, colour = habitat)) +
  geom_pointrange(
    aes(ymin = asymp.LCL, ymax = asymp.UCL),
    position = position_dodge(width = 0.42),
    linewidth = 0.65,
    fatten = 2.3
  ) +
  scale_drmtmb_colour +
  labs(
    title = "Estimated marginal means across an interaction grid",
    subtitle = "The same fixed-effect mu bridge, evaluated at named temperature slices",
    x = "Temperature slice",
    y = "Expected growth",
    colour = "Habitat"
  ) +
  theme_drmtmb_gallery()

Point interval plot of estimated marginal mean growth for reef and kelp habitats at three temperature values.

Empirical marginal summaries use a different rule. prediction_grid(..., margin = "empirical") crosses a focal value with fitted rows, and marginal_parameters() averages the resulting parameter predictions. The helper itself is still a plug-in point summary, so the interval bars below are a display approximation: they average row-wise Wald prediction limits from predict_parameters(conf.int = TRUE) within each habitat.

empirical_habitat_grid <- prediction_grid(
  fit_growth,
  focal = "habitat",
  at = list(habitat = levels(fish$habitat)),
  margin = "empirical"
)
empirical_habitat <- marginal_parameters(
  fit_growth,
  newdata = empirical_habitat_grid,
  dpar = "mu",
  by = "habitat"
)
empirical_habitat_rowwise <- predict_parameters(
  fit_growth,
  newdata = empirical_habitat_grid,
  dpar = "mu",
  conf.int = TRUE
)
empirical_habitat_wald <- aggregate(
  cbind(conf.low, conf.high) ~ habitat,
  data = empirical_habitat_rowwise,
  FUN = mean
)
empirical_habitat_plot <- empirical_habitat
empirical_habitat_plot$conf.low <- empirical_habitat_wald$conf.low[
  match(empirical_habitat_plot$habitat, empirical_habitat_wald$habitat)
]
empirical_habitat_plot$conf.high <- empirical_habitat_wald$conf.high[
  match(empirical_habitat_plot$habitat, empirical_habitat_wald$habitat)
]
empirical_habitat_plot$conf.status <- "wald"
empirical_habitat_plot$interval_source <- "averaged_rowwise_wald"

plot_parameter_surface(
  empirical_habitat_plot,
  x = "habitat",
  colour = "habitat",
  facet = NULL,
  line = FALSE,
  point = TRUE
) +
  scale_drmtmb_colour +
  labs(
    title = "Empirical marginal mu summary",
    subtitle = "Bars average row-wise Wald prediction limits over fitted rows",
    x = NULL,
    y = "Average expected growth",
    colour = "Habitat",
    caption = "marginal_parameters() returns plug-in point averages; bars are a display approximation."
  ) +
  theme_drmtmb_gallery()

Point-and-interval plot of empirical marginal expected growth for reef and kelp habitats, where points are plug-in marginal means and bars average row-wise Wald prediction limits.

Unsupported targets should be visible before a reader tries to plot them. The current emmeans bridge is for fixed-effect univariate mu; other targets should use parameter prediction, corpairs(), profile intervals, or a future feature-specific design.

emmeans_boundary <- data.frame(
  target = c(
    "Fixed-effect univariate mu",
    "sigma",
    "Bivariate responses",
    "Zero-inflated or hurdle means",
    "Ordinal expected scores",
    "Random-effect targets"
  ),
  status = c(
    "Supported emmeans route",
    "Unsupported boundary",
    "Unsupported boundary",
    "Unsupported boundary",
    "Unsupported boundary",
    "Unsupported boundary"
  )
)
emmeans_boundary$target <- factor(
  emmeans_boundary$target,
  levels = rev(emmeans_boundary$target)
)
emmeans_boundary$status <- factor(
  emmeans_boundary$status,
  levels = c("Supported emmeans route", "Unsupported boundary")
)
emmeans_boundary$status_label <- ifelse(
  emmeans_boundary$status == "Supported emmeans route",
  "supported\nmu EMM",
  "unsupported\nboundary"
)
emmeans_boundary$status_text_colour <- ifelse(
  emmeans_boundary$status == "Supported emmeans route",
  "white",
  "grey10"
)

ggplot(emmeans_boundary, aes(status, target, fill = status)) +
  geom_tile(width = 0.84, height = 0.72, colour = "white", linewidth = 0.6) +
  geom_text(
    aes(label = status_label, colour = status_text_colour),
    size = 3.1,
    lineheight = 0.9
  ) +
  scale_fill_manual(values = c(
    "Supported emmeans route" = "#0072B2",
    "Unsupported boundary" = "#E69F00"
  )) +
  scale_colour_identity() +
  labs(
    title = "Current emmeans support boundary",
    subtitle = "Only fixed-effect univariate mu is advertised as an emmeans target",
    x = "Current status",
    y = NULL,
    fill = "Status"
  ) +
  theme_drmtmb_gallery() +
  theme(
    panel.grid = element_blank(),
    legend.position = "none"
  )

Tile plot showing fixed-effect univariate mu as the supported emmeans target and sigma, bivariate, zero-inflated, hurdle, ordinal, and random-effect targets as unsupported boundaries.

Continuous by continuous interaction

For two continuous predictors, choose scientifically meaningful slices of one predictor and show fitted slopes across the other. The values should be named in the legend so the reader can interpret the conditioning.

set.seed(2593)
n_cc <- 180
soil <- data.frame(
  temperature = runif(n_cc, -1.5, 1.5),
  moisture = runif(n_cc, -1.4, 1.4)
)
soil$growth <- rnorm(
  n_cc,
  mean = 0.3 +
    0.55 * soil$temperature +
    0.45 * soil$moisture -
    0.35 * soil$temperature * soil$moisture,
  sd = 0.45
)

fit_cont_cont <- drmTMB(
  bf(growth ~ temperature * moisture, sigma ~ 1),
  family = gaussian(),
  data = soil
)

grid_cont_cont <- expand.grid(
  temperature = seq(-1.5, 1.5, length.out = 80),
  moisture = c(-1, 0, 1)
)
grid_cont_cont$moisture_level <- factor(
  grid_cont_cont$moisture,
  levels = c(-1, 0, 1),
  labels = c("low moisture", "average moisture", "high moisture")
)

pred_cont_cont <- predict_parameters(
  fit_cont_cont,
  newdata = grid_cont_cont,
  dpar = "mu",
  conf.int = TRUE
)

ggplot(
  pred_cont_cont,
  aes(temperature, estimate, colour = moisture_level, fill = moisture_level)
) +
  geom_point(
    data = soil,
    aes(temperature, growth),
    inherit.aes = FALSE,
    alpha = 0.18,
    size = 1,
    colour = "grey35"
  ) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.14, colour = NA) +
  geom_line(linewidth = 0.8) +
  scale_colour_manual(values = c(
    "low moisture" = gallery_colour("blue"),
    "average moisture" = gallery_colour("green"),
    "high moisture" = gallery_colour("orange")
  )) +
  scale_fill_manual(values = c(
    "low moisture" = gallery_colour("blue"),
    "average moisture" = gallery_colour("green"),
    "high moisture" = gallery_colour("orange")
  )) +
  labs(
    title = "Continuous by continuous interaction",
    subtitle = "Moisture is held at three named slices",
    x = "Temperature anomaly",
    y = "Growth",
    colour = "Moisture slice",
    fill = "Moisture slice"
  ) +
  theme_drmtmb_gallery()

Continuous by continuous interaction plot showing raw observations in grey and fitted temperature slopes at low, average, and high moisture with confidence bands.

The raw point cloud is not coloured by the sliced predictor because the fitted lines are conditional displays, not clusters in the data. This keeps the visual message honest: the model is being read at three moisture values, while the observations span the whole moisture gradient.

Fitted correlation summaries

Correlation displays should keep layers separate. Residual rho12, ordinary-group correlations, phylogenetic correlations, spatial correlations, and future animal or relmat() layers are not the same parameter even when all are correlations. A fitted workflow should start from corpairs(fit); this gallery uses a compact compatible table so the page builds quickly while preserving the same required columns.

pair_table <- data.frame(
  level = c("residual", "group", "phylogenetic"),
  class = c("residual", "mean-slope", "mean-mean"),
  parameter = c(
    "rho12",
    "cor((Intercept),temperature | site)",
    "cor(mu1,mu2 | phylo | species)"
  ),
  estimate = c(-0.32, 0.44, 0.61),
  modelled = c(TRUE, FALSE, FALSE),
  conf.low = c(-0.54, 0.12, 0.29),
  conf.high = c(-0.05, 0.68, 0.79),
  conf.status = c("wald", "profile", "profile"),
  layer = c("Residual coscale", "Ordinary group", "Phylogenetic")
)
pair_table$layer <- factor(
  pair_table$layer,
  levels = c("Residual coscale", "Ordinary group", "Phylogenetic")
)
pair_table$display <- factor(
  c(
    "Residual\nrho12",
    "Group\nintercept-slope",
    "Phylogenetic\nlocation-location"
  ),
  levels = c(
    "Residual\nrho12",
    "Group\nintercept-slope",
    "Phylogenetic\nlocation-location"
  )
)

correlation_palette <- c(
  "Residual coscale" = "#0072B2",
  "Ordinary group" = "#009E73",
  "Phylogenetic" = "#D55E00"
)

plot_corpairs(pair_table, colour = "layer", facet = "layer", label = "display") +
  scale_colour_manual(values = correlation_palette) +
  labs(
    title = "Implemented correlation rows stay in their own layers",
    subtitle = "Intervals are drawn only when the pair table has finite bounds",
    colour = "Layer"
  ) +
  theme_drmtmb_gallery() +
  theme(
    axis.text.y = element_text(size = 10, lineheight = 0.95),
    legend.position = "none"
  )

Faceted point interval plot showing separate residual, group, and phylogenetic correlation layers with short row labels.

The residual row is the fitted within-observation coscale rho12. The group row is a latent random-effect correlation, here an intercept-slope row from an ordinary grouped model. The phylogenetic row is a species-level structured location-location correlation. These rows can live in the same corpairs() table, but the figure keeps them apart because they answer different biological questions.

layer_boundary <- data.frame(
  layer = c(
    "Residual rho12",
    "Ordinary group",
    "Phylogenetic",
    "Spatial",
    "Animal",
    "relmat"
  ),
  status = c(
    "Fitted corpairs row",
    "Fitted corpairs row",
    "Fitted corpairs row",
    "Planned boundary",
    "Planned boundary",
    "Planned boundary"
  )
)
layer_boundary$layer <- factor(
  layer_boundary$layer,
  levels = rev(layer_boundary$layer)
)
layer_boundary$status <- factor(
  layer_boundary$status,
  levels = c("Fitted corpairs row", "Planned boundary")
)
layer_boundary$status_label <- ifelse(
  layer_boundary$status == "Fitted corpairs row",
  "fitted\nrow",
  "planned\nboundary"
)
layer_boundary$status_text_colour <- ifelse(
  layer_boundary$status == "Fitted corpairs row",
  "white",
  "grey10"
)

ggplot(layer_boundary, aes(status, layer, fill = status)) +
  geom_tile(width = 0.82, height = 0.72, colour = "white", linewidth = 0.6) +
  geom_text(
    aes(label = status_label, colour = status_text_colour),
    size = 3.1,
    lineheight = 0.9
  ) +
  scale_fill_manual(values = c(
    "Fitted corpairs row" = "#0072B2",
    "Planned boundary" = "#E69F00"
  )) +
  scale_colour_identity() +
  labs(
    title = "Correlation-layer support boundary",
    subtitle = "Spatial, animal, and relmat rows are reserved until fitted correlation-pair evidence exists",
    x = "Current status",
    y = NULL,
    fill = "Status"
  ) +
  theme_drmtmb_gallery() +
  theme(
    panel.grid = element_blank(),
    legend.position = "none"
  )

Tile plot showing implemented residual, group, and phylogenetic correlation-pair layers separately from planned spatial, animal, and relmat boundaries.

The status strip is not an estimate plot. Coordinate-spatial Gaussian mu slopes are fitted elsewhere, but spatial corpairs() rows remain planned. animal() and relmat() are also reserved structured-effect markers until the likelihood, diagnostics, profile-target, recovery-test, and example evidence exist.

Simulation operating characteristics

Simulation figures answer different questions from model-interpretation figures. They should show bias, error, coverage, runtime, and failures beside the model surface being tested.

sim_summary <- data.frame(
  surface = rep(
    c("Gaussian location-scale", "Proportion", "Count"),
    each = 4
  ),
  estimand = rep(c("beta_x", "sigma", "sd_intercept", "rho12"), 3),
  bias = c(
    0.01, -0.02, 0.04, 0.00,
    -0.01, 0.03, NA, NA,
    0.02, -0.01, -0.04, NA
  ),
  rmse = c(
    0.07, 0.09, 0.13, 0.11,
    0.08, 0.10, NA, NA,
    0.10, 0.09, 0.16, NA
  ),
  coverage = c(
    0.95, 0.94, 0.91, 0.93,
    0.94, 0.92, NA, NA,
    0.93, 0.95, 0.89, NA
  )
)
sim_summary$n_rep <- 500L
sim_summary$bias_mcse <- sim_summary$rmse / sqrt(sim_summary$n_rep)
sim_summary$bias_low <- sim_summary$bias - 1.96 * sim_summary$bias_mcse
sim_summary$bias_high <- sim_summary$bias + 1.96 * sim_summary$bias_mcse
sim_summary$coverage_mcse <- sqrt(
  sim_summary$coverage * (1 - sim_summary$coverage) / sim_summary$n_rep
)
sim_summary$coverage_low <- pmax(
  0,
  sim_summary$coverage - 1.96 * sim_summary$coverage_mcse
)
sim_summary$coverage_high <- pmin(
  1,
  sim_summary$coverage + 1.96 * sim_summary$coverage_mcse
)
sim_summary$surface <- factor(
  sim_summary$surface,
  levels = names(surface_palette)
)
estimand_labels <- c(
  "beta_x" = "beta_x fixed slope",
  "sigma" = "sigma residual SD",
  "sd_intercept" = "sd_intercept group SD",
  "rho12" = "rho12 residual correlation"
)
sim_summary$estimand_label <- factor(
  estimand_labels[sim_summary$estimand],
  levels = rev(unname(estimand_labels))
)
display_replicates <- 45
sim_rain <- do.call(
  rbind,
  lapply(seq_len(nrow(sim_summary)), function(i) {
    row <- sim_summary[i, ]
    if (is.na(row$bias)) {
      return(NULL)
    }
    data.frame(
      surface = row$surface,
      estimand_label = row$estimand_label,
      bias_replicate = row$bias +
        qnorm(ppoints(display_replicates)) * row$rmse
    )
  })
)
sim_rain$surface <- factor(sim_rain$surface, levels = levels(sim_summary$surface))
sim_rain$estimand_label <- factor(
  sim_rain$estimand_label,
  levels = levels(sim_summary$estimand_label)
)

bias_position <- position_dodge(width = 0.65)

ggplot() +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60") +
  geom_violin(
    data = sim_rain,
    aes(
      estimand_label,
      bias_replicate,
      group = interaction(estimand_label, surface),
      fill = surface
    ),
    position = bias_position,
    width = 0.62,
    alpha = 0.16,
    colour = NA,
    scale = "width",
    trim = FALSE
  ) +
  geom_point(
    data = sim_rain,
    aes(estimand_label, bias_replicate, colour = surface),
    position = position_jitterdodge(jitter.width = 0.08, dodge.width = 0.65),
    alpha = 0.22,
    size = 0.75,
    show.legend = FALSE
  ) +
  geom_linerange(
    data = sim_summary,
    aes(
      estimand_label,
      ymin = bias_low,
      ymax = bias_high,
      colour = surface
    ),
    position = bias_position,
    linewidth = 0.95,
    na.rm = TRUE
  ) +
  geom_point(
    data = sim_summary,
    aes(estimand_label, bias, colour = surface),
    position = bias_position,
    size = 2.6,
    na.rm = TRUE
  ) +
  coord_flip() +
  scale_colour_manual(values = surface_palette, drop = FALSE) +
  scale_fill_manual(values = surface_palette, drop = FALSE) +
  labs(
    title = "Bias by estimand",
    subtitle = "Rainclouds show replicate errors; large points are mean bias with 95% MCSE intervals",
    x = NULL,
    y = "Estimate minus truth",
    colour = "Surface"
  ) +
  guides(fill = "none") +
  theme_drmtmb_gallery()

Simulation summary plots showing illustrative replicate-level bias clouds with mean and MCSE intervals, plus interval coverage bars with binomial MCSE whiskers, by estimand across Gaussian location-scale, proportion, and count surfaces.


coverage_position <- position_dodge2(width = 0.78, preserve = "single")

ggplot(sim_summary, aes(estimand_label, coverage, fill = surface)) +
  geom_hline(yintercept = 0.95, linetype = "dashed", colour = "grey60") +
  geom_col(
    position = coverage_position,
    width = 0.65,
    colour = "white",
    linewidth = 0.25,
    na.rm = TRUE
  ) +
  geom_errorbar(
    aes(ymin = coverage_low, ymax = coverage_high, colour = surface),
    position = coverage_position,
    width = 0.22,
    linewidth = 0.45,
    na.rm = TRUE
  ) +
  coord_flip(ylim = c(0.75, 1)) +
  scale_fill_manual(values = surface_palette, drop = FALSE) +
  scale_colour_manual(values = surface_palette, drop = FALSE, guide = "none") +
  labs(
    title = "Interval coverage",
    subtitle = "Bars are means; whiskers are approximate 95% binomial MCSE intervals",
    x = NULL,
    y = "Empirical coverage",
    fill = "Surface"
  ) +
  theme_drmtmb_gallery()

These figures use a small illustrative table with the same columns that simulation summaries should carry. The bias panel shows the intended raincloud-style grammar: real reports should draw the replicate-level estimates, then overplot the mean and MCSE interval for each defined surface-estimand cell. The panel still leaves estimands disconnected because beta_x, sigma, sd_intercept, and rho12 are categorical targets, not steps along one trajectory. The comprehensive simulation programme should not be count-only. It should add continuous, meta-analysis, proportion, count, ordinal, bivariate, random-slope, and structured-dependence surfaces as each surface becomes ready. Real simulation reports should also show the run manifest and warning/error ledger so missing or failed fits are not hidden behind clean-looking panels.

Gallery source map

The gallery is easier to audit when every figure names its source object, extractor, interval source, and support boundary. Use this table as a quick check before adapting a plot to a manuscript, tutorial, or simulation report.

gallery_source_map <- data.frame(
  figure = c(
    "Mean with confidence band",
    "mu and sigma surfaces",
    "nu, zi, and rho12 panels",
    "Variance-component dot plot",
    "Conditional random-slope deviations",
    "sd(site) surface",
    "Discrete habitat comparison",
    "Categorical interaction",
    "Simple emmeans display",
    "Factor-conditioned emmeans grid",
    "Interaction emmeans grid",
    "Empirical marginal mu summary",
    "emmeans support boundary",
    "Continuous interaction",
    "Correlation estimate rows",
    "Correlation support boundary",
    "Simulation bias and coverage"
  ),
  source = c(
    "fit_growth",
    "fit_growth",
    "fit_tail, fit_zi, fit_pair",
    "fit_re_slopes",
    "fit_re_slopes",
    "fit_site_sd",
    "fit_growth",
    "fit_cat_cat",
    "fit_growth",
    "fit_cat_cat",
    "fit_growth",
    "fit_growth",
    "emmeans_boundary",
    "fit_cont_cont",
    "pair_table",
    "layer_boundary",
    "sim_summary"
  ),
  extractor_or_plotter = c(
    "predict_parameters(dpar = \"mu\") + ggplot2",
    "predict_parameters() + plot_parameter_surface()",
    "predict_parameters() + plot_parameter_surface()",
    "summary(fit)$parameters + ggplot2",
    "ranef() + ggplot2",
    "prediction_grid() + predict_parameters() + plot_parameter_surface()",
    "predict_parameters(dpar = \"mu\") + plot_parameter_surface()",
    "predict_parameters(dpar = \"mu\") + ggplot2",
    "emmeans::emmeans(type = \"response\") + ggplot2",
    "emmeans::emmeans(~ habitat | season) + ggplot2",
    "emmeans::emmeans(~ habitat | temperature) + ggplot2",
    "prediction_grid(margin = \"empirical\") + marginal_parameters() + predict_parameters()",
    "ggplot2 status strip",
    "predict_parameters(dpar = \"mu\") + ggplot2",
    "plot_corpairs()",
    "ggplot2 status strip",
    "ggplot2 fixture display"
  ),
  interval_source = c(
    "Wald intervals from predict_parameters()",
    "Wald intervals when finite",
    "Wald intervals from each fitted parameter table",
    "Point estimates only",
    "Conditional modes only",
    "conf.status = \"wald_unavailable\"",
    "Wald intervals from predict_parameters()",
    "Wald intervals from predict_parameters()",
    "asymptotic emmeans intervals",
    "asymptotic emmeans intervals",
    "asymptotic emmeans intervals",
    "Display approximation from averaged row-wise Wald limits",
    "No intervals; status only",
    "Wald intervals from predict_parameters()",
    "Finite bounds supplied in compatible pair table",
    "No intervals; status only",
    "Illustrative summaries, not final evidence"
  ),
  support_boundary = c(
    "fixed-effect univariate mu",
    "fixed-effect mu and sigma",
    "fixed-effect shape, inflation, and residual rho12",
    "ordinary Gaussian random-slope fit",
    "shrunken deviations, not variance intervals",
    "direct random-effect SD surface without Wald bands",
    "fixed-effect univariate mu",
    "fixed-effect univariate mu",
    "fixed-effect univariate mu only",
    "fixed-effect univariate mu only",
    "fixed-effect univariate mu only",
    "empirical marginalisation, not emmeans",
    "unsupported targets remain unsupported",
    "fixed-effect univariate mu",
    "implemented residual, group, and phylogenetic rows",
    "spatial, animal, and relmat rows remain planned",
    "display grammar for later result articles"
  ),
  stringsAsFactors = FALSE
)

knitr::kable(gallery_source_map)

figure	source	extractor_or_plotter	interval_source	support_boundary
Mean with confidence band	fit_growth	predict_parameters(dpar = “mu”) + ggplot2	Wald intervals from predict_parameters()	fixed-effect univariate mu
mu and sigma surfaces	fit_growth	predict_parameters() + plot_parameter_surface()	Wald intervals when finite	fixed-effect mu and sigma
nu, zi, and rho12 panels	fit_tail, fit_zi, fit_pair	predict_parameters() + plot_parameter_surface()	Wald intervals from each fitted parameter table	fixed-effect shape, inflation, and residual rho12
Variance-component dot plot	fit_re_slopes	summary(fit)$parameters + ggplot2	Point estimates only	ordinary Gaussian random-slope fit
Conditional random-slope deviations	fit_re_slopes	ranef() + ggplot2	Conditional modes only	shrunken deviations, not variance intervals
sd(site) surface	fit_site_sd	prediction_grid() + predict_parameters() + plot_parameter_surface()	conf.status = “wald_unavailable”	direct random-effect SD surface without Wald bands
Discrete habitat comparison	fit_growth	predict_parameters(dpar = “mu”) + plot_parameter_surface()	Wald intervals from predict_parameters()	fixed-effect univariate mu
Categorical interaction	fit_cat_cat	predict_parameters(dpar = “mu”) + ggplot2	Wald intervals from predict_parameters()	fixed-effect univariate mu
Simple emmeans display	fit_growth	emmeans::emmeans(type = “response”) + ggplot2	asymptotic emmeans intervals	fixed-effect univariate mu only
Factor-conditioned emmeans grid	fit_cat_cat	emmeans::emmeans(~ habitat \| season) + ggplot2	asymptotic emmeans intervals	fixed-effect univariate mu only
Interaction emmeans grid	fit_growth	emmeans::emmeans(~ habitat \| temperature) + ggplot2	asymptotic emmeans intervals	fixed-effect univariate mu only
Empirical marginal mu summary	fit_growth	prediction_grid(margin = “empirical”) + marginal_parameters() + predict_parameters()	Display approximation from averaged row-wise Wald limits	empirical marginalisation, not emmeans
emmeans support boundary	emmeans_boundary	ggplot2 status strip	No intervals; status only	unsupported targets remain unsupported
Continuous interaction	fit_cont_cont	predict_parameters(dpar = “mu”) + ggplot2	Wald intervals from predict_parameters()	fixed-effect univariate mu
Correlation estimate rows	pair_table	plot_corpairs()	Finite bounds supplied in compatible pair table	implemented residual, group, and phylogenetic rows
Correlation support boundary	layer_boundary	ggplot2 status strip	No intervals; status only	spatial, animal, and relmat rows remain planned
Simulation bias and coverage	sim_summary	ggplot2 fixture display	Illustrative summaries, not final evidence	display grammar for later result articles

Future simulation result articles

Specialised simulation result articles should come later, after the comprehensive simulation surfaces are ready. Those articles can be organised by scientific question: power, bias, coverage, runtime, convergence, and failure patterns. They should also be organised across data types, not only counts: continuous location-scale models, meta-analysis with known V, proportions, counts, ordinal responses, bivariate rho12, random-slope surfaces, and structured dependence.

This figure gallery is broader and lighter. It shows the visual language users should see across tutorials, examples, model interpretation, and later simulation reports.