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 1 0 0 1 0 1 fn = 33971.76 At return eval: 372 fn: 33514.059 par: 0.196500 -0.0155714 0.0332990 0.0405932 0.131812 -0.0771914 -0.0501426 0.123680 -0.0635997 0.0553785

Code
summ(mod_ml_rs)
Observations 10395
Dependent variable violent_intent
Type Mixed effects linear regression
AIC 33602.06
BIC 33921.02
Pseudo-R² (fixed effects) 0.48
Pseudo-R² (total) 0.52
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 2.10 0.08 25.47 40.20 0.00
age_cw -0.01 0.00 -6.09 10209.18 0.00
education_num_cw -0.01 0.01 -1.37 10266.25 0.17
moral_neutralization_cw 0.42 0.03 14.68 48.30 0.00
anger_cw -0.01 0.02 -0.34 10301.42 0.73
negative_affect_cw 0.07 0.02 4.10 10280.44 0.00
past_activism_cw 0.20 0.04 4.76 9610.87 0.00
collective_relative_deprivation_cw 0.02 0.01 1.79 10306.65 0.07
perceived_discrimination_cw 0.24 0.03 7.86 41.22 0.00
ingroup_superiority_cw 0.06 0.01 5.67 10306.70 0.00
obsessive_passion_cw 0.29 0.02 12.23 60.94 0.00
commitment_passion_cw -0.09 0.01 -6.49 10264.75 0.00
identity_fusion_cw 0.02 0.01 1.77 10304.52 0.08
age_country_mean_c 0.01 0.01 1.15 27.25 0.26
education_num_country_mean_c 0.08 0.10 0.82 28.84 0.42
moral_neutralization_country_mean_c 0.46 0.26 1.75 26.94 0.09
anger_country_mean_c 0.71 0.24 2.91 29.13 0.01
negative_affect_country_mean_c 0.22 0.40 0.54 28.53 0.59
past_activism_country_mean_c 0.43 0.47 0.92 28.11 0.37
collective_relative_deprivation_country_mean_c 0.12 0.19 0.61 28.72 0.55
perceived_discrimination_country_mean_c -0.25 0.13 -1.90 29.98 0.07
ingroup_superiority_country_mean_c 0.11 0.16 0.65 28.05 0.52
obsessive_passion_country_mean_c 0.23 0.21 1.08 29.07 0.29
commitment_passion_country_mean_c 0.06 0.20 0.29 30.57 0.77
identity_fusion_country_mean_c -0.21 0.23 -0.91 29.19 0.37
noRel_vs_christian -0.00 0.04 -0.01 9768.22 0.99
noRel_vs_muslim 0.00 0.07 0.05 5463.20 0.96
noRel_vs_buddhist -0.01 0.12 -0.07 911.89 0.95
noRel_vs_jewish -0.13 0.13 -1.05 652.36 0.29
noRel_vs_other -0.01 0.08 -0.14 10245.59 0.89
left_vs_right -0.04 0.09 -0.39 26.76 0.70
left_vs_religious 0.08 0.16 0.48 28.14 0.63
female_vs_male 0.03 0.03 1.14 10314.83 0.25
p values calculated using Satterthwaite d.f.
Random Effects
Group Parameter Std. Dev.
country (Intercept) 0.23
country moral_neutralization_cw 0.16
country perceived_discrimination_cw 0.18
country obsessive_passion_cw 0.13
Residual 1.19
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 1 0 0 1 0 1 fn =  33891.4 
At return
eval: 469 fn:      33365.176 par: 0.163827 -0.0195790 0.0396625 0.0466030 0.129302 -0.0747822 -0.0474363 0.119988 -0.0662339 0.0462818
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 44 33453.18 33772.14 -16682.59 33365.18 322.0501 9 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 functions to harmonise axes across plots
source("../scripts/functions/meta-plot-helpers.R")

forest_limits <- meta_results$.results %>%
  purrr::map(extract_meta_axis_values, include_prediction_interval = TRUE) %>%
  purrr::flatten_dbl() %>%
  expand_axis_range(buffer = 0.05)

funnel_limits <- meta_results$.results %>%
  purrr::map(function(obj) {
    tf_obj <- tryCatch(metafor::trimfill(obj), error = function(e) obj)
    extract_meta_axis_values(
      obj,
      include_prediction_interval = TRUE,
      trimfill_obj = tf_obj
    )
  }) %>%
  purrr::flatten_dbl() %>%
  expand_axis_range(buffer = 0.05)

plot_meta_suite_tabbed(meta_results,
                       add_funnel = TRUE,
                       highlight = "Jihadist",
                       highlight_color = "red",
                       forest_limits = forest_limits,
                       funnel_limits = funnel_limits)
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