Extracts and tidies meta-analytic model-estimated effect-sizes and generates forest plots for *Eucalyptus* and blue tit datasets

generate_ManyAnalysts_forest_plots.R

# remotes::install_github("egouldo/ManyEcoEvo") # install latest ManyEcoEvo
library(tidyverse)
library(ManyEcoEvo)
library(metafor)
library(remotes)


# ------ define data extraction / plotting functions ------

get_forest_plot_data <- function(model){
  model %>% 
    broom.mixed::tidy(conf.int = TRUE, include_studies = TRUE) %>% 
    dplyr::mutate(
      point_shape = 
        ifelse(stringr::str_detect(term, "overall"), 
               "diamond", 
               "circle"),
      term = 
        forcats::fct_reorder(term, 
                             estimate) %>% 
        forcats::fct_reorder(., 
                             point_shape,
                             .desc = TRUE),
      parameter_type = case_when(str_detect(term, "overall") ~ "mean",
                                 TRUE ~ "study"),
      meta_analytic_mean =  pull(., estimate, type) %>% 
        keep_at(at = "summary")) %>% 
    select(-type, Parameter = term, everything())
}


plot_forest <- function(data, intercept = TRUE, MA_mean = TRUE){
  if (MA_mean == FALSE){
    data <- filter(data, Parameter != "overall")
  }
  
  p <- ggplot(data, aes(y = estimate, 
                        x =  Parameter, 
                        ymin = conf.low, 
                        ymax = conf.high,
                        shape = point_shape,
                        colour = parameter_type)) +
    geom_pointrange(fatten = 2) +
    ggforestplot::theme_forest() +
    theme(axis.line = element_line(linewidth = 0.10, colour = "black"),
          axis.line.y = element_blank(),
          text = element_text(family = "Helvetica")#,
          # axis.text.y = element_blank()
    ) +
    guides(shape = guide_legend(title = NULL), 
           colour = guide_legend(title = NULL)) +
    coord_flip() +
    ylab(bquote(Standardised~Effect~Size~Z[r])) +
    xlab(element_blank()) +
    # scale_y_continuous(breaks = c(-4,-3,-2,-1,0,1),
    # minor_breaks = seq(from = -4.5, to = 1.5, by = 0.5)) +
    NatParksPalettes::scale_color_natparks_d("Glacier")
  
  if(intercept == TRUE){
    p <- p + geom_hline(yintercept = 0)
  }
  if(MA_mean == TRUE){
    p <- p + geom_hline(aes(yintercept = meta_analytic_mean), 
                        data = data,
                        colour = "#01353D", 
                        linetype = "dashed")
  }
  
  return(p)
}

# ----- Load Data -----

data("ManyEcoEvo_viz")
ManyEcoEvo_viz %>% dim()
ManyEcoEvo_viz %>% head()

# ----- Extract Data and Plot -----

# dplyr::filter() arguments

filter_subsets <- rlang::exprs(
  exclusion_set == "complete",
  expertise_subset == "All",
  collinearity_subset == "All",
  publishable_subset == "All",
  estimate_type == "Zr",
  model_name == "MA_mod"
)

ManyEcoEvo_viz %>% 
  ungroup %>% 
  dplyr::filter(!!!filter_subsets) %>% 
  pull(model, dataset) %>% 
  map(get_forest_plot_data) %>% 
  map(plot_forest)