Extract effect from brms model

Helpers.R

extract_region_effects <- function(fit, region_var = "ROI", parameter = "Intercept",
                                   probs = c(0.025, 0.05, 0.5, 0.95, 0.975),
                                   digits = 8, scale = 1) {
  # Check if the fit object is from brms
  if (!inherits(fit, "brmsfit")) {
    stop("The 'fit' argument must be a brmsfit object")
  }
  
  # Extract posterior samples for population-level effects
  aa <- fixef(fit, summary = FALSE) / scale
  
  # Extract group-level (random) effects
  bb <- ranef(fit, summary = FALSE)
  
  # Check if the specified parameter exists
  if (!parameter %in% colnames(aa)) {
    stop(paste0("Parameter '", parameter, "' not found in the model"))
  }
  
  # Check if the specified region variable exists
  if (!region_var %in% names(bb)) {
    stop(paste0("Region variable '", region_var, "' not found in the model random effects"))
  }
  
  # Calculate number of posterior samples
  ns <- nrow(aa)
  
  # Calculate posterior samples at each region for the parameter
  ps <- apply(bb[[region_var]][, , parameter], 2, "+", aa[, parameter])
  
  # Create initial result data frame
  result <- data.frame(
    mean = apply(ps, 2, mean),
    SD = apply(ps, 2, sd)
  )
  
  # Add P+ column (probability of effect being positive)
  # Using quote() to preserve the "+" in the column name
  p_plus <- apply(ps, 2, function(x) sum(x > 0) / ns)
  result <- cbind(result, "P+" = p_plus)
  
  # Add quantiles
  quantiles <- t(apply(ps, 2, quantile, probs = probs))
  colnames(quantiles) <- paste0(100 * probs, "%")
  
  # Combine results
  result <- cbind(result, quantiles)
  
  # Round values to specified number of digits
  result <- round(result, digits)
  result <- result %>%  tibble::rownames_to_column(var = "label")
  # Return the summarized results
  return(result)
}