effect-size.qmd

# Calculating and interpreting effect sizes

```{r}
#| label: setup-packages-size
#| echo: false
#| message: false
#| warning: false

# Load packages
pacman::p_load(metafor, flextable, tidyverse, orchaRd, pander, mathjaxr, equatags)
options(digits=3)
```

## **Introduction**

Meta-analysis quantitatively aggregates effect size data collected from existing research. The effect size one uses for a study question needs to be chosen carefully. Each make certain assumptions about the underlying data. It's important to be weary of these assumptions and identify when something might have gone wrong. 

The effect size chosen must satisfy two important features. First, the effect size must be comparable across studies. This means that they must be placed on the same scale so that they can be aggregated. Some effect size estimates can be inter-converted to each other (e.g., Zr, Hedges' g), but this isn't the case for all of them. Second, we must be able to derive or approximate the sampling variance for a given effect size. As [already indicated](introduction-to-meta.qmd), the sampling variance is very important if one wishes to use meta-analytic models that account for sampling variance and is needed to report certain measures of [heterogeneity](het.qmd). 

Meta-analyses in comparative physiology, and indeed, ecology and evolution more generally, often use highly heterogeneous data [@Noble2022; @Senior2016]. This poses other challenges in aggregating effect size data that make their comparability questionable. We'll discuss some of these issues in later [tutorials](nuis-het.qmd) and identify some ways in which some of this heterogeneity can be dealt with when deriving the effect size itself. 

By the end of this tutorial, we hope that you will feel comfortable:

1. Calculating different effect sizes from your data.
2. Interpreting the meaning of different effect size statistics
3. Transforming effect size estimates back to their original scales when needed

## **Common Effect Sizes**

What is an effect size? The answer to this question depends on what kind of study one wishes to synthesise. For example, an effect size could be a contrast between experimental treatments. We might extract mean differences from experiments manipulating, say, bisphenol-A (BPA) and looking at the effect of BPA on aquatic ectotherm phenotype relative to some control [@WuSeebacher2020]. Alternatively, the effect of interest might be the magnitude and direction of a correlation between two observed variables, such as metabolism and behaviour [@Holtman2016]. At times, we might even just be interested in analysing the mean or variance itself. 

> "An effect size is a statistical parameter that can be used to compare, on the same scale, the results of different studies in which a common effect of interest has been measured" [@Rosenberg2013]

<br>
There are many types of effect sizes that are commonplace in ecological and evolutionary research (@tbl-tablea describes many of the more common ones from @Noble2022). 

```{r}
#| label: tbl-tablea
#| echo: false
#| tbl-cap: "**Common effect sizes used throughout meta-analyses in comparative physiology, their associated sampling variances and examples on when they might be used**. M: mean; SD: standard deviation; N: sample size; CTmin: critical thermal minimum; CTmax: critical thermal maximum; LC_50: median lethal concentration; LT_50: median lethal temperature; BPA: bisphenol A."

FitFlextableToPage <- function(ft, pgwidth = 6){

  ft_out <- ft %>% flextable::autofit()

  ft_out <- flextable::width(ft_out, width = dim(ft_out)$widths*pgwidth /(flextable::flextable_dim(ft_out)$widths))
  return(ft_out)
} 

# I'm just listing a bunch here, but we should think about what we want to include. Here's an example of how we can make a table with all equations. Fun! Really useful website here: https://math.meta.stackexchange.com/questions/5020/mathjax-basic-tutorial-and-quick-reference. Something wrong with character coding when rendering word document. Has to do with the log odds ratio equation. Renders fine in html, but just stops when word
names <- c("Mean",
           "Log Standard deviation, lnSD",
           "Log Response Ratio, lnRR",
           # leave this for now??? - we introduce this later
           "Standardised Mean Difference, SMD",
           "Zr (Fisher Transformation of \n Correlation Coefficient, r)")
           #"Log Odds Ratio",
           #"")

# the index i is not really applied consistently so I am adding it more (to do)
# we shoudl discuss what kinds of notations we want to use consistely
           
eqns_yi <- c("M",
             "\\ln{SD} + \\frac{1}{2(N - 1)}",
             # Senior will have a better version but this is good for now
             "\\ln \\left( \\frac{M_{1}}{M_{2}}\\right)",
             "\\frac{\\left( M_{2} - M_{1}\\right)}{SD_{p}}J",
             # we need to spell out SD_p somewhere
             "\\frac{1}{2}log \\left( \\frac{1+r}{1-r}\\right)")
             #"log\\left(\\frac{y_{i2}+0.5}{n_{i2}-y_{i2}+0.5}\\right) ",
             #"- log\\left(\\frac{y_{i1}+0.5}{n_{i1}-y_{i1}+0.5}\\right)")
             
# Just placeholders for now
eqns_vi <- c("\\frac{SD^2}{N}",
             "\\frac{1}{2(N - 1)}",
             "\\frac{SD_{1}^2} {M^2_{1}N_{1}} + \\frac{SD_{2}^2} {M^2_{2}N_{2}}",
             "\\frac{N_{1} + N_{2}}{N_{1}N_{2}} + \\frac{SMD^2}{2(N_{1} + N_{2})}",
             # see my commenst
             "\\frac{1}{N - 3}")
             #"\\frac{s_{i}^2}{n_{i}}",
             #"")

# Maybe add some examples for each effect size in a new column. SHould help
eg <- c("CTmin, CTmax, LC50, LT50, Metabolic Rate (MR)",
        "Variability in CTmin, CTmax, Metabolic Rate (MR)",
        "Ratio between pollutant exposed treatment (e.g., BPA-exposed group) and control (no pollutant)",
        "Difference in immune response between males and females, performance difference in the presence of stressor compared to absence of stressor",
        "Relationship between sex hormones and immune responses or metabolic rate and behaviour")

table_fin <- data.frame(names, eqns_yi, eqns_vi, eg)

footnotes <- c("Notes: J = 1 - \\frac{3}{4\\left(N_{1} + N_{2} - 2\\right) - 1}",
               "SD_{p} = \\sqrt \\frac{\\left( N_{1}-1 \\right)SD_{1}^2 + \\left( N_{2}-1 \\right)SD_{2}^2}{N_{1} + N_{2} - 2}")

# Note: It is esssential that you be clear that the function 'compose' is exported from the flextable package. That's because purr also has the same function, so it gets messy.
tablea <- flextable(table_fin)                                                                                %>% 
          flextable::compose(j = 1:4, part = "header", value = as_paragraph(as_b(c("Effect Measure", "Definition", "Sampling Variance", "Examples")))) %>% 
          flextable::compose(j = 2:3, value = as_paragraph(as_equation(c( eqns_yi, eqns_vi))))             %>%
          flextable::compose(j = 4, value = as_paragraph(eg))             %>%
          flextable::width(j = 1, 3)                                                                          %>% 
          flextable::width(j = 2, 2) %>% 
          flextable::width(j = 4, 5) %>% 
          flextable::footnote(i = 4, j = 1, value = as_paragraph(as_equation(footnotes)),
               ref_symbols = c("a"),
               part = "body", inline = TRUE) %>% 
          flextable::font(part = "all", fontname = "Times New Roman")
          
FitFlextableToPage(tablea, pgwidth = 8)
```

We'll use `metafor` to demonstrate how to calculate commonly used effect size estimates in the field of comparative physiology, including Z-transformed correlation coefficients (Zr), log response ratios and Hedge's g. 

## **Setting up**

To start, well need to load some packages and data. The data is all stored on GitHub. It's downloaded directly from the web. You won't need to have these on hand. We'll first load the packages that we need.

```{r}
#| label: load packages
#| eval: true
#| message: false
#| warning: false

# Load packages
#install.packages("pacman") # If you haven't already installed the pacman package do so with this code
pacman::p_load(metafor, flextable, tidyverse, orchaRd, pander, mathjaxr, equatags)
# To use mathjaxr you need to run equatags::mathjax_install()
```

## **Correlation between physiology and movement patterns: Z-transformed Correlation Coefficient (Zr)**

![](./figs/WuSeebacher.png)

**Background**: Physiology is expected to drive movement, including short-term activity, exploration of unfamiliar environments, and larger scale dispersal which can influence species distributions in an environmentally sensitive manner. To test whether there is a relationship between physiology and movement pattern (separated as activity, exploration, and dispersal), we extracted data from the literature on common physiological traits measured (e.g., metabolic rate, hormone levels, body condition) with either activity, exploration, and dispersal at the individual level. As we are interested in relationships at the individual level, we extracted correlation data in the form of [correlation coefficient](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient) and sample size.

For complete detail on the study aim and design, please see [@WuSeebacher2022](https://doi.org/10.1038/s42003-022-03055-y)

Download the data. The following code illustrates how to calculate the Fisher z-transformed effect size (Zr) and sampling variance (V_Zr) from correlation coefficients (ri) and sample size (ni) using 10 rows of data from the main dataset.

```{r}
#| label: Zr
#| message: false

# Download and clean data
zr_data <- read.csv("https://raw.githubusercontent.com/daniel1noble/meta-workshop/gh-pages/data/ind_disp_raw_data.csv") %>%
  dplyr::select(study_ID, taxa, species, trait, response, response_unit, disp_trait, disp_unit, corr_coeff, sample_size) %>% # remove irrelevant columns for this tutorial
  dplyr::top_n(10) # select first 10 effect sizes to illustrate escalc function

# Calculate Fisher's r-to-z transformed correlation coefficient (ZCOR) as yi = effect size and vi = sampling variances, where ri = raw correlation coefficients, and ni = sample size.
zr_data <- metafor::escalc(measure = "ZCOR", ri = corr_coeff, ni = sample_size, data = zr_data, var.names = c("Zr", "V_Zr"))

zr_data
```

We want to convert to Zr instead of using the correlation (r) itself because Zr will satisfy assumptions of normality (see below) that are inherent to meta-analytic models. We can, however, convert back to the original correlation coefficient quite easily as follows:

```{r}
#| label: Zrescalcbackcalc
#| echo: true

# We can easily convert back to r as follows
zr_data$r <- tanh(zr_data$Zr)

zr_data %>% select(corr_coeff, r)
```

The correlation coefficient itself is a unit less measure that captures both the strength and direction of relationship between two variables. The raw correlation coefficient ranges from -1 to 1. If a correlation is positive it indicates that, as the value of one variable increases so to does the value of the second variable. The magnitude of the value indicates how strong they co-vary. For example, a correlation of 1 indicates that variable 1 perfectly co-varies with variable 2.

## **Difference in physiology between populations from range core and edge: Log response ratio (lnRR)**

**Background**: If movement relied on physiological capacities, it may be expected that the variation in physiological characteristics introduces differences in the tendency to move among individuals of the same populations. Such individual differences can have consequences for gene flow and genetic diversification within a species. Individuals with a greater tendency to move may have particular physiological characteristics, leading to core-edge differences. Populations at the range edge is therefore hypothesised to have different physiological phenotype than populations from the core range. 

To compare differences in the physiological characteristics between individuals from an established range core to their expanding range edge we collected data from studies that reported physiological response from the range core (core) and range edge (front) in the form of mean, sample size, and variance (SD, SE, or 95% CI). This information is used to calculate the log response ratio (lnRR).

For complete detail on the study aim and design, please see [@WuSeebacher2022](https://doi.org/10.1038/s42003-022-03055-y)

Download the data. The following code illustrates how to calculate the log response ratio (`yi`) and sampling variance (vi) from the mean (`mi`), sample size (`ni`) and standard deviation (`sdi`) from the edge population (1), and core population (2) using 10 row of data from the main dataset.

```{r}
#| label: lnRR
#| message: false

# Download and clean data
ratio_data <- read.csv("https://raw.githubusercontent.com/daniel1noble/meta-workshop/gh-pages/data/pop_disp_raw_data.csv") %>%
  dplyr::select(study_ID, taxa, species, trait, response, response_unit, mean_core, sd_core, n_core, mean_front, sd_front, n_front) %>% # remove irrelevant columns for this tutorial
  dplyr::top_n(10) # select first 10 effect sizes to illustrate escalc function

# Calculate log response ratio (lnRR) as yi = effect size and vi = sampling variances, where m1i = mean of edge population, n1i = sample size of edge population, sd1i = standard deviation of edge population, m2i = mean of core population, n2i = sample size of core population, sd2i = standard deviation of core population.
ratio_data <- metafor::escalc(measure = "ROM", 
                              m1i = mean_front, n1i = n_front, sd1i = sd_front,
                              m2i = mean_core, n2i = n_core, sd2i = sd_core, 
                              data = ratio_data)
ratio_data
```

Now that we have calculated the log response ratio lets just make sure we understand what is being done. We'll take this subset of data and do some back calculations to make sure we understand how to inter convert. To demonstrate, we'll just convert (and back convert) the first row of data. 

#### **Testing your Understanding of $lnRR$**

::: {.panel-tabset}

## Task!

>**Using the first row of the data, calculate the log response ratio for this observation by hand. What does the direction of the effect size tell you about the mean at the marginal compared to core?**

## Answer!

>The effect size is positive indicating that the mean of the marginal (`mean_front`) is larger than the mean of the core poulaton (`mean_core`). We can see that by calculating the effect size in two ways, the second of which is easier to see why.
 
```{r}
#| label: ConvertLnRR
#| echo: true

### Let's back calculate effects to make sure we understand why they are interpreted as percentage differences

# First, lets make sure we understand how it's calculated
with(ratio_data, log(mean_front[1] / mean_core[1]))

# Alternatively....
with(ratio_data, log(mean_front[1]) - log(mean_core[1]))
```

## Task 2!

>**What does the value of the effect size tell you about how much the mean differences between marginal and core populations?**

## Answer 2!

>What this tells us is that the numerator is 1.15 times the denominator, or, that the marginal (front) mean is 15% larger compared to the core mean. We can see that this is true as follows:

```{r}
#| label: ConvertLnRR2
#| echo: true

# Now lets back-transform to odds. 
with(ratio_data, exp(log(mean_front[1]) - log(mean_core[1])))

# Interpretation: What this tells us is that the numerator is 1.15 times the denominator, or, that the marginal (front) mean is 15% larger compared to the core mean. We can see that this is true as follows:
with(ratio_data, exp(log(mean_front[1]) - log(mean_core[1])) * mean_core[1])

# This value should now match the marginal mean value, which it does
with(ratio_data, mean_front[1])
```

>A simple way to remember how to interpret this is to think that a converted *lnRR* = 1 means the two are the same. If you are 'subtracting` the marginal population from the core population and the effect size is positive, than 1.15 indicates that the marginal population is ~15% larger, whereas if the effect size is -1.15 than it indicates that the core population is 15% larger. 

<br>

Hopefully this is useful to interpret the meaning of the effect size data at the level of individual observations. The same principles will apply (for the most part) when converting back to the original scale from the mean estimates from meta-analytic models. 

:::

## **Cohen's d and Hedges g -- Visualising Their Meaning**

We could also use [@WuSeebacher2022](https://doi.org/10.1038/s42003-022-03055-y)'s data to calculate Hedge's g. Hedge's g is very similar to Cohen's d, but it has a small sample correction. This can also be done with `metafor::escalc()` by setting measure to "SMD".  We can use the same data here as all the summary statistics that are needed for calculating lnRR are the same as Hedge's g. We can calculate as follows. 

```{r hedges}
# Calculate hedges g as yi = effect size and vi = sampling variances, where m1i = mean of edge population, n1i = sample size of edge population, sd1i = standard deviation of edge population, m2i = mean of core population, n2i = sample size of core population, sd2i = standard deviation of core population.
hedge_data <- metafor::escalc(measure = "SMD", 
                              m1i = mean_front, n1i = n_front, sd1i = sd_front,
                              m2i = mean_core, n2i = n_core, sd2i = sd_core, 
                              data = ratio_data, var.names = c("d", "d_vi"))
hedge_data
```

Visualising and interpreting Hedge's g can be slightly less intuitive compared to log response ratios. There's a great [online resource](https://rpsychologist.com/cohend/) by Kristoffer Magnusson with some visuals that can help you decipher the meaning of Cohen's d effects, which is the same as Hedge's g. 

```{=html}
<iframe width="800px" height="500px" src="https://rpsychologist.com/cohend/?embed=1&CER=20&M0=100&M0lab=Control&M1=112&M1lab=Treatment&SD=15&c0=r-48_g-57_b-79_a-1&c1=r-0_g-0_b-0_a-1&c2=r-106_g-206_b-235_a-1&d=0.8&slider=1&xlab=Outcome" frameborder="0"></iframe>
```

## **References**

<div id="refs"></div>

<br>

## **Session Information**

```{r}
#| label: sessioninfo
#| echo: false

# Here, we'll create some R code
pander(sessionInfo(), locale = FALSE)
```