orchaRd_vignette.qmd

---
title: "**orchaRd 2.0: An R package for drawing 'orchard' plots (and 'caterpillars' plots) from meta-analyses and meta-regressions**"
author: "**S. Nakagawa, M. Lagisz, R. E. O’Dea, P. Pottier, J. Rutkowska, Y. Yang, A. M. Senior & D.W.A. Noble**"
date: "`r Sys.Date()`"
format:
  html:
    toc: true
    toc-location: left
    toc-depth: 3
    toc-title: "**Table of Contents**"
    output-file: "index.html"
    theme: simplex
    embed-resources: true
    code-fold: show
    code-tools: true
    number-sections: true
    bibliography: ./bib/ref.bib
    fontsize: "12"
    max-width: "10"
    code-overflow: wrap
crossref: 
  fig-title: Figure     # (default is "Figure")
  tbl-title: Table     # (default is "Table")
  title-delim: —     # (default is ":")
  fig-prefix: Fig.   # (default is "Figure")
  tbl-prefix: Tab.    # (default is "Table")
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = FALSE, tidy = TRUE)
options(digits=2)
```


# **Introduction**

`orchaRd` 2.0 allows users to create pretty orchard plots that contain both confidence intervals (CIs) and prediction intervals (PIs) around point estimates for different categorical moderators, plots the effect size distribution over top such estimates and weights effect sizes by their precision (1/standard error, SE) or sample size. `orchaRd` takes a `metafor` object of class `rma.mv` or `rma` [@Viechtbauer2010] and plots the results for the meta-analytic or meta-regression model. 

`orchaRd` 2.0 now allows users to take meta-regression models with a multitude of moderator variables that are categorical or continuous along with heterogeneous residual variance models. This is achieved by allowing the user to produce marginalised mean estimates using the `emmeans` package. For plotting, `orchaRd` uses `ggplot2` [@Wickham2009], and as such, layers can be added directly to make plots customizable to the users needs. Also, the function `orchard_plot` uses the principles of bee swarm plots to visualise individual effect sizes [@Eklund2012; see also @Clarke2017; @Sherrill-Mix2017].

In addition to orchard plots, we have another plotting function `caterpillars` to draw caterpillar(s) plots. Furthermore, the package `orchaRd` comes with two functions that are useful for multilevel meta-analysis and meta-regression (i.e., `rma.mv`): 1) `i2_ml` to calculate the 'total' I^2^ along with I^2^ for each level [*sensu* @NakagawaSantos2012; see also @Senior2016] and 2) `r2_ml` to calculate marginal and conditional R^2^ for mixed models [*sensu* @Nakagawa2013].

# **Citing orchaRd**

To cite `orchaRd` 2.0 in publications one can use the following reference:

Shinichi Nakagawa, Malgorzata Lagisz, Rose E. O'Dea, Patrice Pottier, Joanna Rutkowska, Alistair M. Senior, Yefeng Yang, Daniel W.A. Noble. 2023. orchaRd 2.0: An R package for visualizing meta-analyses with orchard plots. **Methods in Ecology and Evolution**, https://doi.org/10.1111/2041-210X.14152 (preprint = **EcoEvoRxiv**, https://doi.org/10.32942/X2QC7).  

For earlier versions please cite:

Nakagawa, S., Lagisz, M., O'Dea, R. E., Rutkowska, J., Yang, Y., Noble, D. W., & Senior, A. M. (2021). The Orchard Plot: Cultivating Forest Plots for Use in Ecology, Evolution and Beyond. **Research Synthesis Methods**, 12: 4-12, https://doi.org/10.1002/jrsm.1424 (preprint = *EcoEvoRxiv* https://doi.org/10.32942/osf.io/epqa7)

# **Installation**

To install `orchaRd` use the following code in R:


:::{.column-margin}
:::{.callout-note}
If you want to download the development branch, which contains newer features you can execute the following code:

```{r}
#| echo: true
#| eval: false
#| results: hide
#| code-fold: true
rm(list = ls())
devtools::install_github("daniel1noble/orchaRd", ref = "dev", force = TRUE)
pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, orchaRd, emmeans, ape, phytools, flextable)
```

It's important to note that this version may not be as stable as the master branch.

:::
:::

```{r echo = TRUE, eval = TRUE, results = 'hide', warning=FALSE, message=FALSE,  class.source='klippy'}
#install.packages("pacman")
rm(list = ls())
devtools::install_github("daniel1noble/orchaRd", ref = "main", force = TRUE)
pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, orchaRd, emmeans, ape, phytools, flextable)
```

Installation will make the primary functions accessible to users along with their help files. You will also need the `tidyverse` and `metafor` packages.

# **Examples of orchard and caterpillars plots**

In this vignette we'll walk the reader through a number of case studies and show you how you can create beautiful looking orchard plots. We overview three different case studies that make use of different effect sizes and moderators. The data sets associated with each case study come as part of the `orchaRd` package.

## **Example 1: Dietary Restriction and Lifespan**

@EnglishUller2016 performed a systematic review and meta-analysis on the effects of early life dietary restriction (a reduction in a major component of the diet without malnutrition; e.g. caloric restriction) on average age at death, using the standardised mean difference (often called *d*). They found that, across the whole data set, there was little evidence for an effect of dietary restriction on mean age at death. Here we'll use the data set to first calculate the effect size measures and then fit an intercept only, meta-analytic model.

```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}

data(english)

# We need to calculate the effect sizes, in this case d
english <- escalc(measure = "SMD", n1i = NStartControl, sd1i = SD_C, m1i = MeanC, n2i = NStartExpt, sd2i = SD_E, m2i = MeanE, 
                  var.names=c("SMD","vSMD"),
                  data = english)

english_MA <- rma.mv(yi = SMD, V = vSMD, random = list( ~ 1 | StudyNo, ~ 1 | EffectID), data = english)
summary(english_MA)
```

We have fit a meta-analytic model and thus the only estimate we see is the overall effect size on the effects of caloric restriction on mean death across all studies examined. Now that we have fit our meta-analytic model we can get the confidence intervals and prediction intervals with a few functions in the `orchaRd` package. If one is interested in getting the table of results we can use the `mod_results` function. This will allow users to make nice tables of the results if needed. We can do that as follows:

:::{.column-margin}
:::{.callout-important}
In `orchaRd` version 2.1 and above, the data argument has now been removed from the `mod_results` and `orchard_plot`, `bubble_plot` and `catepillars` functions. This is because the data is already contained within the `metafor` object.
:::
:::

```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}

model_results <- orchaRd::mod_results(english_MA, mod = "1", at = NULL,  group = "StudyNo")
model_results
 
```

If we instead want to create an orchard plot and visualise the results we can do so quite simply as:

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure1
#| fig-cap: "Orchard plot of the impact caloric restriction using standardised mean difference"
orchaRd::orchard_plot(english_MA, mod="1", group = "StudyNo", xlab = "Standardised mean difference", transfm = "none", twig.size = 0.5, trunk.size = 1) 
```

In @fig-Figure1, we simply add in the metafor model and it will create a default orchard plot. Alternatively, we could also add in the table of results. Also. some people may wish to present the heterogeneity index, I^2^ and, at the same time, change the colours of the 'trunk and fruit' by adding a few arguments.

```{r , echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure2
#| fig-cap: "Orchard plot of the impact caloric restriction using standardised mean difference but instead of using the metafor model, using the model results table"
#| fig-width: 6
#| fig-height: 4

# use i2_sn function to obtain the total I^2

I2 <- orchaRd::i2_ml(english_MA)

orchaRd::orchard_plot(model_results, mod="1", xlab = "Standardised mean difference") + 
  annotate(geom="text", x= 0.80, y= 1, label= paste0("italic(I)^{2} == ", round(I2[1],4), "*\"%\""), 
              color="black", parse = TRUE, size = 5) +
  scale_fill_manual(values="grey") +
  scale_colour_manual(values="grey")
```

Figure @fig-Figure2 and @fig-Figure1 above show that overall estimate from a random-effects meta-analysis of 77 effect sizes is centered on zero, with a 95% CI that spans the line of no-effect. The prediction intervals clearly demonstrate the high level of heterogeneity (over 75%), with effects size less than -0.5 and greater than 0.5 predicted to be observed. 

We can now draw a comparable caterpillar plot by using the function `caterpillars`.

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure3
#| fig-cap: "Catterpilar plot of the impact caloric restriction using standardised mean difference"

# a caterpillar plot (not a caterpillars plot)
orchaRd::caterpillars(model_results, mod="1", xlab = "Standardised mean difference") 
```

It is now interesting to compare how the same data set can be visualised different by look at @fig-Figure2 and @fig-Figure3. It is good that the caterpillar plot can show CIs for each effect size. An equivalent function was done with the size of fruit in the orchard plot. 

In a subsequent publication, @Senior2017 analysed this data set for effects of dietary-restriction on among-individual variation in the age at death using the log coefficient of variation ratio. A major prediction in both @EnglishUller2016 and @Senior2017 was that the type of manipulation, whether the study manipulated quality of food versus the quantity of food would be important. As such, we can fit a meta-regression model to test whether the moderator "Manipulation Type" `ManipType` impacts our results on the mean and variance

```{r, warning=FALSE, message=FALSE, fig.width = 6, fig.height = 8}
#| label: fig-Figure4
#| fig-cap: "Orchard plot of diet qualities impact on SMD (top) and log coefficient of variation (bottom)"

# First we need to calculate the lnCVR
english <- escalc(measure = "CVR", n1i = NStartControl, sd1i = SD_C, m1i = MeanC, 
                        n2i = NStartExpt, sd2i = SD_E, m2i = MeanE,
                        var.names=c("lnCVR","vlnCVR"),
                        data = english)

# Now we can fit the meta-regression model (contrast)
english_MR0 <- rma.mv(yi = SMD, V = vSMD, mods = ~ ManipType, 
                     random = list(~ 1 | StudyNo, ~ 1 | EffectID), data = english)
summary(english_MR0)

# Again, we can create a table of results
res2 <- orchaRd::mod_results(english_MR0, mod = "ManipType", group = "StudyNo")
res2

# we fit a meta-regression (contrast) for the log coefficient of variation
senior_MR0 <- rma.mv(
  yi = lnCVR, V = vlnCVR, mods = ~ ManipType, random = list(~ 1 | StudyNo, ~ 1 | EffectID),
  data = english
)
summary(senior_MR0)


# creating a table of results
res3 <- orchaRd::mod_results(senior_MR0, mod = "ManipType",  group = "StudyNo", at = list(ManipType = "Quality"), subset = TRUE)
res3

# We can now plot SMD and lnCVR beside each other and compare the results
p1 <- orchaRd::orchard_plot(res2, 
                   mod = "ManipType", group = "StudyNo", xlab = "Standardised mean difference")

p2 <- orchaRd::orchard_plot(res3,
  mod = "ManipType", group = "StudyNo", xlab = "log(CV ratio) (lnCVR)")

(p1 / p2)
```

Our orchard plot for the log coefficient of variation demonstrates that, while restrictions on dietary quality and quantity do not affect the average age at death (top of @fig-Figure4), among-individual variation may be altered by quality restrictions (bottom of @fig-Figure4). The effect is negative suggesting that the coefficient of variation in the control group is lower than that in the treatment group, and the 95% CI does not span zero. Again though, the effect is heterogeneous; a substantial number of positive effects are still predicted. 

Let's draw comparable caterpillars plots to compare these two results again. 

```{r, echo = TRUE, eval = TRUE,  warning=FALSE, message=FALSE}
#| label: fig-Figure5
#| fig-cap: "Catterpilar plot of diet qualities impact on the log coefficient of variation (lnCVR)"

# We can now plot SMD and lnCVR beside each other and compare the results
p1c <- orchaRd::caterpillars(res2, 
                   mod = "ManipType", group = "StudyNo",  xlab = "Standardised mean difference")

p2c <- orchaRd::caterpillars(res3,
  mod = "ManipType",  group = "StudyNo", xlab = "log(CV ratio) (lnCVR)")

p1c / p2c
```

Now compare between @fig-Figure4 and @fig-Figure5. We think they both look nice. 

## **Example 2: Predation and Invertebrate Community**

@Eklof2012 evaluated the effects of predation on benthic invertebrate communities. Using the log response ratio they quantified differences in abundance and/or biomass of gastropods and Amphipods in groups with and without predation in an experimental setting. 

Here again, we can create orchard plots of the model results, but we'll show how a few simple things can be modified. Again, we can fit the meta-analytic model first: 

```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}

data(eklof)

# Calculate the effect size
eklof <- escalc(
  measure = "ROM", n1i = N_control, sd1i = SD_control, m1i = mean_control,
  n2i = N_treatment, sd2i = SD_treatment, m2i = mean_treatment,
  var.names = c("lnRR", "vlnRR"),
  data = eklof
)

# Add the observation level factor
eklof$Datapoint <- as.factor(seq(1, dim(eklof)[1], 1))

# Also, we can get the sample size, which we can use for weighting if we would like
eklof$N <- rowSums(eklof[, c("N_control", "N_treatment")])

# fit a meta-regression with the intercept (contrast)
eklof_MR0 <- rma.mv(
  yi = lnRR, V = vlnRR, mods = ~ Grazer.type, random = list(~ 1 | ExptID, ~ 1 | Datapoint),
  data = eklof
)

summary(eklof_MR0)


```

Above we have fit a meta-regression model using "Grazer Type" as a moderator which is predicted to explain variation in log response ratios. We can demonstrate a few simple changes users can make, but we note here that users can make far more complex changes down the line if needed, but we'll save explaining these more complex changes for the last example. The first is the angle at which the y-axis labels are positioned (bottom of @fig-Figure6):

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure6
#| fig-cap: "Orchard plots of the effects of predation on benthic invertebrate communities compared using the log response ratio. Top panel is the default plot and bottom panel a plot containing changes in label axes"

#| 
f <- orchaRd::mod_results(eklof_MR0, mod = "Grazer.type",  group = "ExptID") 

p3 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type",  group = "ExptID", xlab = "log(Response ratio) (lnRR)", 
  transfm = "none") 

p4 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type",  group = "ExptID", xlab = "log(Response ratio) (lnRR)", 
  transfm = "none", angle = 45, g = FALSE) 

p3/p4
```

The other thing we can change is the type of scaling we wish to use. Let's say we are interested in scaling the effect size by the total sample size of the study we use a vector of N, sample size (bottom of @fig-Figure7:

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure7
#| fig-cap: "Orchard plots of the effects of predation on benthic invertebrate communities compared using the log response ratio. Top panel is the default plot and bottom panel a plot containing changes in label axes and scaling with sample size instead of precision"

p5 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type",  group = "ExptID", xlab = "log(Response ratio) (lnRR)") 

p6 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type",  group = "ExptID", xlab = "log(Response ratio) (lnRR)", angle = 45, N = "N", g = FALSE) 

p5/p6

```

Overall, our orchard plot shows the results of a re-analysis of their data. The estimated mean effects are negative for both gastropods and amphipods suggesting that mean abundance/biomass in the control group is lower than in the treatment groups, although the effect is largest, and is statistically significant, for amphipods. In both cases the prediction intervals reveal the extent of heterogeneity, with positive effects predicted to be observed.

Again, we can also draw the caterpillars-plot counterpart for this data set (bottom of @fig-Figure8):

```{r , echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure8
#| fig-cap: "Caterpillars plot of the effects of predation on benthic invertebrate communities compared using the log response ratio"

eklof_MR_result <- orchaRd::mod_results(eklof_MR0, mod = "Grazer.type", group = "ExptID")
orchaRd::caterpillars(eklof_MR_result, mod = "Grazer.type", group = "ExptID", xlab = "log(Response ratio) (lnRR)") 

```
Like last time, we compare between @fig-Figure7 and @fig-Figure8. Again, they both look quite informative. 

The last thing we can do is transform these data to percentage. You can also make other transformation depending on what effect size that you are working with, see `?orchard_plot`.

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure7perc
#| fig-cap: "Orchard plots of the effects of predation on benthic invertebrate communities using the percentage difference."

p6_per <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type",  group = "ExptID", xlab = "Percentage Change (%)", angle = 45, N = "N", g = FALSE, transfm = "percentr") 

p6_per

```

## **Example 3: Maternal-Offspring Morphological Correlations**

Finally, we also look at the case discussed by @Lim2014, who meta-analysed the strength of correlation between maternal and offspring size within-species, across a very wide range of taxa. They found, that typically, there is a moderate positive correlation between maternal size and offspring size within species (i.e. larger mothers have larger offspring). However, they also found evidence for relatively strong phylogenetic effects suggesting that the strength of the association is dependent on evolutionary lineage.

Here we have used an orchard plot to represent the results obtained when meta-analysing the data from @Lim2014 by taxonomic Phylum. 

```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
	data(lim)

# Add in the sampling variance
	lim$vi<-(1/sqrt(lim$N - 3))^2

# Lets fit a meta-regression - I will do Article non-independence. The phylogenetic model found phylogenetic effects, however, instead we could fit Phylum as a fixed effect and explore them with an Orchard Plot
  lim_MR<-metafor::rma.mv(yi=yi, V=vi, mods=~Phylum-1, random=list(~1|Article, ~1|Datapoint), data=lim)
  summary(lim_MR)
  
```

Noe we can plot a default orchard plot, scaling each effect size by N. Also, because we are using Zr, we can use transfm = "tanh" and it will do the conversions for us:

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure9
#| fig-cap: "Orchard plot of the the strength of correlation between maternal and offspring size within-species"

  #Plot the meta-regression model
  orchaRd::orchard_plot(lim_MR, mod = "Phylum", group = "Article",  xlab = "Correlation coefficient", 
    alpha = 0.5, transfm = "tanh", angle = 45, N = "N", cb = FALSE) 
```

Now that we have @fig-Figure9, we can do some small changes to make it pretty and add some more details, such as between study heterogeneity and R^2^ for phylum. Currently, the cb argument is "FALSE", we can change this to "TRUE" to use colour blind friendly colours. Additionally, because we are using `ggplot` we can add element to the figure to make it look pretty. 

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure10
#| fig-cap: "Orchard plot of the the strength of correlation between maternal and offspring size within-species"

  # Lets add the R2 on the figure as well as little modifications:
  R2 <- orchaRd::r2_ml(lim_MR)

  # We can draw an orchard plot with R2
  orchaRd::orchard_plot(lim_MR, mod = "Phylum", group = "Article",  xlab = "Correlation coefficient (r)", 
                        alpha = 0.5, transfm = "tanh", angle = 45, N = "N", cb = TRUE) +
theme(legend.position= c(0.05, 0.99), legend.justification = c(0,1), legend.key.size = unit(1, "mm")) +
theme(legend.direction="horizontal", legend.title = element_text(size =8),legend.text = element_text(size = 10)) + 
annotate(geom="text", x=0.8, y=0.7, label=paste0("italic(R)^{2} == ", round(R2[1],4)*100), color="black", parse = TRUE, size = 4)

```

As in @fig-Figure10, new elements can be added to the `orchard_plot` to modify it as one sees fit. It will overwrite existing elements. From our orchard plots above, it is clear that the analysis is dominated by data from Chordates and Arthropods, with the other Phyla being much more poorly represented. Second, there is a difference between the strength of a typical correlation within these two well represented groups (the correlation is stronger in Chordates), which arguably would explain the phylogenetic signals detected by @Lim2014. Lastly, although there are differences within the typical correlation between Chordates and Arthropods, there remains a large overlap in predicted range of individual effect sizes; individual species within Phyla are still highly variable.

Finally, we make a comparable caterpillars plot.
```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure11
#| fig-cap: "Caterpillars plot of the the strength of correlation between maternal and offspring size within-species"

  #Plot the meta-regression model
  lim_MR_results <- orchaRd::mod_results(lim_MR, mod = "Phylum", group = "Article")
  orchaRd::caterpillars(lim_MR_results, mod = "Phylum", group = "Article",  xlab = "Correlation coefficient", transfm = "tanh", g = FALSE)
```

Compared to the orchard plot @fig-Figure10, @fig-Figure11 does not look as elegant, unlike the previous comparisons. This is because the groups with small sample sizes are packed at the top of the figure. 

From @fig-Figure10 we can see that many taxonimic groups have very little data. We may want instead to just focus the readers attention to means of groups with sufficient data. The new `mod_results` function allows us to do this using the maginalised means approach. We can do that as follows:

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure121
#| fig-cap: "Caterpillars plot of the the strength of correlation between maternal and offspring size within-species"
 
  #Plot the meta-regression model
  lim_MR_results <- orchaRd::mod_results(lim_MR, mod = "Phylum", group = "Article", 
                                           at = list(Phylum = c("Chordata", "Arthropoda", "Mollusca")),
                                           subset = TRUE)
  orchaRd::orchard_plot(lim_MR_results, xlab = "Correlation coefficient", transfm = "tanh", g = TRUE, angle = 45)
```

We can now see that @fig-Figure121 limits what it plots to only the three levels relevant for presentation. 

We may also at times want to revert back to a normal plot with catagorical variables on the x-axis. We can do this by setting `flip` argument to `FALSE`. 

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-Figure1212
#| fig-cap: "Caterpillars plot of the the strength of correlation between maternal and offspring size within-species"
 
  #Plot the meta-regression model
  lim_MR_results <- orchaRd::mod_results(lim_MR, mod = "Phylum", group = "Article", 
                                           at = list(Phylum = c("Chordata", "Arthropoda", "Mollusca")),
                                           subset = TRUE)
  orchaRd::orchard_plot(lim_MR_results, xlab = "Correlation coefficient", transfm = "tanh", g = TRUE, flip = FALSE)
```

# **Orchard Plots with Meta-regression Models: Marginalised and conditional means with heteroscedastic residual variances**

In `orchaRd` 2.0 we've made some major updates to the types of models that can be used with the package. Meta-analyses in ecology and evolution often fit a multitude of moderators and random effects within a single model. This is even desirable in many cases to deal with high heterogeneity and disentangle the independent contributions of moderators to changes in mean effect size. The challenge with the previous version of the package was that it was difficult to accommodate these models. However, the `emmeans` package is an excellent solution to the problem. It will produce marginsalised or conditional mean estimates for a given level of moderator. This includes models with categorical and continuous moderators. 

An additional issue that is of importance with meta-regression models is the need to deal with heteroscedastic variance across levels of a moderator. Homogeneity of variance assumptions are often violated in meta-analysis. `orchaRd` 2.0 allows the user to fit models that explicitly estimate different residual variances reducing the probability of type I errors. 

In this section, we show users: 1) how to fit these models using `metafor` and 2) how these models interface with the `orchaRd` package to produce pretty plots that capture the model predictions. 

## **Example 4: Meta-regression and heterogeneous variance models** {#sec-sec4}

Meta-regression models with single categorical moderators are extremely common. Meta-analyst's want to estimate the mean effect size within levels of a moderator and test whether these means are significantly different from each other. One inherent assumption to fitting these models is that we are assuming that the within study (or residual) variance within these groups is the same (i.e., "homogeneity of variance assumption"). This assumption will likely be violated in many situations in meta-regression models. Unless one conducts a sub-group analysis for each level of the moderator (although this is less powerful and makes it hard to test whether levels are significantly different) different residual variances need to be explicitly modeled in `metafor`. Fortunately, this can be done. `orchaRd` 2.0 now allows plotting of these models with their corresponding CI's and PI's.

We'll work with a meta-analysis by @ODea2019. In this meta-analysis, the authors were interested in testing whether developmental temperature in fish had an effect on mean and variance in a host of phenotypic traits. They used experimental studies manipulating developmental temperatures using the log response ratio (`lnrr`) and the log coefficient of variation (`lncvr`) as their effect sizes. We'll only focus on `lnrr` here for demonstration purposes. Importantly, they categorised effect sizes as having come from one of 4 trait categories: Behaviour, Life-history, Morphology and Physiology. They were interested in whether developmental temperature impact trait types differently. We'll first explore how to fit a meta-regression model that explicitly estimates residual variance in each of these trait levels. 

First, load the data. 
```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}

  # Load the dataset that comes with orchaRd
    data(fish)

  # Subset data for demonstration purposes.
    warm_dat <- fish
```

To fit a heterogeneous variance model in `metafor` we need to use some additional arguments and specify the random effect structure differently (i.e., specify inner argument). To do this, we need to add the `struc = HCS` to the argument list, which specifies the variance structure of an ~inner|outer formula of a particular random effect. In this case, we're assuming a heteroscedastic compound symmetry structure. To avoid estimating the correlation / covariance between levels we need to also add the `rho` argument and set this to 0. This will assume that the covariance between the levels is set to 0 (although it could be estimated). 

Once we have these additional arguments we need to modify our random effect structure. Note that `metafor` does not estimate a residual / within study variance. We need to explicitly add an observation level random effect, in this case `es_ID`. If we have more than one effect size / study (which is common) then we want to have this random effect in our models. 

```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}

model_het <- metafor::rma.mv(    yi = lnrr, V = lnrr_vi, 
                               mods = ~ trait.type, method = "REML", test = "t", 
                             random = list(~1 | group_ID, ~1 + trait.type| es_ID), 
                                rho = 0, struc = "HCS", 
                               data = warm_dat, 
                            control = list(optimizer="optim", optmethod="Nelder-Mead"))
model_het

```

We can see from above that we now have four variances; one estimated for each level of the `trait.type` moderator. Now that we have fit our model and estimated the relevant parameters we can produce the orchard plot  using the `marginal_means` function:

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-figure12
#| fig-cap: "Orchard plot with heteroscedastic error for trait typess"


het_model <- orchaRd::mod_results(model_het, 
                                     group = "group_ID", 
                                       mod = "trait.type")
  
orchaRd::orchard_plot(het_model, xlab = "lnRR")
```

We can now see from @fig-figure12 that the variance in effects for the morphological and physiological trait categories are smaller than for life-history and behavioural trait categories; likely because of the very small sample sizes in the latter trait categories.

As a second example, we can turn to @Pottier2021 who did a meta-analysis comparing acclimation differences between males and females. @Pottier2021 found that aquatic and terrestrial systems differed in their variability of effects. Here, we'll also show how phylogenetic correlation matrices can be incorporated within the models and work with orchaRd.

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
# Clear working space
rm(list = ls())

# Load the data
data("pottier")

# Load the phylogenetic correlation matrix
data("phylo_matrix")

# Load the sampling variance matrix
data("VCV_dARR")
```

Now that we have all the data and matrices loaded, we can fit the multi-level mete-regression model. 

```{r, pot_model, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
##### Heteroscedasticity modeled at the effect size level
mod.habitat_het <-  rma.mv(yi=dARR, V=VCV_dARR,mods= ~habitat,method="REML",test="t",dfs="contain",
                           random=list(~1|species_ID,~1|phylogeny,~habitat|es_ID),struct="HCS", rho=0, R =list(phylogeny = phylo_matrix),data=pottier)
mod.habitat_het
```

Now we can produce an `orchaRd` plot using this model to show the unequal variance in the two habitat types. 

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-figsure13a
#| fig-cap: "Orchard plot with heteroscedastic error for habitat type"

orchard_plot(mod.habitat_het,  group = "species_ID", mod = "habitat", xlab = "dARR", angle = 45)

```

We can also create the model results table (@tbl-tab1). 

:::{.column-body}
```{r}
#| label: tbl-tab1
#| tbl-cap: Model results table for a model with heteroscedastic error for habitat type
#| 
flextable::flextable(mod_results(mod.habitat_het,  group = "species_ID", mod = "habitat")$mod_table)
```
:::

<!-- We can check calculation for, say,aquatic  is correct: 0.2087 + qt(0.975, 136)*sqrt(sum(0.0655^2 + 0.0716 + 0.0129 + 0.0090)) -->


We might also think about heterogeneous variance existing at multiple random effect levels. These models require a lot of data, and need to be carefully thought about, but we `orchard` 2.0 can also handle these models. 

```{r, pot_model2, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
##### Heteroscedasticity modeled at the effect size level
mod.habitat_het2 <-  rma.mv(yi=dARR, V=VCV_dARR,mods= ~habitat,method="REML",test="t",dfs="contain",
                           random=list(~habitat|species_ID,~1|phylogeny,~habitat|es_ID),struct="HCS", rho=0, phi = 0, R =list(phylogeny = phylo_matrix),data=pottier)
mod.habitat_het2

```

Again, we can plot this new model normally.

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-figsure14a
#| fig-cap: "Orchard plot with heteroscedastic error for habitat type at the residual and species level"
orchard_plot(mod.habitat_het2, group = "species_ID", mod = "habitat", xlab = "dARR", angle = 45)

```

We can also view the model results table (@tbl-2)

:::{.column-body}
```{r}
#| label: tbl-2
#| tbl-cap: Model results table for a model with heteroscedastic error for habitat type at the residual and species level
flextable::flextable(mod_results(mod.habitat_het2, group = "species_ID", mod = "habitat")$mod_table)
```
:::

## **Example 5: Bubble plots with continuous moderators**

`orchaRd` 2.0 can now also create bubble plots. Bubble plots are extremely useful when one has a continuous moderator variable. For example, this might be publication year or sampling variance if one is testing for various forms of publication bias. We'll demonstrate using two examples.

First, we'll return to the @Lim2014 dataset. Here, we can fit a model that assessing time-lag effects for studies on captive and wild populations. Say, for the sake of argument, that we expect the time-lag effect to vary by lab/wild conditions. We can fit an interaction model and create bubble plots for these situations (@fig-bubbleslim):

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-bubbleslim
#| fig-cap: "Bubble plot with @Lim2014 dataset assessing publication bias (time-lag bias) in wild and captive studies."

data(lim)
lim[, "year"] <- as.numeric(lim$year)
lim$vi<- 1/(lim$N - 3)
model<-metafor::rma.mv(yi=yi, V=vi, mods= ~Environment*year,
                       random=list(~1|Article,~1|Datapoint), data=na.omit(lim))


lim_bubble <- orchaRd::mod_results(model, mod = "year", group = "Article",
                    weights = "prop", by = "Environment")

orchaRd::bubble_plot(lim_bubble, group = "Article",  mod = "year", xlab = "Year", legend.pos = "top.left")
```

As second example, we'll return back to @Pottier2021, who collected data on body mass for species in the data set where the acclimation response ratio (ARR) was calculated. We might expect the ARR to vary by body mass because small species might be able to acclimate faster than larger species. Using a simplified version of @Pottier2021's model, we can look at this effect (@fig-bubblespatrice):

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-bubblespatrice
#| fig-cap: "Bubble plot with @Pottier2021 dataset assessing whether body mass affects the acclimation response ratio (ARR)."
#| 
  # extract complete cases
      data(pottier)
      pottier2 <- pottier[complete.cases(pottier$body_mass),]

  # Model body mass effects
mod.body_mass <-  metafor::rma.mv(yi = dARR, V = Var_dARR,mods= ~body_mass, method="REML", test="t", dfs="contain",
                           random=list(~1|species_ID, ~1|es_ID), data = pottier2)

orchaRd::bubble_plot(mod.body_mass, mod = "body_mass", group = "species_ID", xlab = "Body mass (g)", ylab = "Acclimation Response Ratio (ARR)", legend.pos = "top.left")

```

We may want to colour code the points based on some grouping variable, such as habitat. As indicated through the vignette, we can override layers of the orchard plot to do this fairly easily (@fig-bubblespatrice22).

```{r}
#| label: fig-bubblespatrice22
#| fig-cap: "Bubble plot with @Pottier2021 dataset assessing whether body mass affects the acclimation response ratio (ARR), with points coloured by habitat type. 'Blue' are terrestrial species, 'red' are aquatic species."
#| message: FALSE
#| warning: FALSE
orchaRd::bubble_plot(mod.body_mass, mod = "body_mass", group = "species_ID", xlab = "Body mass (g)", ylab = "Acclimation Response Ratio (ARR)", legend.pos = "top.left") + geom_point(data = pottier, aes(x = body_mass, y = dARR, color = habitat, size = 1/sqrt(Var_dARR)), alpha = 0.6) + scale_color_discrete()
```

<br>

## **Example 6: Meta-regression and conditional mean effect sizes**

We'll now demonstrate how multi-moderator meta-regression models can be used with `orchaRd` 2.0, by plotting conditional mean effects within categories of a single moderator. 

To demonstrate how to do this we'll again turn to @ODea2019 as our example data set. We'll first load the data set, if you haven't already done so:

```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}

  # Load the dataset that comes with orchaRd
    data(fish)

  # Subset data for demonstration purposes.
    warm_dat <- fish
```

Now that we have the data set assume we are interested in modelling multiple moderators that we hypothesize will explain effect size variance. We want to provide a model that will allow us to understand the effects of each multiple moderators on mean effect size. We know that `experimental_design` and the temperature difference between treatments (`deg_diff`) are expected to impact mean effect size. We also expect experiments that last longer to have stronger effects (`treat_end_days`). We'll now extend our model from @sec-sec4 above to include these moderators so we can look at the effects of each controlling for each other. 

```{r echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}

 # Fit the metaregerssion model
model <- metafor::rma.mv(yi = lnrr, V = lnrr_vi, 
                         mods = ~ experimental_design + trait.type + deg_dif + treat_end_days,
                         method = "REML", test = "t", 
                         random = list(~1 | group_ID, ~1 + trait.type| es_ID), 
                                rho = 0, struc = "HCS", 
                         data = warm_dat, 
                         control=list(optimizer="optim", optmethod="Nelder-Mead"))

model
```

Our model is complicated. It not only estimates the effects of trait type, but type of experimental design, degree difference and treatment duration. But, effects of moderators are now conditional on each other. What happens if we just want to plot mean effects for trait type, but we want to condition these effects on a given temperature different? We can do this using calculating conditional meta-analytic means. We can see that there's not much going on for experimental design, so we will marginalize the means across each level of this moderator. Degree difference is continuous. To deal with continuous moderators we need to specify the levels we want predictions to be made. 

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-figure13
#| fig-cap: "Orchard plot of meta-anlaytic mean estimates conditioned on a 5, 10, and 15 degree Celcius temperture difference"
HetModel <- orchaRd::mod_results(model, group = "group_ID", mod = "trait.type", at = list(deg_dif = c(5, 10, 15)), by = "deg_dif", weights = "prop" )

orchaRd::orchard_plot(HetModel, xlab = "lnRR" , angle = 45, g = FALSE, legend.pos = "top.left", condition.lab = "Temperature Difference") + 
  theme(legend.direction = "vertical")

```

We can see the orchard plot changes. Here we have a 'condition' label, which we have renamed to "Temperature Difference" using the `condition.lab` argument. This new legend specifies the degree difference for each mean predicted. Note that these means effectively average over all levels of experimental design **and** all possible values of `treat_end_days` (i.e., marginalize). We have also made some small changes to the plot. First, we set `g = FALSE` which will suppress the total groups on the plot (e.g. studies), moved the position of legends and angled the labels. 

We can also specify the specific value of `treat_end_day`, further conditioning our mean estimates. We can do that as follows:

```{r, echo = TRUE, eval = TRUE, warning=FALSE, message=FALSE}
#| label: fig-figure14
#| fig-cap: "Orchard plot of meta-anlaytic mean estimates conditioned on a 5, 10, and 15 degree Celcius temperture difference"
HetModel2 <- orchaRd::mod_results(model, group = "group_ID", mod = "trait.type", at = list(deg_dif = c(5, 10, 15), treat_end_days = c(10)), by = "deg_dif", weights = "prop")

orchaRd::orchard_plot(HetModel2, xlab = "lnRR", angle = 45, g = FALSE, legend.pos = "top.left", condition.lab = "Temperature Difference") + 
  theme(legend.direction = "vertical")


```

Now the orchard plot depicts the mean effect sizes (plus 95% CI and PI) for each trait category, at 5, 10 and 15 $^{\circ}$C when the treatment duration is 10 days. 


\newpage
# **References**