# Running intercept-only models for MHI and MABSHI using the maximum clade credibility (MCC) phylogeny.
# These models can be adapted based on the specifications outlined in Supplementary Table 1
# to reproduce all analyses described in the main text.


# Load packages
library(ape)
library(brms)
library(ggplot2)

# Load pre-processed data and tree
load("phylo_meta_analysis_input.RData")

# Fit MHI intercept-only model
model_intercept_mhi <- brm(
  mean_mhi | se(st_err_mhi) ~ 1 +
    (1 | gr(phylo, cov = Acov)) +
    (1 | species) +
    (1 | obs),
  data = data,
  family = gaussian(),
  data2 = list(Acov = mcc_vcv),
  prior = c(
    prior(normal(0, 1), class = "Intercept"),
    prior(cauchy(0, 0.05), class = "sd", group = "phylo", lb = 0),
    prior(cauchy(0, 0.05), class = "sd", group = "species", lb = 0),
    prior(cauchy(0, 0.05), class = "sd", group = "obs", lb = 0)
  ),
  chains = 4, cores = 4,
  warmup = 4000, iter = 12000,
  control = list(adapt_delta = 0.99, max_treedepth = 11)
)

# Summaries and checks
summary(model_intercept_mhi)
pp_check(model_intercept_mhi, ndraws = 100) + ggtitle("Posterior Predictive Check: MHI Intercept Model")
plot(model_intercept_mhi, ask = FALSE)

# Fit MABSHI intercept-only model
model_intercept_mabshi <- brm(
  mean_mabshi | se(st_err_mabshi) ~ 1 +
    (1 | gr(phylo, cov = Acov)) +
    (1 | species) +
    (1 | obs),
  data = data,
  family = gaussian(),
  data2 = list(Acov = mcc_vcv),
  prior = c(
    prior(normal(0, 1), class = "Intercept"),
    prior(cauchy(0, 0.05), class = "sd", group = "phylo", lb = 0),
    prior(cauchy(0, 0.05), class = "sd", group = "species", lb = 0),
    prior(cauchy(0, 0.05), class = "sd", group = "obs", lb = 0)
  ),
  chains = 4, cores = 4,
  warmup = 4000, iter = 12000,
  control = list(adapt_delta = 0.99, max_treedepth = 11)
)

# Summaries and checks
summary(model_intercept_mabshi)
posterior_summary(model_intercept_mabshi)
pp_check(model_intercept_mabshi) + ggtitle("Posterior Predictive Check: MABSHI Intercept Model")
plot(model_intercept_mabshi, ask = FALSE)