Nuisance heterogeneity
Daniel W.A. Noble, Essie Rodgers & Patrice Pottier
2022-07-04
Introduction
In the previous tutorial we over-viewed a few of the common effect size statistics likely to be used in comparative physiology. However, studies can vary in ways that can make the interpretation of these statistics challenging.
Nuisance heterogeneity – factors causing variation in effects that are blatantly obvious, like temperature and dosage differences applied across studies – are common in studies used in comparative physiology meta-analyses (Borenstein, 2019; Noble et al., 2022). Nuisance heterogeneity can be a big driver of effect size magnitude. For example, a contrast-based effect size, such as log response ratio, calculated from a study that has manipulated treatments to have a 10\(^{\circ}\)C temperature difference between treatments is likely going to give you a much larger effect size than a study that manipulated a treatment temperature difference of 2\(^{\circ}\)C. Amalgamating such effects within a meta-analysis will make the interpretation of effect magnitude complicated because it’s difficult to know what magnitude of, say temperature, our effect size should be interpreted over.
There are a few different ways to deal with nuisance heterogeneity. The most common approach is to simply include a moderator in a meta-regression model that controls for the temperature difference between treatments (Noble et al., 2022). You can re-scale this continuous moderator (center the variable on the mean so ‘0’ equates to the average difference) so that intercepts can be interpreted more easily (Noble et al., 2022). Alternatively, we can ‘correct’ our effect sizes to account for the nuisance variable. In this tutorial, we’ll overview the latter approach as it can simplify modelling.
Effect Sizes Corrected for Nuisance Heterogeneity
Temperature Coefficient, \(Q_{10}\)
One common and well known effect size that is used in comparative physiology is the temperature coefficient, \(Q_{10}\). As discussed, effect sizes need to be: 1) comparable across studies and 2) we need to know, or be able to approximate, the sampling variance in order to make use of meta-analytic models. While \(Q_{10}\) satisfies 1) we do now know what the sampling variance should be for this effect size. Noble et al. (2022) fixes that issue by deriving it’s sampling variance.
To begin, recall that \(Q_{10}\), is defined as:
\[ \begin{equation} Q_{10} = \left( \frac{R_{2}}{R_{1}} \right)^{ \left(\ \frac{10^{\circ}C}{T_{2}-T_{1}} \right) } \tag{1} \end{equation} \]
Here, \(R_{1}\) and \(R_{2}\) are mean physiological rates and \(T_{1}\) and \(T_{2}\) are the temperatures that these rates are measured. Natural log transformation of equation (1) leads to the following log transformed \(Q_{10}\):
\[ \begin{equation} lnRR_{Q_{10}} = ln\left( \frac{R_{2}}{R_{1}} \right){ \left(\ \frac{10^{\circ}C}{T_{2}-T_{1}} \right) } \tag{2} \end{equation} \]
Noble et al. (2022) call this transformed effect size \(LnRR_{Q_{10}}\) because of its similarities with the log response ratio. Equation (2) is essentially a temperature-corrected equivalent of lnRR when the numerator and denominator are measured at different temperatures. \(lnRR_{Q_{10}}\) is more likely to satisfy the assumption of residual normality than \(Q_{10}\) and given it’s similarities to lnRR we can now approximate it’s sampling variance (see Hedges et al., 1999). The recognition of this equivalence means that we can calculate the sampling variance for equation (2) as follows:
\[ \begin{equation} s^2_{lnRR_{Q_{10}}} = \left( \frac{SD_{2}^2}{R^2_{2}N_{2}} + \frac{SD_{1}^2}{R^2_{1}N_{1}} \right){ \left(\ \frac{10^{\circ}C}{T_{2}-T_{1}} \right) }^2 \tag{3} \end{equation} \]
\(lnRR_{Q_{10}}\): Diving vertebrates and climate change
We will demonstrate how to apply \(lnRR_{Q_{10}}\) to a data set on diving vertebrates by Rodgers et al. (2021). Rodgers et al. (2021) already analysed \(Q_{10}\), but at the time they did not have \(lnRR_{Q_{10}}\), meaning it was not possible to analyse these data using meta-analytic models – they resorted to meta-regression approaches instead with standard effect sizes, which work equally as well. Here, we will show you how to calculate \(lnRR_{Q_{10}}\) and interpret it.
Download Data and Functions for \(lnRR_{Q_{10}}\)
We’ll download the data set used to analyse \(Q_{10}\) and load a bunch of packages that are necessary.
q10_dat <- read.csv("https://osf.io/download/fb3ht/")
# install.packages('pacman')
pacman::p_load(devtools, pander, tidyverse, metafor)We’ll now need a function for calculating \(lnRR_{Q_{10}}\) because this effect size is not part of the metafor package.
#' @title lnRR_Q10
#' @description Calculates the log Q10 response ratio. Note that temperature 2 is placed in the numerator and temperature 1 is in the denominator
#' @param t1 Lowest of the two treatment temperatures
#' @param t2 Highest of the two treatment temperatures
#' @param r1 Mean physiological rate for temperature 1
#' @param r2 Mean physiological rate for temperature 2
#' @param sd1 Standard deviation for physiological rates at temperature 1
#' @param sd2 Standard deviation for physiological rates at temperature 2
#' @param n1 Sample size at temperature 1
#' @param n2 Sample size at temperature 2
#' @param name Character string for column name
#' @example
#' lnRR_Q10(20, 30, 10, 5, 1, 1, 30, 30)
#' lnRR_Q10(20, 30, 10, 5, 1, 1, 30, 30, name = 'acclim')
#' @export
lnRR_Q10 <- function(t1, t2, r1, r2, sd1, sd2, n1, n2, name = "acute") {
lnRR_Q10 <- (10/(t2 - t1)) * log(r2/r1)
V_lnRR_Q10 <- (10/(t2 - t1))^2 * ((sd1^2/(n1 * r1^2)) + (sd2^2/(n2 * r2^2)))
dat <- data.frame(lnRR_Q10, V_lnRR_Q10)
colnames(dat) <- c(paste0("lnRR_Q10", name), paste0("V_lnRR_Q10", name))
return(dat)
}Calculating and Interpreting \(lnRR_{Q_{10}}\)
We can now use the function above to calculate \(lnRR_{Q_{10}}\).
# Calculate lnRRQ10 and it's associated sampling variance
q10_dat <- cbind(q10_dat, with(q10_dat, lnRR_Q10(t1 = t1, t2 = t2, r1 = mean_t1,
r2 = mean_t2, sd1 = sd_t1, sd2 = sd_t2, n1 = n_t1, n2 = n_t2, name = "")))
head(q10_dat)## study_ID study_name year species order life_stage
## 1 119 119_Udyawer_etal_2016 2016 Hydrophis curtus serpentes unspecified
## 2 119 119_Udyawer_etal_2016 2016 Hydrophis curtus serpentes unspecified
## 3 119 119_Udyawer_etal_2016 2016 Hydrophis curtus serpentes unspecified
## 4 119 119_Udyawer_etal_2016 2016 Hydrophis elegans serpentes unspecified
## 5 119 119_Udyawer_etal_2016 2016 Hydrophis elegans serpentes unspecified
## 6 119 119_Udyawer_etal_2016 2016 Hydrophis elegans serpentes unspecified
## animal_source body_mass respiration_mode replication_level dive.type
## 1 wild 160 bimodal individual voluntary
## 2 wild 160 bimodal individual voluntary
## 3 wild 160 bimodal individual voluntary
## 4 wild 1035 bimodal individual voluntary
## 5 wild 1035 bimodal individual voluntary
## 6 wild 1035 bimodal individual voluntary
## summary_stat acclimation_temp t1 t2 mean_t delta_t t_magnitude mean_t1 sd_t1
## 1 mean 27 21 24 22 3 plus3 40 16.1
## 2 mean 27 24 27 26 3 plus3 28 15.4
## 3 mean 27 27 30 28 3 plus3 16 4.9
## 4 mean 27 21 24 22 3 plus3 27 13.3
## 5 mean 27 24 27 26 3 plus3 22 10.0
## 6 mean 27 27 30 28 3 plus3 14 4.7
## n_t1 se_t1 mean_t2 sd_t2 n_t2 se_t2 shared_control Q10 source lnRR_Q10
## 1 11 4.8 28 15.4 11 4.63 1 0.33 Fig1B -1.12
## 2 11 4.6 16 4.9 11 1.47 2 0.14 Fig1B -1.96
## 3 11 1.5 11 3.5 11 1.05 3 0.31 Fig1B -1.16
## 4 10 4.2 22 10.0 10 3.16 1 0.46 Fig1B -0.78
## 5 10 3.2 14 4.7 10 1.47 2 0.24 Fig1B -1.43
## 6 10 1.5 11 2.0 10 0.63 3 0.43 Fig1B -0.84
## V_lnRR_Q10
## 1 0.46
## 2 0.39
## 3 0.20
## 4 0.50
## 5 0.36
## 6 0.16# Now we can back-calculate to put on the original Q10 scale. We can than check
# that this matches the Q10 already in the data
head(q10_dat %>%
mutate(exp_lnRR_Q10 = exp(lnRR_Q10)) %>%
select(exp_lnRR_Q10, Q10))## exp_lnRR_Q10 Q10
## 1 0.33 0.33
## 2 0.14 0.14
## 3 0.31 0.31
## 4 0.46 0.46
## 5 0.24 0.24
## 6 0.43 0.43Now that we have our \(lnRR_{Q_{10}}\) effect size estimates which we back-transformed how do we interpret them?
Task
How do we interpret \(lnRR_{Q_{10}}\)?
Answer
Given these are simply \(Q_{10}\)’s we can interpret these normally. For the first row of our data we can see that dive duration when temperatures increase by 10\(^{\circ}\)C is expected to decrease by ~67% (1-0.33). In other words, the dive duration 10\(^{\circ}\)C above the first temperature is 33% of the rate 10\(^{\circ}\)C lower.
Acclimation Response Ratio, ARR
Pottier et al. (2021) used the acclimation response ratio, ARR, to investigate variation in the relationship between acclimation temperature and thermal tolerance (i.e. variation in slope). Because heat tolerance is consistently measured in degrees Celsius (CTmax) in the studies compiled, the slope between acclimation temperature and heat tolerance can be calculated, which is termed as the acclimation response ratio (ARR) in the literature. It can be defined as follows:
\[ \text{ARR} = \frac{M_{1} - M_{2}}{T_2 - T_1}, \] where, \(M_{i}\) is the mean for treatment 1 and 2 while \(T_{i}\) is the temperature of treatment 1 and 2.
When the two groups are independent than the sampling variance for \(ARR\) is:
\[ \sigma^2_{ARR} = \left( \frac{1} {T_2 - T_1} \right)^2 \left( \frac{SD^2_{1}} {N_{1}} + \frac{SD^2_{2}} {N_{2}} \right) \]
In this example, we will use a simplified data set from Pottier et al. (2021), and demonstrate how to investigate the mean effect of acclimation temperature on heat tolerance.This effect size is calculated as the difference between the heat tolerance at the higher (mean_high) and the lower (mean_low) acclimation temperatures, divided by the difference between acclimation temperatures (acc_temp_high - acc_temp_low).
Download Data and Functions for \(ARR\)
asr_dat <- read.csv("https://osf.io/qn2af/download")We’ll also need a function for calculating ARR and its sampling varinace because these don’t exist in any current packages.
#' @title arr
#' @description Calculates the acclimation response ratio (ARR).
#' @param t2_l Lowest of the two treatment temperatures
#' @param t1_h Highest of the two treatment temperatures
#' @param x1_h Mean trait value at high temperature
#' @param x2_l Mean trait value at low temperature
#' @param sd1_h Standard deviation of mean trait value at high temperature
#' @param sd2_l Standard deviation of mean trait value at low temperature
#' @param n1_h Sample size at high temperature
#' @param n2_l Sample size at low temperature
arr <- function(x1_h, x2_l, sd1_h, sd2_l, n1_h, n2_l, t1_h, t2_l) {
ARR <- (x1_h - x2_l)/(t1_h - t2_l)
V_ARR <- ((1/(t1_h - t2_l))^2 * (sd2_l^2/n2_l + sd1_h^2/n1_h))
return(data.frame(ARR, V_ARR))
}Calculating and Interpreting \(ARR\)
Using the arr function above we can now calculate the acclimation response ratio (ARR)
# Calculate the effect sizes
asr_dat <- asr_dat %>%
mutate(ARR = arr(x1_h = mean_high, x2_l = mean_low, t1_h = acc_temp_high, t2_l = acc_temp_low,
sd1_h = sd_high, sd2_l = sd_low, n1_h = n_high_adj, n2_l = n_low_adj)[, 1],
V_ARR = arr(x1_h = mean_high, x2_l = mean_low, t1_h = acc_temp_high, t2_l = acc_temp_low,
sd1_h = sd_high, sd2_l = sd_low, n1_h = n_high_adj, n2_l = n_low_adj)[,
2])
head(asr_dat %>%
select(ARR, V_ARR))## ARR V_ARR
## 1 0.060 1.6e-04
## 2 0.089 1.1e-04
## 3 0.086 1.8e-04
## 4 0.090 9.3e-05
## 5 0.082 1.0e-04
## 6 0.062 7.8e-05Now try to interpret \(ARR\) effect size estimates
Task
How do we interpret \(ARR\)? for the first row of the data
Answer
If we look at the first row of data, we can see that for each degree increase in acclimation temperature, the heat tolerance of the animals included in this data set increases by 0.06\(^{\circ}\)C.
References
Session Information
R version 4.0.5 (2021-03-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: metafor(v.3.5-4), metadat(v.1.2-0), Matrix(v.1.4-0), forcats(v.0.5.1), stringr(v.1.4.0), dplyr(v.1.0.9), purrr(v.0.3.4), readr(v.2.1.2), tidyr(v.1.2.0), tibble(v.3.1.7), ggplot2(v.3.3.6), tidyverse(v.1.3.1), pander(v.0.6.4), devtools(v.2.4.3) and usethis(v.2.1.5)
loaded via a namespace (and not attached): nlme(v.3.1-152), fs(v.1.5.2), lubridate(v.1.8.0), httr(v.1.4.3), rprojroot(v.2.0.2), tools(v.4.0.5), backports(v.1.4.1), bslib(v.0.3.1), utf8(v.1.2.2), R6(v.2.5.1), DBI(v.1.1.3), colorspace(v.2.0-3), withr(v.2.5.0), tidyselect(v.1.1.2), prettyunits(v.1.1.1), processx(v.3.5.3), curl(v.4.3.2), compiler(v.4.0.5), cli(v.3.3.0), rvest(v.1.0.2), formatR(v.1.11), pacman(v.0.5.1), xml2(v.1.3.3), desc(v.1.4.0), bookdown(v.0.24), sass(v.0.4.1), scales(v.1.2.0), callr(v.3.7.0), digest(v.0.6.29), rmarkdown(v.2.14), pkgconfig(v.2.0.3), htmltools(v.0.5.2), sessioninfo(v.1.2.2), dbplyr(v.2.2.0), fastmap(v.1.1.0), rlang(v.1.0.2), readxl(v.1.4.0), rstudioapi(v.0.13), jquerylib(v.0.1.4), generics(v.0.1.2), jsonlite(v.1.8.0), magrittr(v.2.0.3), Rcpp(v.1.0.8.3), munsell(v.0.5.0), fansi(v.1.0.3), lifecycle(v.1.0.1), stringi(v.1.7.6), yaml(v.2.3.5), mathjaxr(v.1.6-0), pkgbuild(v.1.3.1), grid(v.4.0.5), crayon(v.1.5.1), lattice(v.0.20-45), haven(v.2.5.0), hms(v.1.1.1), knitr(v.1.39), ps(v.1.6.0), klippy(v.0.0.0.9500), pillar(v.1.7.0), pkgload(v.1.2.4), reprex(v.2.0.1), glue(v.1.6.2), evaluate(v.0.15), remotes(v.2.4.2), modelr(v.0.1.8), vctrs(v.0.4.1), tzdb(v.0.3.0), testthat(v.3.1.1), cellranger(v.1.1.0), gtable(v.0.3.0), assertthat(v.0.2.1), cachem(v.1.0.6), xfun(v.0.31), broom(v.0.8.0), memoise(v.2.0.1) and ellipsis(v.0.3.2)