---
title: "Location–scale model grid"
---
::: {.workflow}
::: {}
**Input**
Effect sizes, sampling covariance, and phylogenetic matrix.
:::
::: {}
**Action**
Read saved model summaries or fit models in a separate run.
:::
::: {}
**Output**
Model status, coefficient, and diagnostic tables.
:::
::: {}
**Next**
Open the chapter for each fitted model.
:::
:::
```{r setup}
source(here::here("Scripts", "00_packages.R"))
source(here::here("Scripts", "01_paths.R"))
source(here::here("Scripts", "05_phylogeny.R"))
source(here::here("Scripts", "06_model_registry.R"))
source(here::here("Scripts", "07_model_formulas.R"))
source(here::here("Scripts", "08_fit_or_read_model.R"))
source(here::here("Scripts", "09_model_summaries.R"))
source(here::here("Scripts", "10_model_diagnostics.R"))
source(here::here("Scripts", "15_render_guard.R")) # read-only: never refit / rebuild phylogeny
# Change to TRUE only when you want to refit all models from scratch
refit_models <- FALSE
# The grid status + combined effect tables are pre-computed (precompute_grid_summaries.R)
# into Rdata/tables/grid_summary_cache.rds so the book never loads a 400 MB model.
# Set TRUE only to regenerate them by loading every model (slow).
rebuild_grid <- FALSE
```
## Overview
This chapter fits (or reads) one independent location–scale meta-regression for each moderator in the registry. No multi-moderator model is fitted. The registry is defined in `Scripts/06_model_registry.R`.
::: {.callout-important title="Render and fitting are separate tasks"}
**Normal book render:** keep `refit_models <- FALSE` and
`rebuild_grid <- FALSE`. The chapter reads saved results.
**Model-fitting run:** use the fitting scripts on a compute machine. After the
fits finish, run `precompute_grid_summaries.R`, then render the book. A normal
render must not start model fitting.
:::
```{r show-registry}
moderator_grid
```
## Input: load data and phylogeny
```{r load-data}
dat_es <- readRDS(here::here("Rdata", "effect_sizes", "proceed_lnm_safe.rds"))
cat("Total contrasts:", nrow(dat_es), "\n")
```
```{r load-phylogeny}
# build_or_read_phylogeny() returns list(A = ..., name_map = ...)
# It uses prepR4pcm for name normalisation + 4-stage reconciliation before rotl.
phylo_out <- tryCatch({
A_path <- here::here("Rdata", "phylogeny", "proceed_A_matrix.rds")
name_map_path <- here::here("Rdata", "phylogeny", "proceed_name_map.rds")
if (file.exists(A_path) && file.exists(name_map_path)) {
list(A = readRDS(A_path),
name_map = readRDS(name_map_path))
} else {
message("No cached phylogeny — building via prepR4pcm + rotl ...")
build_or_read_phylogeny(
species_vec = unique(as.character(dat_es$sp_ncbi)),
recompute = FALSE
)
}
}, error = function(e) {
message("Phylogeny load/build failed: ", conditionMessage(e))
NULL
})
A_full <- phylo_out$A
name_map <- phylo_out$name_map
if (!is.null(A_full)) {
cat("Phylogenetic matrix:", nrow(A_full), "×", ncol(A_full), "\n")
} else {
message("Phylogenetic matrix unavailable — models will omit the phylogenetic term.")
}
```
```{r apply-name-map}
# Apply the prepR4pcm name map to dat_es ONCE here, adding sp_ncbi_canonical.
# prepare_phylo_and_data() will use this column to match rows to A.
# Rows whose sp_ncbi had no canonical match get NA in sp_ncbi_canonical and
# are excluded from the complete primary-model dataset before fitting.
if (!is.null(name_map)) {
dat_es <- apply_phylo_name_map(dat_es, name_map)
n_resolved <- sum(!is.na(dat_es$sp_ncbi_canonical))
n_na_orig <- sum(is.na(dat_es$sp_ncbi) | dat_es$sp_ncbi == "")
n_unresolved <- sum(
!is.na(dat_es$sp_ncbi) & dat_es$sp_ncbi != "" &
is.na(dat_es$sp_ncbi_canonical)
)
cat("sp_ncbi_canonical coverage:\n")
cat(" Resolved to canonical name :", n_resolved, "\n")
cat(" Original sp_ncbi was NA/blank:", n_na_orig, "\n")
cat(" Named but unresolved :", n_unresolved, "\n")
} else {
# No name map — fall back to using sp_ncbi directly
dat_es$sp_ncbi_canonical <- dat_es$sp_ncbi
message("No name map available; using sp_ncbi as-is for phylogeny matching.")
}
```
## Verify moderators exist in the data
```{r check-moderators}
grid_status <- moderator_grid |>
dplyr::mutate(
in_data = moderator %in% names(dat_es),
n_nonmissing = purrr::map_int(moderator, function(m) {
if (m %in% names(dat_es)) sum(!is.na(dat_es[[m]])) else 0L
})
)
grid_status
```
The availability table is evaluated before phylogenetic name matching: all
7,237 eligible contrasts have both elapsed years and elapsed generations. Known generation time is
part of the inclusion filter. The fitted-model status
table below reports the analysis sample after matching to the phylogenetic
covariance matrix. Consequently, the elapsed-years and elapsed-generations
models share the maximum fitted sample size (7,186 contrasts); the 51 unmatched
contrasts are excluded from every phylogeny-adjusted primary model.
## Action: fit or read models
### Covariance matrices used by the model
Two group-level terms in the location submodel use covariance matrices supplied through `data2` in `brms`:
```r
(1 | gr(sp_ncbi, cov = A)) +
(1 | gr(es_id_model, cov = V))
```
`A` is the phylogenetic correlation matrix. It is calculated from the reconciled, pruned phylogeny and then subset (or expanded with zero off-diagonal correlations for unmatched species) so its row and column names match the `sp_ncbi` factor levels used by a given model:
```r
# Construct the full phylogenetic correlation matrix
A_full <- prepR4pcm::pr_phylo_cor(tree_final)
# Equivalent fallback used if pr_phylo_cor() fails:
# A_full <- ape::vcv.phylo(tree_final, corr = TRUE)
# Align A with the species retained for this moderator model
species_model <- levels(dat_model$sp_ncbi)
A <- A_full[species_model, species_model, drop = FALSE]
```
`V` is the known sampling-error covariance matrix. SAFE supplies one sampling variance per effect size (`vi_lnM_safe`), so the variances occupy the diagonal and the off-diagonal covariances are zero. Its dimension names must match the rebuilt `es_id_model` factor:
```r
dat_model <- dat_model |>
dplyr::mutate(es_id_model = factor(seq_len(dplyr::n())))
V <- diag(dat_model$vi_lnM_safe)
rownames(V) <- colnames(V) <- as.character(dat_model$es_id_model)
```
The aligned matrices are passed to `brms` as `data2 = list(A = A, V = V)`. The standard deviation of the `es_id_model` term is fixed to 1 with a `constant(1)` prior, ensuring that `V` contributes the known SAFE sampling uncertainty rather than an additional estimated variance component.
The loop below iterates over every row of the registry. For each moderator it:
1. Checks the moderator exists and has enough data.
2. Drops rows missing that moderator.
3. Rebuilds the sequential `es_id_model` and sampling-variance matrix `V`.
4. Subsets the phylogenetic matrix `A` to the species present.
5. Builds the formula and priors.
6. Calls `fit_or_read_model()`, which reads the cached `.rds` if it exists.
7. Extracts summaries and diagnostics.
```{r fit-grid, eval=rebuild_grid, results='hide'}
# Heavy path (rebuild_grid = TRUE only): load every model and extract summaries.
# Normal renders skip this and read the precomputed cache below.
all_fe <- list()
all_diag <- list()
status_rows <- list()
for (i in seq_len(nrow(moderator_grid))) {
row <- moderator_grid[i, ]
mod_id <- row$model_id
moderator <- row$moderator
mod_label <- row$label
mod_type <- row$type
model_file <- paste0(mod_id, "_ls_", moderator)
message("\n--- ", mod_id, ": ", moderator, " ---")
# ── 1. Check moderator exists ──────────────────────────────────────────────
if (!moderator %in% names(dat_es)) {
status_rows[[i]] <- tibble::tibble(
model_id = mod_id, moderator = moderator, label = mod_label,
type = mod_type, n_rows = 0L, n_levels = NA_integer_,
n_missing_excluded = NA_integer_,
model_path = NA_character_, status = "skipped: variable absent",
notes = "Moderator not found in dataset.",
max_rhat = NA_real_, min_bulk_ess = NA_real_, min_tail_ess = NA_real_,
n_divergent = NA_integer_, max_treedepth_hits = NA_integer_
)
next
}
# ── 2. Subset data for this moderator ─────────────────────────────────────
n_before <- nrow(dat_es)
dat_model <- dat_es[!is.na(dat_es[[moderator]]), ]
n_missing_ex <- n_before - nrow(dat_model)
if (mod_type == "categorical") {
dat_model[[moderator]] <- droplevels(factor(dat_model[[moderator]]))
}
# Check sample size per level for categorical moderators
n_levels <- NA_integer_
level_counts <- NULL
if (mod_type == "categorical") {
level_counts <- table(dat_model[[moderator]])
n_levels <- length(level_counts)
sparse <- level_counts[level_counts < 5]
if (length(sparse) > 0) {
message("Sparse levels (<5 obs) in ", moderator, ": ",
paste(names(sparse), collapse = ", "))
}
if (all(level_counts < 3)) {
status_rows[[i]] <- tibble::tibble(
model_id = mod_id, moderator = moderator, label = mod_label,
type = mod_type, n_rows = nrow(dat_model), n_levels = n_levels,
n_missing_excluded = n_missing_ex,
model_path = NA_character_, status = "skipped: too sparse",
notes = "All levels have fewer than 3 observations.",
max_rhat = NA_real_, min_bulk_ess = NA_real_, min_tail_ess = NA_real_,
n_divergent = NA_integer_, max_treedepth_hits = NA_integer_
)
next
}
}
if (nrow(dat_model) < 10) {
status_rows[[i]] <- tibble::tibble(
model_id = mod_id, moderator = moderator, label = mod_label,
type = mod_type, n_rows = nrow(dat_model), n_levels = n_levels,
n_missing_excluded = n_missing_ex,
model_path = NA_character_, status = "skipped: insufficient data",
notes = paste("Only", nrow(dat_model), "rows after missingness exclusion."),
max_rhat = NA_real_, min_bulk_ess = NA_real_, min_tail_ess = NA_real_,
n_divergent = NA_integer_, max_treedepth_hits = NA_integer_
)
next
}
# ── 3 & 4. Align data to phylogeny, build V ─────────────────────────────────
# prepare_phylo_and_data() subsets A to species present in this model's data,
# filters dat_model to the intersection (so every sp_ncbi level is in A),
# rebuilds es_id_model sequentially, and constructs the diagonal V matrix.
dat_model <- dat_model |>
dplyr::mutate(es_id_model = factor(seq_len(dplyr::n())))
phylo <- prepare_phylo_and_data(dat_model, A_full, label = mod_id)
dat_model <- phylo$dat_model
A_mod <- phylo$A_mod
V <- phylo$V
has_phylo <- phylo$has_phylo
# ── 5. Build formula and priors ───────────────────────────────────────────
formula <- build_ls_formula(moderator, has_phylogeny = has_phylo)
priors <- build_ls_priors(formula, dat_model, V, A = A_mod)
# ── 5b. Verify constant(1) is correctly set before any fitting ────────────
prior_ok <- verify_esid_prior(priors, model_name = mod_id)
if (!prior_ok) {
status_rows[[i]] <- tibble::tibble(
model_id = mod_id, moderator = moderator, label = mod_label,
type = mod_type, n_rows = nrow(dat_model), n_levels = n_levels,
n_missing_excluded = n_missing_ex,
model_path = NA_character_, status = "skipped: prior check failed",
notes = "constant(1) was not confirmed at coef=Intercept level for es_id_model.",
max_rhat = NA_real_, min_bulk_ess = NA_real_, min_tail_ess = NA_real_,
n_divergent = NA_integer_, max_treedepth_hits = NA_integer_
)
next
}
# ── 6. Fit or read ─────────────────────────────────────────────────────────
model_path_full <- file.path(
here::here("Rdata", "models"),
paste0(model_file, ".rds")
)
fit <- tryCatch({
fit_or_read_model(
model_name = model_file,
fit_fun = function() {
fit_ls_model(dat_model, formula, priors,
V = V, A = A_mod,
mcmc_args = default_mcmc_args)
},
model_dir = here::here("Rdata", "models"),
refit = refit_models
)
}, error = function(e) {
message("Model fitting failed for ", moderator, ": ", conditionMessage(e))
NULL
})
# ── 7. Extract summaries and diagnostics ──────────────────────────────────
if (is.null(fit)) {
status_rows[[i]] <- tibble::tibble(
model_id = mod_id, moderator = moderator, label = mod_label,
type = mod_type, n_rows = nrow(dat_model), n_levels = n_levels,
n_missing_excluded = n_missing_ex,
model_path = model_path_full, status = "failed",
notes = "Model returned NULL (fitting error).",
max_rhat = NA_real_, min_bulk_ess = NA_real_, min_tail_ess = NA_real_,
n_divergent = NA_integer_, max_treedepth_hits = NA_integer_
)
next
}
fe <- extract_fixed_effects(fit, mod_id, moderator)
diag <- tryCatch(
extract_diagnostics(fit, mod_id, moderator),
error = function(e) {
message("Diagnostics failed for ", moderator, ": ", conditionMessage(e))
tibble::tibble(
model_id = mod_id, moderator = moderator,
max_rhat = NA_real_, min_bulk_ess = NA_real_,
min_tail_ess = NA_real_, n_divergent = NA_integer_,
max_treedepth_hits = NA_integer_
)
}
)
all_fe[[mod_id]] <- split_location_scale(fe)
all_diag[[mod_id]] <- diag
status_rows[[i]] <- tibble::tibble(
model_id = mod_id,
moderator = moderator,
label = mod_label,
type = mod_type,
n_rows = nrow(dat_model),
n_levels = n_levels,
n_missing_excluded = n_missing_ex,
model_path = model_path_full,
status = "fitted",
notes = "",
max_rhat = diag$max_rhat,
min_bulk_ess = diag$min_bulk_ess,
min_tail_ess = diag$min_tail_ess,
n_divergent = diag$n_divergent,
max_treedepth_hits = diag$max_treedepth_hits
)
}
```
## Output: assemble grid summaries
```{r grid-assemble}
grid_cache <- here::here("Rdata", "tables", "grid_summary_cache.rds")
if (rebuild_grid) {
# Persist what the loop above produced (also refreshes the CSVs + cache).
model_status <- dplyr::bind_rows(status_rows)
readr::write_csv(model_status, here::here("Rdata", "tables", "model_grid_status.csv"))
summ <- save_combined_summaries(all_fe, tables_dir = here::here("Rdata", "tables"))
diag_tbl <- save_combined_diagnostics(all_diag, tables_dir = here::here("Rdata", "tables"))
loc_effects <- summ$location; scl_effects <- summ$scale
saveRDS(list(status = model_status, location = loc_effects,
scale = scl_effects, diagnostics = diag_tbl), grid_cache)
} else {
# Normal render: read the small precomputed cache — no 400 MB model is loaded.
pace_need(grid_cache, "model grid summary cache")
g <- readRDS(grid_cache)
model_status <- g$status; loc_effects <- g$location
scl_effects <- g$scale; diag_tbl <- g$diagnostics
}
model_status <- model_status |>
dplyr::filter(model_id %in% moderator_grid$model_id) |>
dplyr::arrange(model_id)
loc_effects <- loc_effects |>
dplyr::filter(model_id %in% moderator_grid$model_id)
scl_effects <- scl_effects |>
dplyr::filter(model_id %in% moderator_grid$model_id)
```
## Model status table
This table contains only the six primary moderator models.
```{r status-table}
model_status |>
dplyr::select(
model_id, label, n_rows, n_levels, n_missing_excluded, status,
max_rhat, min_bulk_ess, min_tail_ess, n_divergent,
max_treedepth_hits
)
```
## Convergence overview
Convergence criteria: max $\hat{R} \leq 1.01$, bulk and tail ESS $\geq 400$,
and no divergent transitions.
```{r convergence}
model_status |>
dplyr::filter(status == "fitted") |>
dplyr::select(model_id, moderator, max_rhat, min_bulk_ess, min_tail_ess,
n_divergent) |>
dplyr::mutate(
rhat_ok = dplyr::if_else(max_rhat <= 1.01, "ok", "warn"),
ess_ok = dplyr::if_else(min_bulk_ess >= 400 & min_tail_ess >= 400, "ok", "warn"),
div_ok = dplyr::if_else(n_divergent == 0, "ok", "warn")
)
```
## Excluded and dropped models
These models were fitted but are not part of the reported results. Models
m06, m08, m09, and m10 did not meet the convergence criteria; m11
(measurement scale) was dropped.
```{r non-primary-models}
excluded_grid |>
dplyr::select(model_id, label, reason) |>
knitr::kable(col.names = c("Model", "Field", "Recorded status"))
```