orchaRd 2.0: An R package for visualising meta‐analyses with orchard plots

Although meta‐analysis has become an essential tool in ecology and evolution, reporting of meta‐analytic results can still be much improved. To aid this, we have introduced the orchard plot, which presents not only overall estimates and their confidence intervals, but also shows corresponding heterogeneity (as prediction intervals) and individual effect sizes. Here, we have added significant enhancements by integrating many new functionalities into orchaRd 2.0. This updated version allows the visualisation of heteroscedasticity (different variances across levels of a categorical moderator), marginal estimates (e.g. marginalising out effects other than the one visualised), conditional estimates (i.e. estimates of different groups conditioned upon specific values of a continuous variable) and visualisations of all types of interactions between two categorical/continuous moderators. orchaRd 2.0 has additional functions which calculate key statistics from multilevel meta‐analytic models such as I2 and R2. Importantly, orchaRd 2.0 contributes to better reporting by complying with PRISMA‐EcoEvo (preferred reporting items for systematic reviews and meta‐analyses in ecology and evolution). Taken together, orchaRd 2.0 can improve the presentation of meta‐analytic results and facilitate the exploration of previously neglected patterns. In addition, as a part of a literature survey, we found that graphical packages are rarely cited (~3%). We plea that researchers credit developers and maintainers of graphical packages, for example, by citations in a figure legend, acknowledging the use of relevant packages.


| INTRODUC TI ON
Meta-analysis has become an essential synthesis tool across the medical, social and biological sciences (Cooper et al., 2019;Gurevitch et al., 2018;Higgins et al., 2019;Schmid et al., 2021). In fields such as medicine, meta-analytic results are typically shown in a forest plot that presents effect sizes and their 95% confidence intervals (CIs) from each study in the meta-analysis. However, in ecology and evolution, forest plots are infrequently used because meta-analyses in this field often include >100 effect sizes, making a traditional forest plot impractical (Gurevitch et al., 2018;Senior et al., 2016). Instead, researchers use a 'forest-like plot' with the overall mean effect size estimate and their 95% CIs for different levels of a categorical moderator (predictor variable). Such estimates are derived from a metaregression model or from subset/sub-group analyses. For example, such a plot could show estimates from five different taxa, six different geographical areas or three different methods. A recent survey found 72 out of 102 ecological and evolutionary meta-analyses presented forest-like plots . Contributing to the popularity of forest-like plots is the fact that meta-analytic moderators are often categorical rather than continuous variables. Despite their popularity, forest-like plots in ecology and evolution often lack important information such as individual effect sizes and estimates of heterogeneity among effect sizes Schild & Voracek, 2015). Nakagawa et al. (2021) introduced an information-rich version of a forest-like plot, named the 'orchard' plot. Orchard plots provide (1) point estimates (i.e. meta-analytic means); (2) CIs, (3) prediction intervals (PIs; which show heterogeneity among effect sizes); and (4) individual effect sizes scaled by their precision (the inverse of the square root of the sampling variance). Nakagawa et al. (2021) implemented the orchard plot using functions that use the most popular and comprehensive meta-analysis R package, metafor (Viechtbauer, 2010) and ggplot2 graphics (Wickham, 2009). However, the original implementation (orchaRd 1.0) was limited to single moderator metaregression models. In addition, it was only possible to visualise a meta-regression model that assumed homoscedasticity across levels of the single categorical moderator (i.e. all levels have the same variance, which may be unrealistic, e.g. Wilson et al., 2022;Zajitschek et al., 2020).
In this article, we enhance the visualisation capabilities of the orchaRd package by integrating the functionalities of the R package emmeans (Lenth et al., 2018) in four ways. The first three extend orchard plots by allowing visualisation of (I) heteroscedasticity (different variances across levels of a categorical moderator; Section 3.1); (II) marginal estimates (e.g. marginalising all other moderators apart from the one visualised; Section 3.2) and (III) conditional estimates (i.e. estimates of different groups/levels of a categorical variable, conditioned upon specific values of a continuous variable; Section 3.3). The fourth capability allows for 'bubble' plots to be created of (i) a continuous variable, (ii) interactions between a continuous and categorical variable and (iii) interactions between two continuous variables from multi-moderator models (Section 3.4). In addition, we add helper functions to calculate key statistics from multilevel meta-analytic models such as I 2 (Cheung, 2014) and R 2 (Aloe et al., 2010;Nakagawa & Schielzeth, 2013), along with their CIs (Section 3.5). These new functionalities not only better visualise meta-analytic results in ecology and evolution, but also facilitate the exploration of previously neglected patterns, such as heteroscedasticity in meta-analytic data. Throughout we support the motivation for creating these functionalities with a survey of meta-analyses in ecology and evolution (Section 2).
Notably, orchaRd 2.0 improves reporting transparency in a metaanalysis by following the 'Preferred Reporting Items for Systematic reviews and Meta-Analyses in Ecology and Evolution' (PRISMA-EcoEvo;O'Dea et al., 2021). Our package's vignette also provides detailed instructions and examples on how to use all the main functions, as well as how to customise plots (https://danie l1nob le.github. io/orcha Rd/).

| SURVE Y ME THODS
To gauge the potential usefulness of the orchaRd package's extensions, we surveyed 102 meta-analyses in ecology and evolution.
Notably, this dataset was initially collected to quantify reporting quality of ecological and evolutionary meta-analyses to assist in  et al., 2021). We previously explored this dataset to survey the use of forest and forest-like plots in ecology and evolution .
For this study's survey, we asked the following 8 questions. Q6: How many papers that used a multi-moderator regression also modelled interactions? Q7: How many papers, which use R, cite an R software package they used for meta-analysis? Q8: How many papers, which use R, cite an R software package they used for the graphical presentation of meta-analytic results?
We report relevant results below, but the full results of this survey can be found in the Supporting Information.

| NE W SOF T WARE C APAB ILITIE S
The orchaRd 2.0 package has six main functions with three different  Table 1). Among these six functions, the core function is orchard_plot. This function enables users to draw orchard plots from a table created by mod_results, which uses emmeans functionality (Lenth et al., 2018) to process metafor model objects (object classes: rma, rma.mv and robust.rma; Viechtbauer, 2010). Below we first showcase three new capabilities of orchard_plot. Then, we describe a new function, bubble_plot, followed by the other main functions (caterpillars, i2_ml and r2_ml). Notably, the focus of our orchaRd package is to visualise multilevel meta-analytic models, which deal with two different types of non-independence due to (1) correlated effect sizes (e.g. multiple effect sizes per study) and (2) correlated sampling errors (e.g. shared control groups or shared measurements; see Noble et al., 2017). The former requires adding random effects (e.g. study ID), while the latter requires modelling a within-study variance-covariance matrix (note one can use the vcalc function in metafor to create such a matrix; Viechtbauer, 2010).

| Orchard plots: Heteroscedasticity
Categorical variables (moderators) are extremely common in metaanalyses. In our survey, >97% of the papers had at least one categorical variable. The categorical variable was used to subset data for sub-group analyses, where a series of meta-analyses (intercept models) were run, or to fit a meta-regression model (Q1). In many meta-analyses, researchers assumed all levels of a categorical moderator had the same variation (homoscedasticity). Our survey shows that only 5% of papers investigated heteroscedasticity, while others assumed homoscedasticity (Q2). Yet, differences in variances can be as biologically insightful as differences in means among groups. For example, Pottier et al. (2022) found that not only were aquatic ectotherms more thermally plastic than their terrestrial counterparts, but their plastic responses were much more variable than those of F I G U R E 1 Orchard plots (using orchard_plot function) and model result tables (using the mod_result function and extracting the table from the object-See Section 5 of the vignette) for terrestrial and aquatic ectotherm developmental acclimation response ratios (dARR). (a) Model assuming the variance in terrestrial and aquatic ectothermic species is the same (i.e. homogeneity of variance); (b) Model assuming the variance in terrestrial and aquatic ectothermic species is not the same (i.e. heterogeneity of variance), with the lower and upper confidence intervals (CIs) and prediction intervals (PIs) adjusted accordingly for each level of the habitat type moderator (data from Pottier et al., 2022). Siepielski et al., 2017; see also Noble et al., 2018). Absolute estimates can be calculated assuming a 'folded' normal distribution (see Morrissey, 2016;, with the accuracy of mean magnitudes being dependent on within-group variances. As such, it is important that heteroscedasticity is evaluated if such an approach is taken.

| Orchard plots: Marginal means
Many meta-analyses include multiple variables (moderators), and often they are modelled together in a single meta-regression model. In our survey, meta-analytic studies often modelled two or more categorical moderators together (Q3: 41%) and modelled at least one categorical moderator and one continuous modera- when a categorical moderator is unbalanced between groups/levels ( Figure 2). Equal weighting is, for example, useful when your sample is unequal in your dataset, but in the population, it should TA B L E 1 Main functions in the orchard package, their general categorisation and a description of what they can be used for in combination with metafor meta-analytic model objects (rma.mv, rma and robust.rma).

Function Category Description
mod_results Table  mod_results takes multi-level meta-analytic and meta-regression models (with multiple moderatorscontinuous or categorical) of class rma.mv/rma/robust.rma and calculates mean or marginalised mean meta-analytic estimates across all levels of a given moderator or overall (i.e. intercept only). The mod_ results table can then be used with orchard_plot, bubble_plot or caterpillars to plot results graphically. If a multivariate meta-regression model (with many moderators) is provided, users can specify the 'by' and/or 'at' arguments to marginalise over desired levels of other moderators orchard_plot Figure  Modified forest plot that plots the meta-analytic means, CIs, prediction intervals and raw data for each level of a categorical moderator. Users can use a number of arguments for modifying the look of plots including the legend, colour schemes, size and weight of points and lines and angle and naming of text on the axes. Sub-setting allows the users to plot a subset of the levels for a given moderator. Additional modifications can be made by adding and modifying layers of the ggplot object. Plots can be made using either mod_results objects directly or using the rma.mv/rma/robust.rma model object in combination with the raw data. If a multivariate meta-regression model (many moderators) is provided directly users can specify the 'by' and/or 'at' arguments to marginalise over desired levels of other moderators bubble_plot Figure Creates a bubble plot(s) depicting the predicted mean effect size, confidence and prediction interval as a function of a continuous moderator (slope estimate) or a series of separate plots showing predictions across an additional moderator (i.e. interaction plots). Plots can be made using either mod_results objects directly or using the rma.mv/rma/robust.rma model object in combination with the raw data. Raw data are plotted, and point size is adjusted according to effect size precision caterpillars Figure Creates a caterpillar plot from an intercept model or from mean effect size estimates for all levels of a given categorical moderator, their corresponding confidence and prediction intervals. Plots can be made using either mod_results objects directly or using the rma.mv/rma/robust.rma model object in combination with the raw data i2_ml Statistics Calculates heterogeneity statistics using measures of I 2 for a multilevel meta-analytic or meta-regression models. Point estimates can be calculated quickly for each level of random effect along with an estimate of total heterogeneity. Users also have the option of generating 95% CIs for all I 2 estimates using the 'boot' argument (percentile method). This argument will conduct parametric bootstrapping r2_ml Statistics Calculates marginal and conditional R 2 for multilevel meta-analytic or meta-regression models. Point estimates can be calculated quickly using a couple of different methods, but users also have the option of generating 95% CIs for R 2 using the 'boot' argument (percentile method). This argument will conduct parametric bootstrapping be ~50:50%; for example, males and females in many animals (cf. Deffner et al., 2022).

| Orchard plots: Conditional means
As mentioned above, our survey showed that it was not uncommon to have a study with a continuous moderator and a categori-

| Interactions: Orchard, bubbles and bubbleless
In our survey, ~30 (out of 102) meta-analyses modelled some type of interaction (Q5). Three types of interactions might manifest in a meta-analysis, those between (1)

| Other functions
In addition to orchard and bubble plots, the orchaRd package provides 'caterpillar' plots (via the function caterpillars, which is a forest plot without labels for each effect size; see our vignette-https://danie l1nob le.github.io/orcha Rd/). We also present two new non-plot functions to give meta-analysts convenient tools to quantify heterogeneity and variances explained by multilevel meta-analyses. The function i2_ml calculates I 2 , which is the percentage of variation among effect sizes not driven by sampling error (much of which is due to differences in sample sizes across studies; Higgins & Thompson, 2002).
Our function not only calculates the original I 2 (referred to as 'total' F I G U R E 2 Orchard plots of overall meta-analytic mean developmental acclimation response ratios (dARR). (a) Marginalised mean estimate assuming aquatic and terrestrial ectotherms are weighted proportionally to their representation in the sample of data (see Figure 1 for sample sizes for each group); (b) Marginalised mean estimate assuming aquatic and terrestrial ectotherms are weighted equally.
Comparing the mean and 95% confidence intervals shows how estimates affect the mean and the inferences. I 2 ) but heterogeneity explained by each additional random effect in the model (e.g. heterogeneity due to study ID or due to species ID; sensu Nakagawa & Santos, 2012). Furthermore, different sets of I 2 values can be calculated for different groups (levels) for a categorical moderator model with heteroscedasticity. While I 2 is estimated from a meta-analytic (intercept-only) model, R 2 is used to quantify variance (heterogeneity) accounted by moderators. The function r2_ml calculates marginal R 2 , proposed by Nakagawa and Schielzeth (2013) as a pseudo-R 2 for linear mixed-effects models. Notably, both i2_ml and r2_ml can provide 95% CIs, using bootstrapping.

F I G U R E 3
Orchard plots of mean log response ratios across four major trait categories (physiology, morphology, life history and behaviour) in fish. (a) Overall meta-analytic mean log response ratio for each trait category (marginalised means); (b) Predicted overall meta-analytic mean for each trait category for three levels of temperature difference, 5, 10 and 15° (conditional means; data from O'Dea et al., 2019). and 59%, respectively. As one can see, our package takes care of these three items in a single orchard plot (Figures 1-4). It is notable that now orchard plots even visualise different heterogeneities among different groups (i.e. heteroscedasticity) via PIs.

| Plea and proposal
Graphical presentation can facilitate better reporting in metaanalyses. However, in our survey, only 2 papers (3.1%) out of 64 articles which used R, cited any graphical package(s) used for visualising meta-analytic results (e.g. orchaRd; Q8). This figure starkly contrasts with 85% of the papers (55 out of 64; Q7) citing the software packages used for meta-analyses (e.g. metafor). This survey result marks a severe under-recognition of graphical packages. The real-world risk here is that this lack of recognition severely disincentivises developers from maintaining and further developing graphical packages.
We argue that authors should acknowledge graphical packages used for presenting meta-analyses (or any research article, for that matter), just as they do with any statistical package. We propose that graphical packages that were used to make a figure should be listed at the end of the figure legend. This standardised reporting format will mean packages do not necessarily need to be listed in the methods, but they will still be given credit. We note, however, that an R package can have many dependencies (i.e. other required R packages other than 'base' packages). For example, orchaRd 2.0 is dependent on emmeans, ggplot2 and metafor. We freely admit that we do not have a satisfying answer on whether dependencies should also be credited. However, for now we think it is reasonable to suggest that researchers provide the reference (in a figure legend or main text) for the immediate R function and package they used to make their figure.

| CON CLUS IONS
As the presence and influence of meta-analyses grow in the field, it is more important than ever to visualise meta-analytic results in an information-rich manner. Here, we have introduced an expanded version of orchaRd (version 2.0), which enables researchers to readily visualise complex as well as simple meta-analytic results, a task that was previously difficult for many. New functionalities that allow for marginal and conditional means to be plotted will improve model communication by allowing for a holistic visual interpretation of the complex numerical information generated by the analysis (see Figure 1-4).
Also, we introduce functions for calculating I 2 and R 2 for multilevel meta-analytic models, which have become standard in ecological and evolutionary meta-analyses. Finally, we hope our paper also becomes a reminder of the importance of acknowledging graphical packages.
Adequate attribution of credits will create a more sustainable environment for developers and maintainers of graphical packages. Rutkowska and Yefeng Yang conducted the survey. All authors contributed to the design of the study and to editing and commenting on drafts.

ACK N O WLE D G E M ENTS
We thank Wolfgang Viechtbauer for his continuous effort to maintain and develop an amazing package, metafor. We also acknowl-

CO N FLI C T O F I NTER E S T S TATEM ENT
The author reported no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-revie w/10.1111/ 2041-210X.14152.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data and code that are part of the R package can be found on GitHub (https://github.com/danie l1nob le/orchaRd) and Zenodo (https://doi.org/10.5281/zenodo.7928743; Nakagawa et al., 2023).