3. Main Analyses Country Comparison

Load Data

Code
# load data for analysis
load("../data/wrangled_data/dt_ana_full.RData")

Analysis Strategy

We evaluate whether the predictive structure for violent intent differs across samples (countries; with a focal interest in the Jihadist sample) using two complementary strategies:

  • a multilevel model that separates within- from between-country effects and, in a second step, allows country-specific slopes;
  • a two-step meta-analysis that fits the same model within each country, then pools coefficient estimates and tests heterogeneity.

Why both?

The multilevel model performs one-step partial pooling (stabilizing noisy country estimates via shrinkage), while the meta-analysis gives two-step, fully transparent country-by-country estimates with classic heterogeneity statistics (Q, τ², I²). Agreement between the two strengthens credibility.

Prepare Data

We use group-mean centering (suffix _cw) for individual-level predictors to estimate within-country relationships, and include country means (grand-mean centered) to estimate between-country differences. This avoids conflating contextual differences with individual-level effects.

Code
ana_df <- dt_ana_full %>%
  mutate(
    across(
      .cols = matches("^(age|education_num|moral_neutralization|anger|negative_affect|past_activism|collective_relative_deprivation|perceived_discrimination|ingroup_superiority|obsessive_passion|commitment_passion|identity_fusion|violent_intent)_country_mean$"),
      .fns  = ~ .x - mean(.x, na.rm = TRUE),
      .names = "{.col}_c"
    )
  )

Multilevel Models

Random Intercept Model

This model allows country-specific baselines for violent intent while estimating common within-country slopes for the centered predictors and between-country effects for country means.

Code
mod_ml_ri <- lmer(
  violent_intent ~
    # CW (within-country) predictors
    age_cw +
    education_num_cw +
    moral_neutralization_cw +
    anger_cw +
    negative_affect_cw +
    past_activism_cw +
    collective_relative_deprivation_cw +
    perceived_discrimination_cw +
    ingroup_superiority_cw +
    obsessive_passion_cw +
    commitment_passion_cw +
    identity_fusion_cw +

    # Country-mean (grand-mean centered) predictors
    age_country_mean_c +
    education_num_country_mean_c +
    moral_neutralization_country_mean_c +
    anger_country_mean_c +
    negative_affect_country_mean_c +
    past_activism_country_mean_c +
    collective_relative_deprivation_country_mean_c +
    perceived_discrimination_country_mean_c +
    ingroup_superiority_country_mean_c +
    obsessive_passion_country_mean_c +
    commitment_passion_country_mean_c +
    identity_fusion_country_mean_c +

    # Contrast-coded / categorical covariates (as-is)
    noRel_vs_christian +
    noRel_vs_muslim +
    noRel_vs_buddhist +
    noRel_vs_jewish +
    noRel_vs_other +
    left_vs_right +
    left_vs_religious +
    female_vs_male +

    (1 | country),
  data = ana_df,
  REML = TRUE
)

summ(mod_ml_ri)
Observations 10395
Dependent variable violent_intent
Type Mixed effects linear regression
AIC 33905.21
BIC 34158.93
Pseudo-R² (fixed effects) 0.49
Pseudo-R² (total) 0.51
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 2.09 0.09 22.36 33.47 0.00
age_cw -0.01 0.00 -6.12 10341.52 0.00
education_num_cw -0.01 0.01 -1.16 10335.57 0.25
moral_neutralization_cw 0.45 0.01 30.36 10338.76 0.00
anger_cw 0.00 0.02 0.05 10335.35 0.96
negative_affect_cw 0.08 0.02 4.65 10336.69 0.00
past_activism_cw 0.20 0.04 4.98 10336.95 0.00
collective_relative_deprivation_cw 0.01 0.01 1.08 10336.08 0.28
perceived_discrimination_cw 0.25 0.01 19.18 10335.46 0.00
ingroup_superiority_cw 0.06 0.01 5.09 10336.01 0.00
obsessive_passion_cw 0.29 0.01 20.79 10337.95 0.00
commitment_passion_cw -0.09 0.01 -6.25 10336.54 0.00
identity_fusion_cw 0.02 0.01 1.87 10336.08 0.06
age_country_mean_c 0.01 0.01 1.20 27.02 0.24
education_num_country_mean_c 0.11 0.12 0.95 27.53 0.35
moral_neutralization_country_mean_c 0.41 0.31 1.32 26.74 0.20
anger_country_mean_c 0.62 0.29 2.16 27.60 0.04
negative_affect_country_mean_c 0.43 0.47 0.91 27.55 0.37
past_activism_country_mean_c 0.11 0.56 0.20 26.97 0.84
collective_relative_deprivation_country_mean_c 0.28 0.23 1.24 27.30 0.23
perceived_discrimination_country_mean_c -0.23 0.16 -1.47 27.59 0.15
ingroup_superiority_country_mean_c 0.16 0.19 0.82 26.58 0.42
obsessive_passion_country_mean_c 0.19 0.25 0.76 27.45 0.46
commitment_passion_country_mean_c -0.13 0.23 -0.56 28.20 0.58
identity_fusion_country_mean_c -0.11 0.28 -0.40 27.57 0.69
noRel_vs_christian 0.01 0.04 0.33 10332.77 0.74
noRel_vs_muslim 0.03 0.07 0.42 7001.63 0.67
noRel_vs_buddhist 0.07 0.12 0.57 1676.14 0.57
noRel_vs_jewish -0.17 0.13 -1.31 2203.18 0.19
noRel_vs_other 0.02 0.08 0.25 10361.99 0.80
left_vs_right -0.09 0.11 -0.85 26.77 0.41
left_vs_religious 0.11 0.19 0.59 27.32 0.56
female_vs_male 0.03 0.03 1.22 10361.81 0.22
p values calculated using Satterthwaite d.f.
Random Effects
Group Parameter Std. Dev.
country (Intercept) 0.24
Residual 1.22
Grouping Variables
Group # groups ICC
country 42 0.04
Code
# apa_lmer_summary(mod_ml_ri)

Random Slopes Model

We now allow country-specific slopes for the individual-level predictors. This addresses the question: do predictors of violent intent vary across samples? If so, we should see non-zero variance in the random slopes and improved model fit.

Code
mod_ml_rs <- lmer(
  violent_intent ~
    # CW predictors (exclude violent_intent_cw)
    age_cw +
    education_num_cw +
    moral_neutralization_cw +
    anger_cw +
    negative_affect_cw +
    past_activism_cw +
    collective_relative_deprivation_cw +
    perceived_discrimination_cw +
    ingroup_superiority_cw +
    obsessive_passion_cw +
    commitment_passion_cw +
    identity_fusion_cw +

    # Country-mean (grand-mean centered) predictors (exclude violent_intent_country_mean_c)
    age_country_mean_c +
    education_num_country_mean_c +
    moral_neutralization_country_mean_c +
    anger_country_mean_c +
    negative_affect_country_mean_c +
    past_activism_country_mean_c +
    collective_relative_deprivation_country_mean_c +
    perceived_discrimination_country_mean_c +
    ingroup_superiority_country_mean_c +
    obsessive_passion_country_mean_c +
    commitment_passion_country_mean_c +
    identity_fusion_country_mean_c +

    # Contrast-coded covariates
    noRel_vs_christian +
    noRel_vs_muslim +
    noRel_vs_buddhist +
    noRel_vs_jewish +
    noRel_vs_other +
    left_vs_right +
    left_vs_religious +
    female_vs_male +

    # Random intercepts + random slopes for *all CW variables*
    (1 +
       # age_cw +
       # education_num_cw +
       moral_neutralization_cw +
       # anger_cw +
       # negative_affect_cw +
       past_activism_cw +
       # collective_relative_deprivation_cw +
       perceived_discrimination_cw +
       ingroup_superiority_cw +
       obsessive_passion_cw +
       # commitment_passion_cw +
       identity_fusion_cw
     | country),
  data = ana_df,
  REML = TRUE,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl   = list(maxfun = 1e6)
  ),
  verbose = TRUE
)

start par. = 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 fn = 34286.4 At return eval: 3146 fn: 33429.017 par: 0.192815 -0.0175597 -0.0872950 0.0268182 0.00488176 0.0308624 0.0246887 0.130685 0.0549478 -0.0762193 -0.0188497 -0.0444966 0.00904905 0.117669 -0.0715363 0.0455067 -0.00808608 0.0166655 0.105864 -0.00744193 -0.0310848 -0.00145495 0.0516065 -0.0905898 0.0203626 0.0429689 0.00241279 5.48344e-07

boundary (singular) fit: see help('isSingular')
Code
summ(mod_ml_rs)
Observations 10395
Dependent variable violent_intent
Type Mixed effects linear regression
AIC 33553.02
BIC 34002.46
Pseudo-R² (fixed effects) 0.48
Pseudo-R² (total) 0.52
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 2.14 0.08 27.48 44.32 0.00
age_cw -0.01 0.00 -5.53 10096.20 0.00
education_num_cw -0.01 0.01 -1.35 9988.69 0.18
moral_neutralization_cw 0.42 0.03 14.80 48.28 0.00
anger_cw -0.01 0.02 -0.50 10273.20 0.62
negative_affect_cw 0.07 0.02 4.14 10214.58 0.00
past_activism_cw 0.19 0.05 3.83 37.18 0.00
collective_relative_deprivation_cw 0.03 0.01 2.24 10191.90 0.03
perceived_discrimination_cw 0.24 0.03 7.79 40.51 0.00
ingroup_superiority_cw 0.06 0.02 3.35 42.44 0.00
obsessive_passion_cw 0.30 0.03 11.50 44.73 0.00
commitment_passion_cw -0.09 0.01 -6.40 10151.85 0.00
identity_fusion_cw 0.02 0.01 1.37 94.61 0.17
age_country_mean_c 0.01 0.01 1.48 27.37 0.15
education_num_country_mean_c 0.16 0.09 1.73 28.84 0.09
moral_neutralization_country_mean_c 0.53 0.24 2.21 26.50 0.04
anger_country_mean_c 0.57 0.22 2.55 30.00 0.02
negative_affect_country_mean_c 0.31 0.36 0.86 29.03 0.40
past_activism_country_mean_c 0.47 0.42 1.13 23.54 0.27
collective_relative_deprivation_country_mean_c 0.04 0.17 0.25 27.52 0.81
perceived_discrimination_country_mean_c -0.21 0.12 -1.69 29.57 0.10
ingroup_superiority_country_mean_c 0.17 0.15 1.13 28.18 0.27
obsessive_passion_country_mean_c 0.07 0.19 0.38 29.20 0.71
commitment_passion_country_mean_c 0.03 0.18 0.16 30.80 0.87
identity_fusion_country_mean_c -0.10 0.21 -0.48 29.92 0.64
noRel_vs_christian -0.00 0.04 -0.01 9717.84 0.99
noRel_vs_muslim 0.01 0.07 0.14 4142.24 0.89
noRel_vs_buddhist -0.04 0.12 -0.34 806.03 0.73
noRel_vs_jewish -0.12 0.13 -0.99 656.43 0.32
noRel_vs_other -0.02 0.08 -0.21 10199.11 0.83
left_vs_right -0.07 0.09 -0.79 27.00 0.43
left_vs_religious 0.00 0.15 0.02 28.04 0.98
female_vs_male 0.03 0.03 1.35 10192.59 0.18
p values calculated using Satterthwaite d.f.
Random Effects
Group Parameter Std. Dev.
country (Intercept) 0.23
country moral_neutralization_cw 0.16
country past_activism_cw 0.19
country perceived_discrimination_cw 0.18
country ingroup_superiority_cw 0.09
country obsessive_passion_cw 0.14
country identity_fusion_cw 0.04
Residual 1.18
Grouping Variables
Group # groups ICC
country 42 0.04
Code
mod_ml_ri_ml <- update(mod_ml_ri, REML = FALSE)
mod_ml_rs_ml <- update(mod_ml_rs, REML = FALSE)
start par. =  1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 fn =  34219.91 
At return
eval: 3870 fn:      33278.917 par: 0.164438 -0.0237916 -0.0934695 0.0338811 0.00672257 0.0332343 0.0278346 0.127721 0.0511397 -0.0734431 -0.0179031 -0.0428649 0.0114227 0.112828 -0.0665871 0.0486860 -0.00378730 0.0147397 0.105149 -0.00718328 -0.0298592 -0.00330420 0.0473841 -0.0988220 0.0196426 8.22683e-06 -1.31934e-06 7.05612e-07
boundary (singular) fit: see help('isSingular')
Code
model_comp <- anova(mod_ml_ri_ml, mod_ml_rs_ml)  # LRT for added random slopes

df_model_comp <- as.data.frame(model_comp)

df_model_comp %>% 
  kable(., caption = "Comparing fixed and random slope model") %>% 
  kable_styling(full_width = F, latex_options = c("hold_position", "scale-down"))
Comparing fixed and random slope model
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
mod_ml_ri_ml 35 33757.23 34010.94 -16843.61 33687.23 NA NA NA
mod_ml_rs_ml 62 33402.92 33852.36 -16639.46 33278.92 408.3083 27 0

Model fit improves drasticially with random slopes (not suprisingly given the many different countries). We now assess how the Jihadist sample differs from the other samples.

Code
re_slopes <- broom.mixed::tidy(mod_ml_rs, effects = "ran_vals") %>%
  filter(grepl("_cw$", term)) %>% arrange(level, term) %>% 
  kable(., caption = "Random slopes of the multilevel model") %>% 
  kable_styling(full_width = F, latex_options = c("hold_position", "scale-down")) %>%
  scroll_box(width = "100%", height = "500px")

# re_slopes
Code
library(dplyr)
library(purrr)
library(tidyr)
library(broom)

country_models <- ana_df %>%
  group_by(country) %>%
  nest() %>%
  mutate(
    model = map(data, ~ lm(
      violent_intent ~
        age_gmz +
        education_num_gmz +
        moral_neutralization_gmz +
        anger_gmz +
        negative_affect_gmz +
        past_activism_gmz +
        collective_relative_deprivation_gmz +
        perceived_discrimination_gmz +
        ingroup_superiority_gmz +
        obsessive_passion_gmz +
        commitment_passion_gmz +
        identity_fusion_gmz,
      data = .x
    ))
  )

country_effects <- country_models %>%
  mutate(tidy = map(model, broom::tidy)) %>%
  select(country, tidy) %>%
  unnest(tidy) %>%
  filter(term != "(Intercept)") %>%
  select(country, term, estimate, std.error)

The following shrinkage (partial-pooling) contrasts, for each country and each predictor, the stand-alone country OLS (orinary least squares) slope with the multilevel Best Linear Unbiased Prediction (BLUP) slope which comes from the random-slopes model.

  • Hollow circle (OLS): the slope estimated by fitting a separate linear model within that country (no pooling).
  • Filled circle (BLUP): the country’s slope implied by the multilevel model = fixed effect + random deviation (i.e., after partial pooling).
  • Segment connecting them: the amount and direction of shrinkage from OLS → BLUP.
    • Long segments = noisier country estimates (often small n) get pulled more toward the overall mean slope.
    • Short segments = precise country estimates change little.

OLS

For each country, we fit a simple regression model just on that country’s data. The slope we get is the country-specific OLS slope.

Properties:

  • Uses only that country’s sample.
  • Very unbiased, but high variance if the country has a small sample size (noisy, unstable).

BLUB

The multilevel model’s estimate of a country’s slope = the global fixed slope plus that country’s random deviation.

Properties:

  • Shrinks noisy country slopes toward the overall mean slope (partial pooling).
  • Countries with large N or strong signal stay close to their OLS.
  • Countries with small N or noisy estimates get pulled toward the overall slope.
Code
# 3) Shrinkage plot: per-country OLS slopes vs. multilevel BLUPs (highlight Jihadist)
# Map *_gmz terms to *_cw to align labels
# 1) Recode gmz predictors to their CW counterparts (for easier matching)
ols <- country_effects %>%
  mutate(term = str_replace(term, "_gmz$", "_cw"))

# 2) Get fixed effects (all CW terms) from the random-slope model
fixef_df <- broom.mixed::tidy(mod_ml_rs, effects = "fixed") %>%
  filter(str_detect(term, "_cw$")) %>%
  select(term, fix = estimate)

# 3) Get random deviations (all CW terms) per country
re_df <- broom.mixed::tidy(mod_ml_rs, effects = "ran_vals") %>%
  filter(str_detect(term, "_cw$")) %>%
  select(country = level, term, ran = estimate)

# 4) Compute BLUPs = fixed slope + random deviation
blup <- re_df %>%
  left_join(fixef_df, by = "term") %>%
  mutate(blup = fix + ran)

# 5) Combine OLS and BLUPs into a single plotting dataset
plot_df <- ols %>%
  rename(ols = estimate) %>%
  select(country, term, ols) %>%
  inner_join(blup %>% select(country, term, blup), by = c("country","term")) %>%
  mutate(
    group    = if_else(country == "Jihadist", "Jihadist", "Other"),
    term_lab = term %>%
      str_remove("_cw$") %>%
      str_replace_all("_", " ") %>%
      str_to_title()  # nice readable labels
  )

# 6) Plot
# Helper: one predictor's shrinkage plot
shrinkage_plot_one <- function(df, term_label) {
  d <- df %>% dplyr::filter(term_lab == term_label)

  ggplot(d, aes(y = forcats::fct_reorder(country, blup))) +
    geom_segment(aes(x = ols, xend = blup, yend = country), alpha = 0.5, color = "grey60") +
    geom_point(aes(x = ols), shape = 1, size = 2) +
    geom_point(aes(x = blup, color = group), size = 2) +
    geom_vline(xintercept = 0, linetype = 2) +
    labs(
      x = "Slope",
      y = NULL,
      color = NULL,
      title = paste0("Partial pooling: OLS vs BLUP — ", term_label)
    ) +
    theme_minimal() +
    theme(
      axis.text.y  = element_text(size = 6),
      legend.position = "top"
    )
}

# Emit a single tabset with one tab per predictor
render_shrinkage_tabs <- function(df) {
  labs <- sort(unique(df$term_lab))

  # blank line before fenced div is IMPORTANT
  cat("\n\n::: {.panel-tabset}\n\n")

  for (lb in labs) {
    # each heading becomes a tab (no {.tabset} here)
    cat("\n")
    cat("#### ", lb, "\n\n", sep = "")
    print(shrinkage_plot_one(df, lb))
    cat("\n")
  }

  cat(":::\n")
}

# Render the tabs
render_shrinkage_tabs(plot_df)

Two Step Meta-Analysis

To compare effects across countries transparently, we refitted the same country-level linear model using grand-mean standardized predictors (suffix _gmz), extracted per-country coefficients and standard errors. We now meta-analyze them with a random-effects model (REML). This yields an overall average effect and heterogeneity statistics (Q, τ², I²).

Code
meta_results <- country_effects %>%
  tidyr::nest(data = c(country, estimate, std.error)) %>%  # one row per term
  dplyr::mutate(
    .results = purrr::map(data, ~ metafor::rma(
      yi  = .x$estimate,
      sei = .x$std.error,
      method = "REML",
      slab = .x$country
    ))
  ) %>%
  dplyr::select(term, .results)
Code
# 4) Forest + Baujat (meta-analytic echo), adapted to generic slopes/SEs

# Helper: draw ggplot-style forest 
make_forest_generic <- function(rma_obj,
                                main_title = "Meta-analysis",
                                order_by = NULL,
                                show_prediction_interval = TRUE,
                                study_labels = NULL,
                                highlight = "Jihadist",          # <- NEW: what to highlight (character pattern(s) or logical vector)
                                highlight_color = "red") {       # <- NEW: highlight color
  # robust SE: prefer $sei if present, otherwise sqrt(vi)
  sei_vec <- if (length(rma_obj$sei) > 0) as.numeric(rma_obj$sei) else sqrt(as.numeric(rma_obj$vi))
  wi <- 1 / (sei_vec^2)

  slab_vec <- if (!is.null(study_labels)) study_labels else rma_obj$slab

  # build highlight flag (case-insensitive substring match OR logical vector)
  hl_flag <- rep(FALSE, length(slab_vec))
  if (is.logical(highlight) && length(highlight) == length(slab_vec)) {
    hl_flag <- highlight
  } else if (is.character(highlight) && length(highlight) > 0) {
    for (pat in highlight) hl_flag <- hl_flag | grepl(pat, slab_vec, ignore.case = TRUE)
  }

  df <- tibble(
    slab = slab_vec,
    yi   = as.numeric(rma_obj$yi),
    sei  = sei_vec,
    ci_l = yi - qnorm(.975) * sei_vec,
    ci_u = yi + qnorm(.975) * sei_vec,
    w    = wi,
    hl   = ifelse(hl_flag, "highlight", "other")
  )

  if (!is.null(order_by) && order_by %in% names(df)) {
    df <- arrange(df, .data[[order_by]])
  }

  # overall pooled effect (+ CI + optional PI)
  overall <- tibble(
    slab = "Overall (REML)",
    yi   = as.numeric(rma_obj$b),
    ci_l = as.numeric(rma_obj$ci.lb),
    ci_u = as.numeric(rma_obj$ci.ub)
  )

  # prediction interval if available
  pi_band <- NULL
  if (isTRUE(show_prediction_interval)) {
    pr <- tryCatch(predict(rma_obj, transf = NULL), error = function(e) NULL)
    if (!is.null(pr) && !any(is.na(c(pr$pi.lb, pr$pi.ub)))) {
      pi_band <- c(pr$pi.lb, pr$pi.ub)
    }
  }

  g <- ggplot(bind_rows(df %>% select(-w, -sei), overall),
              aes(x = yi, y = reorder(slab, yi))) +
    geom_vline(xintercept = 0, linetype = 2, linewidth = 0.4) +
    annotate("rect",
             xmin = overall$ci_l[1], xmax = overall$ci_u[1],
             ymin = -Inf, ymax = Inf,
             alpha = 0.06) +
    {if (!is.null(pi_band))
      annotate("rect",
               xmin = pi_band[1], xmax = pi_band[2],
               ymin = -Inf, ymax = Inf,
               alpha = 0.04)} +
    # Error bars + points ONLY for study rows, colored by highlight flag
    geom_errorbarh(data = df,
                   aes(xmin = ci_l, xmax = ci_u, color = hl),
                   height = .15) +
    geom_point(data = df,
               aes(color = hl, size = pmin(w / max(w, na.rm = TRUE), 1)),
               show.legend = FALSE) +
    # pooled vertical + CI
    geom_vline(xintercept = overall$yi[1], linewidth = 0.4) +
    geom_vline(xintercept = overall$ci_l[1], linetype = "dashed", linewidth = 0.3) +
    geom_vline(xintercept = overall$ci_u[1], linetype = "dashed", linewidth = 0.3) +
    scale_color_manual(values = c(other = "black", highlight = highlight_color), guide = "none") +
    labs(x = "Effect (95% CI)",
         y = NULL,
         title = main_title,
         subtitle = paste0(
           "Pooled = ",
           formatC(overall$yi[1], digits = 3, format = "f"),
           " [", formatC(overall$ci_l[1], digits = 3, format = "f"),
           ", ", formatC(overall$ci_u[1], digits = 3, format = "f"), "]",
           if (!is.null(pi_band)) paste0(
             " | PI [",
             formatC(pi_band[1], digits = 3, format = "f"), ", ",
             formatC(pi_band[2], digits = 3, format = "f"), "]"
           ) else ""
         )) +
    theme_bw() +
    theme(panel.grid.minor = element_blank())

  return(g)
}

# Helper: draw funnel (generic) with highlighting
make_funnel_generic <- function(rma_obj,
                                main_title = "Funnel plot",
                                refline = 0,
                                levels = c(90, 95, 99),
                                grayscale = c("white", "gray55", "gray75"),
                                legend_pos = "topright",
                                highlight = "Jihadist",         # <- NEW
                                highlight_color = "red") {       # <- NEW
  # Ensure we have yi/vi/sei
  yi <- as.numeric(rma_obj$yi)
  vi <- if (!is.null(rma_obj$vi)) as.numeric(rma_obj$vi) else (as.numeric(rma_obj$sei))^2
  sei <- sqrt(vi)

  slab_vec <- rma_obj$slab
  # highlight flag (same rules as forest)
  hl_flag <- rep(FALSE, length(slab_vec))
  if (is.logical(highlight) && length(highlight) == length(slab_vec)) {
    hl_flag <- highlight
  } else if (is.character(highlight) && length(highlight) > 0) {
    for (pat in highlight) hl_flag <- hl_flag | grepl(pat, slab_vec, ignore.case = TRUE)
  }

  # Trim-and-fill on the fitted model
  tf <- tryCatch(trimfill(rma_obj), error = function(e) rma_obj)

  # PET / PEESE regressions
  pet   <- tryCatch(rma(yi = yi, sei = sei, mods = ~ sei), error = function(e) NULL)
  peese <- tryCatch(rma(yi = yi, sei = sei, mods = ~ I(sei^2)), error = function(e) NULL)

  reml_b <- as.numeric(rma_obj$b)

  funnel(tf, level = levels, shade = grayscale, refline = refline,
         legend = legend_pos,
         main = paste0(
           main_title, "\n",
           "REML = ", formatC(reml_b, digits = 2, format = "f"),
           if (!is.null(pet))   paste0(",  PET = ",   formatC(pet$b[1],   digits = 2, format = "f")) else "",
           if (!is.null(peese)) paste0(",  PEESE = ", formatC(peese$b[1], digits = 2, format = "f")) else "",
           ",  n = ", length(yi)
         ))

  # pooled marker
  abline(v = reml_b, lty = "dashed")
  points(x = reml_b, y = 0, cex = 1.4, pch = 17)

  # PEESE curve + point
  if (!is.null(peese) && length(peese$b) >= 2) {
    sei_grid <- seq(0, max(sei, na.rm = TRUE), length.out = 1000)
    vi_grid  <- sei_grid^2
    yi_grid  <- as.numeric(peese$b[1]) + as.numeric(peese$b[2]) * vi_grid
    lines(x = yi_grid, y = sei_grid)
    points(x = as.numeric(peese$b[1]), y = 0, cex = 1.3, pch = 5)
  }

  # PET line + point + CI segment
  if (!is.null(pet) && length(pet$b) >= 2) {
    abline(a = -as.numeric(pet$b[1]) / as.numeric(pet$b[2]),
           b =  1 / as.numeric(pet$b[2]))
    points(x = as.numeric(pet$b[1]), y = 0, cex = 1.2)
    segments(x0 = as.numeric(pet$ci.lb[1]), y0 = 0,
             x1 = as.numeric(pet$ci.ub[1]), y1 = 0,
             lty = "dashed")
  }

  # OVERDRAW highlighted studies in red so they pop
  if (any(hl_flag, na.rm = TRUE)) {
    points(x = yi[hl_flag], y = sei[hl_flag], pch = 19, col = highlight_color, cex = 1.2)
  }

  invisible(list(trimfill = tf, PET = pet, PEESE = peese))
}
Code
plot_meta_suite_tabbed <- function(meta_results,
                                   add_funnel      = TRUE,
                                   highlight       = "Jihadist",
                                   highlight_color = "red") {

  # Open ONE tabset container
  cat("::: {.panel-tabset}\n\n")

  terms <- unique(meta_results$term)

  for (tm in terms) {
    mobj <- meta_results %>%
      dplyr::filter(term == tm) %>%
      dplyr::pull(.results) %>%
      .[[1]]

    pretty <- tm %>%
      stringr::str_remove("_gmz$") %>%
      stringr::str_replace_all("_", " ") %>%
      stringr::str_to_title()

    # Each of these headings becomes a tab (NO {.tabset} here)
    cat("\n")
    cat("#### ", pretty, "\n\n", sep = "")

    # Forest (ggplot)
    print(make_forest_generic(
      mobj,
      main_title               = paste0("Forest: ", pretty),
      order_by                 = NULL,
      show_prediction_interval = TRUE,
      highlight                = highlight,
      highlight_color          = highlight_color
    ))

    if (isTRUE(add_funnel)) {
      # Make sure the base-graphics funnel is a separate figure
      graphics::plot.new()
      make_funnel_generic(
        mobj,
        main_title      = paste0("Funnel: ", pretty),
        highlight       = highlight,
        highlight_color = highlight_color
      )
    }

    cat("\n") # spacing between tabs
  }

  # Close the tabset container
  cat(":::\n")
}

# Call it
plot_meta_suite_tabbed(meta_results, add_funnel = TRUE)
Warning: `geom_errobarh()` was deprecated in ggplot2 4.0.0.
ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

`height` was translated to `width`.

References