Input Saved prediction and summary files.
Action Read the model record and build its figure.
Output Model table, checks, and plot.
Next Continue to the next model chapter.
This model contains 253 species (7,186 contrasts). The phylogenetic correlation matrix was constructed and aligned to those species as follows; brms uses it in (1 | gr(sp_ncbi, cov = A)).
A_full <- prepR4pcm::pr_phylo_cor(tree_final)
# Fallback if needed: ape::vcv.phylo(tree_final, corr = TRUE)
dat_model <- dat_es
phylo <- prepare_phylo_and_data(dat_model, A_full, fill_unmatched = TRUE)
dat_model <- phylo$dat_model
A <- phylo$A_mod
formula <- brms::bf(
yi_lnM_safe ~ 1 + (1 | ref_id) +
(1 | gr(sp_ncbi, cov = A)) + (1 | gr(es_id_model, cov = V)),
sigma ~ 1
)
V <- phylo$V
fit <- brms::brm(formula, data = dat_model,
data2 = list(A = A, V = V))The intercept-only model estimates the overall mean lnM and the baseline residual heterogeneity (sigma), with no moderators — the reference against which every moderator model is compared.
summary() is read from Rdata/epred_draws/m00.rds; the full brms fit is not loaded at render.
print_cache_summary(cache) Family: gaussian
Links: mu = identity; sigma = log
Formula: yi_lnM_safe ~ 1 + (1 | ref_id) + (1 | gr(sp_ncbi, cov = A)) + (1 | gr(es_id_model, cov = V))
sigma ~ 1
Data: structure(list(es_id = c("es0313", "es0314", "es03 (Number of observations: 7186)
Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup draws = 8000
Multilevel Hyperparameters:
~es_id_model (Number of levels: 7186)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
~ref_id (Number of levels: 254)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.74 0.04 0.66 0.83 1.00 1737 3017
~sp_ncbi (Number of levels: 253)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.75 0.17 0.43 1.11 1.00 810 1287
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.55 0.43 -1.41 0.30 1.00 2489 3073
sigma_Intercept -0.59 0.01 -0.61 -0.56 1.00 3469 4817
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The model summary reports the scale-model intercept on the log-\(\sigma\) scale, because the brms model uses Links: mu = identity; sigma = log. The summary therefore gives
\[ \beta_{\mathrm{Overall}}^{(s)} = -0.59, \]
which is really \(\log(\sigma) = -0.59\), not \(\sigma\) itself.
The orchard plot below explicitly back-transforms that coefficient with exp() before plotting it:
\[ \sigma = \exp(-0.59) \approx 0.55. \]
The plotting code confirms this directly (sig_est <- exp(fe["sigma_Intercept", "Estimate"])), and the scale panel’s x-axis is labelled “residual lnM (SD)”, not log-sigma.
So the positive value shown in the orchard plot’s scale panel is not a residual \(\ln M\) effect size, and it is not showing logged residuals. It is the residual standard deviation of \(\ln M\) — \(\sigma\) — on its natural scale. Because a standard deviation cannot be negative, it appears positive even though the corresponding row in the summary table above is negative. This applies to every scale-submodel coefficient in this book, not only the intercept.
The overall posterior mean is shown on the \(\ln M\) scale and as \(d_{\mathrm{eq}} \approx \sqrt{2\exp(2\ln M)}\). The transformation is applied to the posterior 95% credible-interval endpoints.
source(here::here("Scripts", "18_contrasts_cache.R"))
kable_contrasts(
format_posterior_d_eq(cache$post$b_Intercept, "Overall mean"),
caption = "Overall mean lnM and equivalent d with 95% credible intervals."
)| Estimate | Mean lnM | lnM lower 95% CrI | lnM upper 95% CrI | deq | deq lower 95% CrI | deq upper 95% CrI |
|---|---|---|---|---|---|---|
| Overall mean | -0.553 | -1.415 | 0.305 | 0.814 | 0.344 | 1.918 |
Location–scale \(I^2\) from this intercept-only baseline, partitioned into between-study (reference), phylogenetic, and within-study (residual) components following the location–scale meta-analysis framework of Nakagawa et al. Values are read from Rdata/summaries/heterogeneity_m00.rds.
kable_contrasts(format_heterogeneity(read_heterogeneity()),
caption = "Location–scale heterogeneity (I², %) with 95% credible intervals.")| Component | I² (%) | Lower 95% CrI | Upper 95% CrI |
|---|---|---|---|
| Total | 96.6 | 95.5 | 97.7 |
| Between-study (reference) | 37.3 | 23.8 | 52.4 |
| Phylogeny | 38.1 | 16.2 | 58.9 |
| Within-study (residual) | 21.0 | 14.6 | 28.4 |
For each posterior draw, the heterogeneity denominator is:
\[ D = \sigma_u^2 + \sigma_{\mathrm{phylo}}^2 + \bar{\sigma}_e^2 + \bar{V}, \]
where \(\sigma_u^2\) is the between-reference variance, \(\sigma_{\mathrm{phylo}}^2\) is the phylogenetic variance, \(\bar{\sigma}_e^2 = \exp(2 b_{\sigma,0})\) is the residual heterogeneity variance from the scale submodel, and \(\bar{V}\) is the typical known sampling variance. Component-specific \(I^2\) values are:
\[ I_x^2 = 100 \times \frac{\sigma_x^2}{D}, \]
with total heterogeneity:
\[ I_{\mathrm{total}}^2 = 100 \times \frac{\sigma_u^2 + \sigma_{\mathrm{phylo}}^2 + \bar{\sigma}_e^2}{D}. \]
The typical sampling variance is computed from the diagonal sampling-variance matrix \(V\) as:
\[ \bar{V} = \frac{n - p}{\operatorname{tr}(P)}, \qquad P = W - W X (X^\top W X)^{-1} X^\top W, \qquad W = V^{-1}. \]
source(here::here("Scripts", "18_contrasts_cache.R"))
kable_contrasts(format_heterogeneity_components(read_heterogeneity()),
caption = "Precomputed variance components used in the m00 heterogeneity denominator.")| Quantity | Value |
|---|---|
| Vbar (typical sampling variance) | 0.050 |
| sigma2_u (between-study, ref_id) | 0.547 |
| sigma2_phylo (phylogeny, sp_ncbi) | 0.559 |
| sigma2bar_e (within-study residual) | 0.309 |
The orchard-style figure is read from Rdata/figures/publication/orchard/. The complete editable code used to build this plot is shown below.
knitr::include_graphics(
here::here("Rdata", "figures", "publication", "orchard", "m00_orchard_combined.png")
)
library(here)
library(tidyverse)
library(brms)
library(ggbeeswarm)
library(patchwork)
out_dir <- here::here("Rdata", "figures", "publication", "orchard")
out_dir_pdf <- here::here("Figures", "publication", "orchard")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(out_dir_pdf, recursive = TRUE, showWarnings = FALSE)
dat_es <- readRDS(here::here("Rdata", "effect_sizes", "proceed_lnm_safe.rds"))
fit <- readRDS(here::here("Rdata", "models", "m00_ls_intercept_only.rds"))
theme_orchard <- function() {
ggplot2::theme_classic(base_size = 13) +
ggplot2::theme(
axis.text.y = ggplot2::element_text(size = 11),
axis.title.x = ggplot2::element_text(size = 12),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.x = ggplot2::element_line(colour = "grey92"),
legend.position = "bottom",
legend.title = ggplot2::element_text(size = 10),
plot.title = ggplot2::element_text(size = 13, face = "bold"),
plot.caption = ggplot2::element_text(size = 9, colour = "grey50")
)
}
save_plot <- function(p, stem, width = 9, height = 6) {
ggplot2::ggsave(file.path(out_dir_pdf, paste0(stem, ".pdf")), p,
width = width, height = height, device = cairo_pdf)
ggplot2::ggsave(file.path(out_dir, paste0(stem, ".png")), p,
width = width, height = height, dpi = 300, type = "cairo")
}
raw <- dat_es |>
dplyr::mutate(precision = 1 / sqrt(vi_lnM_safe))
fe <- brms::fixef(fit)
int_est <- fe["Intercept", "Estimate"]
int_lo <- fe["Intercept", "Q2.5"]
int_hi <- fe["Intercept", "Q97.5"]
sig_est <- exp(fe["sigma_Intercept", "Estimate"])
sig_lo <- exp(fe["sigma_Intercept", "Q2.5"])
sig_hi <- exp(fe["sigma_Intercept", "Q97.5"])
vc <- brms::VarCorr(fit)
caption_txt <- sprintf(
"n = %d effect sizes | study SD = %.2f | phylogeny SD = %.2f",
nrow(raw), vc$ref_id$sd["Intercept", "Estimate"],
vc$sp_ncbi$sd["Intercept", "Estimate"]
)
p_loc <- ggplot2::ggplot() +
ggplot2::geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50") +
ggbeeswarm::geom_quasirandom(
data = raw,
ggplot2::aes(x = yi_lnM_safe, y = 1, size = precision),
alpha = 0.20, shape = 21, fill = "#88CCEE", colour = "#0072B2",
groupOnX = FALSE
) +
ggplot2::geom_linerange(
ggplot2::aes(y = 1, xmin = int_lo, xmax = int_hi),
linewidth = 1.5, colour = "grey15"
) +
ggplot2::geom_point(
ggplot2::aes(x = int_est, y = 1),
size = 5, shape = 21, fill = "white", colour = "grey10", stroke = 1.2
) +
ggplot2::scale_size_continuous(name = "Precision (1/SE)", range = c(0.3, 4)) +
ggplot2::scale_y_continuous(breaks = NULL) +
ggplot2::labs(
x = "Location effect (lnM)",
y = NULL,
title = "Location -- Overall baseline",
caption = caption_txt
) +
theme_orchard()
raw_sig <- raw |>
dplyr::mutate(abs_resid = abs(yi_lnM_safe - int_est))
p_scl <- ggplot2::ggplot() +
ggplot2::geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50") +
ggbeeswarm::geom_quasirandom(
data = raw_sig,
ggplot2::aes(x = abs_resid, y = 1, size = precision),
alpha = 0.20, shape = 21, fill = "#88CCEE", colour = "#0072B2",
groupOnX = FALSE
) +
ggplot2::geom_linerange(
ggplot2::aes(y = 1, xmin = sig_lo, xmax = sig_hi),
linewidth = 3, colour = "white"
) +
ggplot2::geom_linerange(
ggplot2::aes(y = 1, xmin = sig_lo, xmax = sig_hi),
linewidth = 1.5, colour = "grey15"
) +
ggplot2::geom_point(
ggplot2::aes(x = sig_est, y = 1),
size = 5, shape = 21, fill = "white", colour = "grey10", stroke = 1.2
) +
ggplot2::scale_size_continuous(name = "Precision (1/SE)", range = c(0.3, 4)) +
ggplot2::scale_y_continuous(breaks = NULL) +
ggplot2::labs(
x = "residual lnM (SD)",
y = NULL,
title = "Scale -- Overall baseline",
caption = sprintf(
"residual SD = %.3f [%.3f, %.3f]",
sig_est, sig_lo, sig_hi
)
) +
theme_orchard()
p_comb <- patchwork::wrap_plots(p_loc, p_scl, ncol = 1)
save_plot(p_loc, "m00_orchard_location", width = 9, height = 4)
save_plot(p_scl, "m00_orchard_scale", width = 9, height = 4)
save_plot(p_comb, "m00_orchard_combined", width = 9, height = 8)
p_comb