For a chosen effect size, it’s point estimate ($z_j$) and sampling variance ($\nu_j$) can be calculated using function escalc() in metafor package.

Let’s load the dataset of our first worked example from [1], who meta-analysed the intraspecific change in seven morpho‐ecophysiological leaf traits along global elevational gradients.

dat_Midolo_2019 <- read.csv(here("data","Midolo_2019_Global Change Biology.csv"))

For simplicity, we will only keep the columns that are necessary to demonstrate the meta-analytic techniques of interest in this tutorial.

Table S1
The variables in the first worked example ([1]).

Table S1 includes the coded variables that were extracted from published papers.

You can find more details about the variables in the supplementary information of the paper by [1], Tables S1, if interested. Here we classify the variables into three categories types:

1. Bibliographic variables

Variables that distinguish between articles, such as the title (study_name), or publication year.

2. Descriptive statistics

Variables that are used to compute effect sizes, such as mean (treatment and control), standard deviation (sd_treatment and sd_control) and sample size (n_treatment and n_control).

3. Study-related characteristics

Study-level variables that can contribute to random variability among the effect sizes, such as study identity (Study_ID), or systematic variability between the effect sizes, such as treatment intensity, or study location characteristics such as (elevation); which may be of interest as a moderator variable.

After having an (primary) understanding of the extracted/coded variables, we can use escalc() to compute effect size (using lnRR as an example) with.

Note that the argument measure is used to specify the effect size metric that we will calculate. For example, by specifying measure = "SMD" , the function will compute the point estimate of SMD (specifically, Hedges’ g) and its sampling variance.

Here are the most commonly used effect sizes and underappreciated effect sizes. For other types, you can look at the help information of escalc by typing ?escalc in the console of R.

• measure = "ROM": lnRR - quantifies the difference between the means of two groups

• measure = "SMD": SMD - quantifies the differences between the means of two groups

• measure = "ZCOR": Zr - quantifies the strength of association between two variables

• measure = "VR": lnVR - quantifies the differences in variances between two groups

• measure = "CVR": lnCVR - quantifies the differences in variances between two groups, accounting for the mean-variance relationship

Note that lnVR and lnCVR are variation/dispersion-based effect sizes, which have practical implications for environmental studies. For example, environmental stressors such as pesticides and eutrophication are likely to increase variability in biological systems because stress accentuates individual differences in environmental responses.

lnRR <- escalc(measure = "ROM",  # "ROM" means ratio of means; lnRR is specified to be calculated (alternative effect sizes: "SMD" – SMD, "CVR" – lnCVR, "VR" – lnVR; see below);
               m1i = treatment, # mean value of of group 1 (e.g., environmental stressor); in our worked example, m1i denotes the mean value of a trait measured at the higher elevation level;
               m2i = control, # mean value of group 2 (e.g., control); in our worked example, m2i denotes the mean of the same trait measured at the lower elevation level;
               sd1i = sd_treatment, # standard deviation of mean of group 1 (e.g., environmental stressor)
               sd2i = sd_control, # standard deviation of group 2 (e.g., control) 
               n1i = n_treatment, # sample size of group 1 (e.g., environmental stressor) 
               n2i = n_control, # sample size of group 2 (e.g., control) 
               data = dat_Midolo_2019, # dataset containing the above information (here is the dataset of our working example)
               )
               

Let’s explore the effect size estimates and their sampling variances.

Table S2

The point estimate of effect size (lnRR) and their sampling variance (lnRRV) for each included study and comparison (higher elevation vs. lower elevation).

i2_ml(mod_ML_lnRR_VCV)

   I2_Total I2_Study_ID    I2_ES_ID 
   99.63261    17.97605    81.65656 

mod_MLMR_lnRR_elevation <- rma.mv(yi = lnRR, 
                                  V = VCV, 
                                  mods = ~ elevation_log, # adding elevation as a moderator variable (for continuous variable, it is a good practice to log-transform them to avoid data skewness);
                                  random = list(~1 | Study_ID, 
                                                ~1 | ES_ID), 
                                  method = "REML", 
                                  test = "t", 
                                  data = dat2_Midolo_2019,
                                  sparse = TRUE
                                  )

Multivariate Meta-Analysis Model (k = 1294; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-103.6834   207.3669   215.3669   236.0227   215.3979   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0098  0.0991     71     no  Study_ID 
sigma^2.2  0.0515  0.2269   1294     no     ES_ID 

Test for Residual Heterogeneity:
QE(df = 1292) = 65188.5301, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 1292) = 8.0348, p-val = 0.0047

Model Results:

               estimate      se     tval    df    pval    ci.lb    ci.ub     
intrcpt         -0.1507  0.0647  -2.3290  1292  0.0200  -0.2776  -0.0238   * 
elevation_log    0.0290  0.0102   2.8346  1292  0.0047   0.0089   0.0491  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod_MLMR_lnRR_trait <- rma.mv(yi = lnRR, 
                              V = VCV, 
                              mods = ~ trait -1, # setting '-1' is a useful strategy to remove intercept and then directly obtain the estimated effect (slope) for each subgroup or for a particular level within the categorical moderator.
                              random = list(~1 | Study_ID, 
                                                ~1 | ES_ID), 
                              method = "REML", 
                              test = "t", 
                              data = dat2_Midolo_2019
                              )

Multivariate Meta-Analysis Model (k = 1294; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-71.6727  143.3454  161.3454  207.7860  161.4863   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0112  0.1058     71     no  Study_ID 
sigma^2.2  0.0486  0.2205   1294     no     ES_ID 

Test for Residual Heterogeneity:
QE(df = 1287) = 63967.9953, p-val < .0001

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 1287) = 11.4776, p-val < .0001

Model Results:

               estimate      se     tval    df    pval    ci.lb    ci.ub      
I(trait)dC13    -0.0136  0.0243  -0.5596  1287  0.5758  -0.0613   0.0341      
I(trait)LA      -0.0672  0.0275  -2.4493  1287  0.0144  -0.1211  -0.0134    * 
I(trait)LMA      0.0617  0.0251   2.4566  1287  0.0142   0.0124   0.1109    * 
I(trait)Narea    0.0820  0.0243   3.3816  1287  0.0007   0.0344   0.1296  *** 
I(trait)Nmass    0.0842  0.0214   3.9409  1287  <.0001   0.0423   0.1261  *** 
I(trait)Pmass    0.1633  0.0332   4.9155  1287  <.0001   0.0981   0.2285  *** 
I(trait)SLA     -0.0258  0.0241  -1.0702  1287  0.2848  -0.0732   0.0215      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod_MLMR_lnRR_trait2 <- rma.mv(yi = lnRR, 
                              V = VCV, 
                              mods = ~ I(trait),
                              random = list(~1 | Study_ID, 
                                            ~1 | ES_ID), 
                              method = "REML", 
                              test = "t", 
                              data = dat2_Midolo_2019
                              )

Multivariate Meta-Analysis Model (k = 1294; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-71.6727  143.3454  161.3454  207.7860  161.4863   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0112  0.1058     71     no  Study_ID 
sigma^2.2  0.0486  0.2205   1294     no     ES_ID 

Test for Residual Heterogeneity:
QE(df = 1287) = 63967.9953, p-val < .0001

Test of Moderators (coefficients 2:7):
F(df1 = 6, df2 = 1287) = 12.9754, p-val < .0001

Model Results:

               estimate      se     tval    df    pval    ci.lb   ci.ub      
intrcpt         -0.0136  0.0243  -0.5596  1287  0.5758  -0.0613  0.0341      
I(trait)LA      -0.0536  0.0304  -1.7666  1287  0.0775  -0.1132  0.0059    . 
I(trait)LMA      0.0753  0.0283   2.6639  1287  0.0078   0.0198  0.1307   ** 
I(trait)Narea    0.0956  0.0265   3.6063  1287  0.0003   0.0436  0.1477  *** 
I(trait)Nmass    0.0978  0.0242   4.0417  1287  <.0001   0.0503  0.1452  *** 
I(trait)Pmass    0.1769  0.0354   5.0013  1287  <.0001   0.1075  0.2463  *** 
I(trait)SLA     -0.0122  0.0285  -0.4285  1287  0.6684  -0.0681  0.0437      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The goodness-of-fit index $R^2$ can be used to quantify the percentage of variance explained by a moderator. This index has implications of the importance of the examined moderator. It also can be used for model section to select a model with the ‘best’ set of moderators. A general and readily implementable form of $R^2$ is the marginal version:

[2] propose to use a general form of $R^2$ - marginal $R^2$, which can be calculated as:

$$R^2_{marginal}=\frac{f^2} {f^2+\tau^2+\sigma^2}, (12) \\ f^2 = Var(\sum \beta_{h}x_{h[i]}), (13)$$ We wrote a function r2_ml() (in orchaRd package), which does the calculation of $R^2_{marginal}$ :

r2_ml(mod_MLMR_lnRR_elevation)

The result given under R2_marginal shows that elevation can explain 0.98% of the variation between effect sizes.

   R2_marginal R2_conditional 
   0.009803104    0.168309094 

forest(mod_ML_lnRR_VCV, # rma.mv object
       xlab = "Effect size lnRR")

Function orchard_plot() in orchaRd package can make a forest-like plot (termed as ‘orchard’ plot) that can accommodate large number of studies.

A basic orchard plot can be made with:

orchard_plot(mod_multilevel_SMD, 
             mod = "1", 
             xlab = "Standardised mean difference (SMD)", 
             group = "Study_ID",  k = TRUE, g = TRUE,
             data = dat2_Midolo_2019) + 
             scale_x_discrete(labels = c("Overall effect (meta-analytic lnRR)"))

Figure S2
Orchard plot (forest-like plot) showing the effect sizes from each study (and their 95% CIs), the overall effect and its 95% CIs and 95% prediction intervals.

We see that an orchard plot is more informative than a classic forest plot. For example, it shows the distribution of the effect sizes, the number of effect sizes (k) and the number of studies (the number in the bracket). An orchard plot not only can visualize the meta-analytic results, but also lets us realize something we are not bale to see from statistical results, such as influential data points and outliers that could threaten the robustness of our results. Further, an orchard plot not only show 95% CIs (thick whiskers in Figure S2) but also 95% prediction intervals (PIs; thin whiskers), through which we can visually check the heterogeneity among effect sizes. Because CIs only include the standard error of the overall effect $\beta_0$ ($SE[\beta_0]$), which PIs include variance of random effects:

$$\text{95%CI} = \beta_{0} \pm t_{df[\alpha=0.05]} SE[\beta_{0}^2], (14)$$

$$\text{95%PI} = \beta_{0} \pm t_{df[\alpha=0.05]} \sqrt{\tau^2+\sigma^2+ SE[\beta_{0}^2]}, (15)$$ If you want to know the detailed value of PIs, you can use predict:

   pred     se   ci.lb  ci.ub   pi.lb  pi.ub 
 0.0268 0.0170 -0.0066 0.0602 -0.4667 0.5203 

• Visualize meta-regression with a categorical moderator

A fabulous function of orchard_plot() is that it also can nicely show results based on meta-regression model (Figure S3). The results based on meta-regression with trait as a moderator can be visualized with:

orchard_plot(mod_MLMR_lnRR_trait, 
             mod = "trait", 
             xlab = "Effect size lnRR", 
             group = "Study_ID",  k = TRUE, g = TRUE, trunk.size = 1.5, 
             data = dat2_Midolo_2019) + 
             scale_x_discrete(labels = c("SLA","Pmass","Nmass","Narea","LMA","LA","dC13"))

Figure S3
Orchard plot (forest-like plot) showing the effect sizes from each study and the overall effect for each subgroup/levels of a given categorical moderator.

• Visualize meta-regression with a continuous moderator

For a meta-regression with a continuous moderator, we can use bubble_plot() to visualize the results. Let’s visualize the meta-regression with elevation as the continuous moderator:

bubble_plot(mod_MLMR_lnRR_elevation, 
            mod = "elevation_log", 
            xlab = "Elevation (log-transformed)", ylab = "Effect size lnRR",
            group = "Study_ID",k = TRUE, g = TRUE,
            data = dat2_Midolo_2019, legend.pos = "top.left") +
            theme(axis.text.x = element_text(size = 12, colour = "black"),
                  axis.text.y = element_text(size = 12, colour = "black"),
                  axis.title.x = element_text(size = 12, colour = "black"),
                  plot.title = element_text(size = 12, colour = "black"))

Figure S4
Bubble plot showing the results of a meta-regression model: relationship between a continuous moderator (in our case, elevation) and effect size magnitude (lnRR).

ess.var_cal <- function(dat){1/dat$n_control + 1/dat$n_treatment} # write a help function to calculate adapted sampling variance based on effective size based  - tilde n
dat2_Midolo_2019$ess.var <- ess.var_cal(dat2_Midolo_2019) # calculate tilde N
dat2_Midolo_2019$ess.se <- sqrt(dat2_Midolo_2019$ess.var) # calculate adapted sampling error based on effective size - tilde square root n

mod_MLMR_lnRR_ess.se <- rma.mv(yi = lnRR, 
                               V = VCV, 
                               mods = ~ ess.se, # add adjusted based sampling error - tilde square root n as a moderator to test small study effect. 
                               random = list(~1 | Study_ID, 
                                             ~1 | ES_ID), 
                               method = "REML", 
                               test = "t", 
                               data = dat2_Midolo_2019
                               )

Multivariate Meta-Analysis Model (k = 1294; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-106.2100   212.4201   220.4201   241.0758   220.4511   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0112  0.1060     71     no  Study_ID 
sigma^2.2  0.0516  0.2272   1294     no     ES_ID 

Test for Residual Heterogeneity:
QE(df = 1292) = 65891.3832, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 1292) = 1.6336, p-val = 0.2014

Model Results:

         estimate      se     tval    df    pval    ci.lb   ci.ub    
intrcpt   -0.0353  0.0515  -0.6865  1292  0.4925  -0.1364  0.0657    
ess.se     0.1231  0.0963   1.2781  1292  0.2014  -0.0658  0.3120    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod_MLMR_lnRR_Year.c <- rma.mv(yi = lnRR, 
                               V = VCV, 
                               mods = ~ Year.c, # add centered year as a moderator to detect decline effect.
                               random = list(~1 | Study_ID, 
                                             ~1 | ES_ID), 
                               method = "REML", 
                               test = "t", 
                               data = dat2_Midolo_2019,
                               sparse = TRUE
                               )

Multivariate Meta-Analysis Model (k = 1294; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-106.5039   213.0079   221.0079   241.6637   221.0390   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0114  0.1067     71     no  Study_ID 
sigma^2.2  0.0517  0.2273   1294     no     ES_ID 

Test for Residual Heterogeneity:
QE(df = 1292) = 65896.5620, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 1292) = 0.8345, p-val = 0.3611

Model Results:

         estimate      se    tval    df    pval    ci.lb   ci.ub    
intrcpt    0.0267  0.0170  1.5668  1292  0.1174  -0.0067  0.0601    
Year.c     0.0020  0.0021  0.9135  1292  0.3611  -0.0022  0.0062    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod_MLMR_lnRR_ess.var <- rma.mv(yi = lnRR, 
                                V = VCV, 
                                mods = ~ ess.var, # adding adjusted sampling error to obtain publication-bias-corrected overall effect - which is potentially regarded as 'true' effect.
                                random = list(~1 | Study_ID, 
                                              ~1 | ES_ID), 
                                method = "REML", 
                                test = "t", 
                                data = dat2_Midolo_2019,
                                sparse = TRUE
                                )

Multivariate Meta-Analysis Model (k = 1294; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-106.4500   212.8999   220.8999   241.5557   220.9310   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0113  0.1065     71     no  Study_ID 
sigma^2.2  0.0516  0.2272   1294     no     ES_ID 

Test for Residual Heterogeneity:
QE(df = 1292) = 65878.9167, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 1292) = 1.1843, p-val = 0.2767

Model Results:

         estimate      se     tval    df    pval    ci.lb   ci.ub    
intrcpt   -0.0029  0.0321  -0.0897  1292  0.9285  -0.0659  0.0602    
ess.var    0.1048  0.0963   1.0883  1292  0.2767  -0.0841  0.2937    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multivariate Meta-Analysis Model (k = 1294; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-67.0224  134.0449  158.0449  219.9377  158.2904   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0099  0.0996     71     no  Study_ID 
sigma^2.2  0.0486  0.2204   1294     no     ES_ID 

Test for Residual Heterogeneity:
QE(df = 1284) = 63164.1842, p-val < .0001

Test of Moderators (coefficients 1:10):
F(df1 = 10, df2 = 1284) = 8.9597, p-val < .0001

Model Results:

               estimate      se     tval    df    pval    ci.lb    ci.ub      
ess.var          0.0119  0.0936   0.1275  1284  0.8985  -0.1717   0.1956      
Year.c           0.0021  0.0021   1.0005  1284  0.3172  -0.0020   0.0061      
elevation_log    0.0284  0.0100   2.8385  1284  0.0046   0.0088   0.0481   ** 
traitdC13       -0.1900  0.0707  -2.6891  1284  0.0073  -0.3286  -0.0514   ** 
traitLA         -0.2446  0.0713  -3.4329  1284  0.0006  -0.3844  -0.1048  *** 
traitLMA        -0.1165  0.0710  -1.6399  1284  0.1013  -0.2559   0.0229      
traitNarea      -0.0962  0.0708  -1.3597  1284  0.1742  -0.2351   0.0426      
traitNmass      -0.0939  0.0699  -1.3432  1284  0.1795  -0.2311   0.0433      
traitPmass      -0.0150  0.0746  -0.2016  1284  0.8403  -0.1614   0.1313      
traitSLA        -0.2035  0.0698  -2.9160  1284  0.0036  -0.3404  -0.0666   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

profile(mod_MLMR_lnRR_whole, sigm2=1, progbar=FALSE) # check the estimates of between-study ($\sigma^2_1$ on the following figure) and within-study variance ($\sigma^2_2$) 

• Sensitivity of assumption of within-study (sampling) correlation

Here, we show how to perform a sensitivity analysis to examine the extent to which the overall effect (e.g., $\beta_{0}$) is sensitive to the assumption of within-study (sampling) correlation $\rho$ values used for constructing variance-covariance matrix (VCV).

First, set a series of $\rho$ (i.e., 0.3, 0.5, 0.7, 0.9) (we assume these values arbitrarily; you can set them based on your expertise in your field):

rho_range <- c(0.3, 0.5, 0.7, 0.9)

Then, we write a function to help repeatedly run the specified model, changing $\rho$ at a time:

ML_VCV_range <- list() # repeatedly run the specified model with varying rho
for (i in 1:length(rho_range)) {

VCV_range <- vcalc(vi = lnRRV, cluster = Study_ID, obs = ES_ID, 
                  rho = rho_range[i],
                  data = dat2_Midolo_2019
                  ) # impute VCV matrix with variying rho
              
ML_VCV_range[[i]] <- rma.mv(yi = lnRR, 
                            V = VCV_range, # VCV matrix with varying values of rho. 
                            random = list(~1 | Study_ID, 
                                          ~1 | ES_ID), 
                            method = "REML", 
                            test = "t", 
                            data = dat2_Midolo_2019
                           )} # run model with different rho values.

From Table S5, we see that the overall effect $\beta_{0}$ does not change with the changing of $\rho$ values, indicating that the model estimates are robust to different assumption of $\rho$.

Table S5 Sensitivity analysis examining the robustness of the overall effect (i.e., $\beta_0$) to the assumption of $\rho$ values.

• Leave-one-out analysis

The sensitivity (robustness of results) can be evaluated using leave-one-out analysis:repeatedly fitting the model to the dataset with one unit/entry (e.g., study) removed at a time. ‘One unit/entry’ denotes a row from the dataset fed to the model.

Leave-one-out analysis is a useful method to diagnose outlier and influential studies (which exert disproportionate influence on model estimates like the overall effect $\beta_0$ and its 95% CIs). As you might imagine, leave-one-out analysis re-fit the model N (the number of study) times. So it may take quite a long time to finish such a analysis if you have a large N.

For fixed and random effects models, the implementation of leave-one-out analysis can be easily performed using the leave1out() function from package metafor. The syntax is just one line:

leave1out(mod_FE_lnRR)

Let’s have a look at the results of the first 10 rows:

   estimate     se     zval   pval   ci.lb   ci.ub          Q     Qp      I2 
1   -0.0144 0.0004 -34.0889 0.0000 -0.0152 -0.0135 41265.3873 0.0000 96.8690 
2   -0.0144 0.0004 -34.0869 0.0000 -0.0152 -0.0135 41264.7609 0.0000 96.8690 
3   -0.0144 0.0004 -34.0954 0.0000 -0.0152 -0.0135 41266.2873 0.0000 96.8691 
4   -0.0144 0.0004 -34.0948 0.0000 -0.0152 -0.0135 41266.5291 0.0000 96.8691 
5   -0.0143 0.0004 -34.0330 0.0000 -0.0152 -0.0135 41257.2615 0.0000 96.8684 
6   -0.0143 0.0004 -33.9758 0.0000 -0.0151 -0.0135 41220.5405 0.0000 96.8656 
7   -0.0144 0.0004 -34.0714 0.0000 -0.0152 -0.0135 41264.3804 0.0000 96.8690 
8   -0.0143 0.0004 -34.0464 0.0000 -0.0152 -0.0135 41246.8031 0.0000 96.8676 
9   -0.0142 0.0004 -33.6629 0.0000 -0.0150 -0.0134 41245.2609 0.0000 96.8675 
10  -0.0147 0.0004 -34.8368 0.0000 -0.0155 -0.0139 41122.8168 0.0000 96.8582 
        H2 
1  31.9392 
2  31.9387 
3  31.9399 
4  31.9400 
5  31.9329 
6  31.9044 
7  31.9384 
8  31.9248 
9  31.9236 
10 31.8288 

For our proposed multilevel model, no existing packages or functions are ready to use. But you can modify the following code to implement leave-one-out analysis in the framework of multilevel model. Later, we will write a corresponding function and wrap it into orchaRd package. The syntax for using study as unit deleted at a time is:

dat2_Midolo_2019$leave_out <- as.factor(dat2_Midolo_2019$study_name) # create the variable that will be left out
leave1out_ES <- list() # create a list to contain model estimates
for(i in 1:length(levels(dat2_Midolo_2019$leave_out))){
  # create the data with one study removed at a time
  dat <- dat2_Midolo_2019 %>% filter(leave_out != levels(dat2_Midolo_2019$leave_out)[i]) # repeatedly run the multilevel model, leaving out one study at a time
  
  VCV_leave1out <- list() 
  VCV_leave1out[[i]] <- impute_covariance_matrix(vi = dat$lnRRV, cluster = dat$Study_ID, r = 0.5)  # create a list of VCV matrices for following model fitting
  
  leave1out_mod[[i]] <- rma.mv(yi = lnRR, 
                                 V = VCV_leave1out[[i]], 
                                 random = list(~1 | Study_ID,
                                               ~1| ES_ID), 
                                 method = "REML", 
                                 test = "t",
                                 data = dat,
                                 sparse = TRUE)} # model fitting


est.func <- function(mod){
  df <- data.frame(est = mod$b, lower = mod$ci.lb, upper = mod$ci.ub)
  return(df)
} # write a simple function to extract intercept and 95% CIs from each of the above fitted models

leave1out_results <- lapply(leave1out_mod, function(x) est.func(x)) %>% bind_rows %>% mutate(left_out = levels(dat2_Midolo_2019$leave_out)) # turn the model estimates in the format of data frame

The results of leave-one-out analysis are shown on Figure S8, where do not observe any outliers or influential studies obviously having disproportionate effects on the model estimates.

Figure S8 Leave-one-out analysis showing overall effects $\beta_0$ and 95% CIs based on dataset with one study left out at a time from model fitting.

The distribution of overall effect ($\beta_0$) after deleting one study at a time also confirms no obvious outlier (Figure S9).

Figure S9 Leave-one-out analysis showing distribution of overall effects $\beta_0$ based on dataset with one study left out at a time from model fitting.

There are also other diagnostic indices can be calculated from leave-one-out analysis, for example, Cook’s distance, the difference between the regression coefficients (DFBETAS), and the hat value:

• Cook’s distance = the Mahalanobis distance between the predicted (average) effect based on the full dataset and those based on the dataset with one study excluded from the model fitting. A study may be considered as ‘outlier’ if the lower-tail area of a chi-square distribution (with $p$ degrees = the number of model coefficients of freedom) cut off by the Cook’s distance is more than 50%.

• DFBETAS = changes in the predicted (average) effect for a specific study (quantified as standard deviations) after excluding the specific study from the model fitting. A study may be considered as ‘outlier’ the value of is more than 1.

• Hat value = the diagonal elements of the hat matrix that can transform the vector of observed effect into the vector of predicted. Halt value can determine the magnitude of residual based on dataset with a deleted study and therefore can help identify outlying studies. A study may be considered as ‘outlier’ if the Halt value is more than $3(\frac {p} {N})$ (where $N$ = the number of studies).

There is point worth noting - these indices can be easily computed in R, but the above cut-offs are arbitrary. Therefore, if you want to use these indices, you really need informed judgments. Below we show how to calculate these indices.

Cook’s distance can be calculated via cooks.distance():

cooks.distance(mod_ML_lnRR_VCV, cluster = Study_ID)

Figure S10 Cook’s distance showing how much all of the predicted effects in the model change when one study is deleted.

The calculation of DFBETAS and halt values can be done with dfbetas() and hatvalues(), respectively. The syntax is just as simple as that in cooks.distance(). Given the extensive computation these methods will take (our example dataset is very large), we do not show the results.

• Quantitative critical appraisal

Besides the qualitative methods to critically appraise the risk of bias of primary studies included in a meta-analysis, we recommend the quantitative methods. When compiling dataset, you need to extract and code variables related risk-of-bias (RoB) characteristics, for example, the employment of blinding, randomization, and selective reporting. Then you can code these variables as moderators to be incorporated into a meta-regression. Unfortunately, at present, no environmental meta-analyses (at least our surveyed papers) have coded relevant moderators. Therefore, we are unable to provide an example to show implementation. The implementation is straightforward - you just need to use mods argument in rma.mv() function and supply RoB variables to the mods using a formula: mods = ~ RoB.

Our team (Shinichi Nakagawa) has published a paper that examined whether the selective reporting inflates the effect size estimates:

Parker T H, Greig E I, Nakagawa S, et al. Subspecies status and methods explain strength of response to local versus foreign song by oscine birds in meta-analysis[J]. Animal Behaviour, 2018, 142: 1-17.

In this paper, we ran a meta-regression with selective reporting as a moderator and found that selective reporting inflated the effect size estimates. This paper exactly demonstrates the need for transparent reporting & open archiving that we contend in our main text. We refer to the relevant code archived at Open Science Framework to reproduce this quantitative critical appraisal if you think you are capable of conducting such an analysis with our illustration.

missing_SD <- dat2_Midolo_2019 # copy the dataset into missing_SD, which is used as the dataset for illustration of imputation method

set.seed(2022) # set the seed (for the random number generator) to make the following results fully reproducible; it is good practice to set the seed to make our results fully reproducible

stdies <- sample(unique(missing_SD$study_name), size = 0.2*(length(unique(missing_SD$study_name)))) # randomly sample 20% of studies
    
missing_SD[which(missing_SD$study_name %in% stdies), c("sd_treatment", "sd_control")] <- NA # create missingness of SD at the study level

Summary of missingness:

        COLUMN PERCENT_MISSINGNESS IMPUTATIONS
      Study_ID                   0           0
    study_name                   0           0
       species                   0           0
         trait                   0           0
     treatment                   0           0
       control                   0           0
  sd_treatment                  12         159
    sd_control                  12         159
   n_treatment                   0           0
     n_control                   0           0
     elevation                   0           0
 elevation_log                   0           0
          lnRR                   0           0
         lnRRV                   0           0
         ES_ID                   0           0
       ess.var                   0           0
        ess.se                   0           0
       n_tilde                   0           0
          Year                   0           0
        Year.c                   0           0
     leave_out                   0           0

Total missingness: 1.2% (318 imputations needed)

                     COLUMN PERCENT_MISSINGNESS IMPUTATIONS
Study_ID           Study_ID             0.00000           0
study_name       study_name             0.00000           0
species             species             0.00000           0
trait                 trait             0.00000           0
treatment         treatment             0.00000           0
control             control             0.00000           0
sd_treatment   sd_treatment            12.28748         159
sd_control       sd_control            12.28748         159
n_treatment     n_treatment             0.00000           0
n_control         n_control             0.00000           0
elevation         elevation             0.00000           0
elevation_log elevation_log             0.00000           0
lnRR                   lnRR             0.00000           0
lnRRV                 lnRRV             0.00000           0
ES_ID                 ES_ID             0.00000           0
ess.var             ess.var             0.00000           0
ess.se               ess.se             0.00000           0
n_tilde             n_tilde             0.00000           0
Year                   Year             0.00000           0
Year.c               Year.c             0.00000           0
leave_out         leave_out             0.00000           0

dat_simple <- impute_SD(missing_SD, 
                        columnSDnames = c("sd_treatment", "sd_control"), # a string or list containing the labels of the column(s) with missing SD;
                        columnXnames = c("treatment", "control"), a string or list containing the labels of the column(s) with mean values for each SD.
                        method = "HotDeck")

Summary of missingness:

        COLUMN PERCENT_MISSINGNESS IMPUTATIONS
      Study_ID                   0           0
    study_name                   0           0
       species                   0           0
         trait                   0           0
     treatment                   0           0
       control                   0           0
  sd_treatment                   0           0
    sd_control                   0           0
   n_treatment                   0           0
     n_control                   0           0
     elevation                   0           0
 elevation_log                   0           0
          lnRR                   0           0
         lnRRV                   0           0
         ES_ID                   0           0
       ess.var                   0           0
        ess.se                   0           0
       n_tilde                   0           0
          Year                   0           0
        Year.c                   0           0
     leave_out                   0           0

Total missingness: 0% (0 imputations needed)

                     COLUMN PERCENT_MISSINGNESS IMPUTATIONS
Study_ID           Study_ID                   0           0
study_name       study_name                   0           0
species             species                   0           0
trait                 trait                   0           0
treatment         treatment                   0           0
control             control                   0           0
sd_treatment   sd_treatment                   0           0
sd_control       sd_control                   0           0
n_treatment     n_treatment                   0           0
n_control         n_control                   0           0
elevation         elevation                   0           0
elevation_log elevation_log                   0           0
lnRR                   lnRR                   0           0
lnRRV                 lnRRV                   0           0
ES_ID                 ES_ID                   0           0
ess.var             ess.var                   0           0
ess.se               ess.se                   0           0
n_tilde             n_tilde                   0           0
Year                   Year                   0           0
Year.c               Year.c                   0           0
leave_out         leave_out                   0           0

dat_improved <- missing_SD %>% 
                mutate(cv_control = na_if(sd_control / control, Inf),
                       cv_treatment = na_if(sd_treatment / treatment, Inf)) # first calculate CV on dataset with missing SDs wherein missing SD will be ignored when calculating

dat_improved <- cv_avg(x = control, sd = sd_control, n = n_control, 
                       group = study_name, 
                       label = "1", # control group
                       data = dat_improved) # calculate the average between-study CV, which will replace missing SD
    
dat_improved <- cv_avg(x = treatment, sd = sd_treatment, n = n_treatment,
                       group = study_name, 
                       label = "2", # treatment group
                       data = dat_improved)

dat_improved <- dat_improved %>%
                mutate(cv2_control_new = if_else(is.na(cv_control), b_CV2_1, cv_control^2),
                       cv2_treatment_new = if_else(is.na(cv_treatment), b_CV2_2, cv_treatment^2)) # use weighted between-study CV to replace missing CV (which is due to missing SD)


dat_improved <- dat_improved %>%
                mutate(lnRR_new = lnrr_laj(m1 = control, m2 = treatment, cv1_2 = cv2_control_new, cv2_2 = cv2_treatment_new, n1 = n_control, n2 = n_treatment),
                       lnRRV_new = v_lnrr_laj(cv1_2 = cv2_control_new, n1= n_control, cv2_2 = cv2_treatment_new, n2 = n_treatment)) # compute the new point estimate of lnRR and and its sampling variance, respectively. This uses either the between-individual CV^2 when missing or normal CV^2 when not missing

• Multiple sources of non-independence

Environmental meta-analytic datasets may have more many ways for studies to be non-independent from one another. Several ‘clustering variables’, other than study identity, may identify studies as dependent. As you could probably imagine, genus (species) could be such a clustering variable in an environmental dataset with multiple species included, as effect sizes derived from the same species are probably more similar to each other than effect sizes from different species due to similar genetics, shared history or phylogenetic relatedness. The multilevel model has a flexible random effects structure handling complex non-independence like, for example, taxonomic dependence, by adding associated clustering variables (e.g.,species ) as different levels of random effects:

$$z_{i} = \beta_{0} + a_{k[i]} + s_{k[i]} + \mu_{j[i]} + e_{i} + m_{i}, (19)$$ Argument random in rma.mv have an elegant solution to construct random effects structure. For example, if we want to account for both study and species (taxonomic) effects, we need to supply them via argument random and use list() to bind all these clustering variables together: random = list(~ 1 | species, ~ 1 | Study_ID, ~ 1 | ES_ID). The complete code will be:

rma.mv(yi = lnRR, 
       V = VCV, 
       random = list(~1 | species, # add species identity as a random effect, which allows effect sizes vary between species; 
                     ~1 | Study_ID, # add study identity as a random effect, which allows effect sizes vary between studies; 
                     ~1 | ES_ID), # add effect size as a random effect, which allows effect sizes vary within studies. 
       method = "REML", 
       test = "t", 
       data = dat2_Midolo_2019
      )

The output is:

Multivariate Meta-Analysis Model (k = 1294; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-107.3514   214.7028   222.7028   243.3617   222.7338   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0064    109     no   species 
sigma^2.2  0.0113  0.1064     71     no  Study_ID 
sigma^2.3  0.0516  0.2272   1294     no     ES_ID 

Test for Heterogeneity:
Q(df = 1293) = 65913.4258, p-val < .0001

Model Results:

estimate      se    tval    df    pval    ci.lb   ci.ub    
  0.0268  0.0170  1.5730  1293  0.1160  -0.0066  0.0601    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We see that results given under Variance Components show variance components for each of 3 random effects, although there seems not heterogeneity between species.

Note

According to our survey, only a few (XX%) environmental meta-analyses used multilevel model. Therefore, environmental meta-analysts may be not familiar with ‘random effects’? Many researchers have attempted to define it ( see Andrew Gelman’s nice summary ). Here, we give an intuitive explanation in the context of a environmental meta-analysis:

When a variable in a meta-analytic model is specified as a random effect, we believe that it contributes noise (variation) to the overall mean and thus has a random effect on the overall mean. For example, treating ‘species’ as a random effect assumes that the true effects are heterogeneous within species and enables us to quantify how much variation there is among species. On the other hand, treating species as a fixed effect means that species levels are identical across different studies and thus have a systematic effect on the overall effect; this is equivalent to ask: do one species responds more to an environmental stressor than others?).

• select a ‘best’ random effects structure

Here, we introduce some rules of thumb to help decide how to specify vrandom effects structures to capture the hierarchical data structure (e.g., non-independence due to ‘clustering’, ‘nesting’, or ‘crossing’), which may be relevant to your meta-analysis.

Rule 1

Consider whether the variable in question is a true source of heterogeneity according to your knowledge in your field - we are against solely using data-driven random effects structure without considering whether the tested variable has true heterogeneity.

Rule 2

A random effect should have > 5 levels, such that the variance to obtain approximately unbiased estimates.

Rule 3

Examining whether specifying the random effect structure improves model for, such as by examining fit statistics (e.g., AIC and BIC), goodness-of-fit (e.g., $R^2_{marginal}$).

Study (Study_ID) and effect size identities (ES_ID) are typical random effects for an environmental meta-analysis. Let’s use [1]‘s dataset empirically test this point and show how to select ’best’ random effects structure using information-theoretic approaches alongside likelihood ratio tests (implemented via function anova.rma()). There are three random effects candidates:

1. Effect size identity (ES_ID) - unique ID for each pairwise comparison for effect size calculation.

2. Study identity (Study_ID) - unique ID for each included primary study.

3. Species identity (species) - name of species included in the primary studies.

Let’s first fit a null model without including the above random effects candidates as the default reduced model:

null.re <- rma.mv(yi = lnRR, 
                  V = VCV, 
                  method = "ML", # Setting method = "ML" rather than "REML" when conduct model selection.
                  test = "t", 
                  data = dat2_Midolo_2019
                  )

Then, add study identity (Study_ID) as a random effect via argument random:

study <- rma.mv(yi = lnRR, 
                V = VCV, 
                random = ~1 | Study_ID,  
                method = "ML", 
                test = "t", 
                data = dat2_Midolo_2019
                )

Let’s conduct a likelihood ratio test to compare the two fitted models’ quality in terms of AIC, BIC, AICc, log-likelihood values:

anova.rma(null.re, study)

        df        AIC        BIC       AICc      logLik        LRT   pval 
Full     2 47642.9931 47653.3240 47643.0024 -23819.4965                   
Reduced  1 60148.2005 60153.3660 60148.2036 -30073.1003 12507.2075 <.0001 
                QE 
Full    65913.4258 
Reduced 65913.4258 

We see that model (Full) with Study_ID as a random effect has a lower AIC value, in contrast to the null model (Reduced). The log-likelihood ratio test indicates that Study_ID significantly improve model fit (pval = < 0.0001). Note that when we recommend to use maximum likelihood (ML) rather than restricted maximum likelihood (REML) when using information criteria, as using likelihood-based methods (and hence information criteria) to compare models having different fixed effects that are fitted by REML will generally yield nonsense. This is an important point when doing model selection (see next section). It is therefore important to remember to set methods = "ML" rather than methods = "REML" when comparing models with different random effects (and also fixed effects; see next section).

In the same vein, we can examine whether ES_ID contributes to model fit:

es <- rma.mv(yi = lnRR, 
             V = VCV, 
             random = ~1 | ES_ID,  
             method = "ML", 
             test = "t", 
             data = dat2_Midolo_2019
             )

anova.rma(null.re, es)

        df        AIC        BIC       AICc      logLik        LRT   pval 
Full     2   270.8693   281.2003   270.8786   -133.4347                   
Reduced  1 60148.2005 60153.3660 60148.2036 -30073.1003 59879.3312 <.0001 
                QE 
Full    65913.4258 
Reduced 65913.4258 

From AIC, LRT and pval given in the above output, we see that model with ES_ID as a random effect (Full) is much better than model with without random effect (Reduced). Therefore, effect size identity is an important random effect that our model should account for.

Next, lets compare whether a model with both Study_ID and ES_ID as random effects is better than a model with only ES_ID as a random effect:

study.es <- rma.mv(yi = lnRR, 
                   V = VCV, 
                   random = list(~1 | Study_ID,  
                                 ~1 | ES_ID),  
                   method = "ML", 
                   test = "t", 
                   data = dat2_Midolo_2019
                   )

anova.rma(study.es, es)

        df      AIC      BIC     AICc    logLik     LRT   pval         QE 
Full     3 221.5487 237.0451 221.5673 -107.7743                65913.4258 
Reduced  2 270.8693 281.2003 270.8786 -133.4347 51.3207 <.0001 65913.4258 

As expected, the full model has a smaller AIC value than that of reduced model, indicating that model defining the nested random effects (i.e., accounting for non-independence) structure using Study_ID and ES_ID has a better model fit (Full), compared to model ignoring non-independence (Reduced). Next, lets’ explore whether animal species identity is an important random effect. First add the coded variable species as a random effect term via argument random:

species.study.es <- rma.mv(yi = lnRR, 
                              V = VCV, 
                              random = list(~1 | species,
                                            ~1 | Study_ID,  
                                            ~1 | ES_ID),  
                             method = "ML", 
                             test = "t", 
                             data = dat2_Midolo_2019
                             )

anova.rma(species.study.es,study.es)

        df      AIC      BIC     AICc    logLik    LRT   pval         QE 
Full     4 223.5484 244.2103 223.5794 -107.7742               65913.4258 
Reduced  3 221.5487 237.0451 221.5673 -107.7743 0.0003 0.9864 65913.4258 

Looking at the above output, we see that adding species as a random effect does not change the AIC value. This suggests that intraspecific leaf traits are consistent across plant species; that there is only a small amount of heterogeneity among species. This is easily corroborated when calculating $I^2$ at species level:

   I2_Total  I2_species I2_Study_ID    I2_ES_ID 
99.63048070  0.04486139 17.45782398 82.12779532 

As said in our main text, model selection is a powerful method to:

• quantify the importance of moderators in explaining heterogeneity, which is useful when look for, for example, global drivers of environmental changes;

• multimodel inference, which can make model inferences about the moderators in the context of models with all possible combinations of moderators rather than a single ‘best’ model;

• multimodel predictions, which can predict (average) effects of a moderator and its CIs based on models with all possible combinations of moderators rather than a single ‘best’ model.

Below we show how to conduct model selection and multimodel inference using an information-theoretic approach (although $R^2$ based model selection is a reasonable option, it is not preferable here). To do so, metafor package needs to borrow functionality from model-selection-dedicated packages like MuMIn and glmulti. We will use MuMIn package for illustration because the corresponding syntax is more straightforward and simpler than that of glmulti. Let’s stay with the dataset of [1]. This dataset only has two moderators, elevation and trait, so we include an additional term - their interaction - for illustration of a more complex model. We would like to emphasize one point again - in your meta-analysis, you need to need to select your moderators based on their a priori plausibility (predefined environmental questions you are going to address in your analyses) rather than including many more moderators or removing them until you get a significant model.

• Selecting a ‘best5’ model

By ‘best model’, we mean the acceptable amount of information loss when we use a fitted model to approximate the real data generating mechanism. To select the best model, first, we need fit a multilevel multi-moderator meta-analytic model with all plausible moderators (Equation 19; full model). Then, dredge the full model to produce models with all possible combinations of moderators from to the full model (note that you need to use maximum likelihood (ML) rather than restricted ML; see early section for explanation):

mod.full <- rma.mv(yi = lnRR, 
                   V = lnRRV, 
                   mods = ~ trait * elevation_log,
                   random = list(~ 1 | Study_ID, ~ 1 | ES_ID), 
                   method="ML", 
                   test = "t",
                   data = dat2_Midolo_2019) # fit Equation 19 - multilevel multi-moderator meta-analytic model (full model with all plausible moderators).

eval(metafor:::.MuMIn) # use eval() function to extract helper functions from MuMIn and make them usable in metafor.

mod.candidate <- dredge(mod.full, beta = "none", evaluate = TRUE, rank = "AICc", trace=2) # dredge to produce all possible models

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |=============================================================         |  88%

From Table S7, We see the fit statistics (log likelihood and AICc) from all possible models.

Table S7 Fit statistics of models with different combinations of moderators included in the full model (Equation 19 - multilevel multi-moderator meta-analytic model).

Next step, let’s select a sets of best model. This can be done by using thumb of rules, for example, delta AICc <= 2:

subset(mod.candidate, delta <= 2, recalc.weights = FALSE)

Global model call: rma.mv(yi = lnRR, V = lnRRV, mods = ~trait * elevation_log, random = list(~1 | 
    Study_ID, ~1 | ES_ID), data = dat2_Midolo_2019, method = "ML", 
    test = "t", sparse = TRUE)
---
Model selection table 
  (Int)  elv_log trt elv_log:trt df  logLik  AICc delta weight
8     + -0.01991   +           + 16 -74.686 181.8     0      1
Models ranked by AICc(x) 

• Model average

We also can make multi-model inference the model average. You might still remember that in the early section, all model inference (test of model slopes, e.g., $\beta_1$) base the classic null hypothesis testing. The advantage of model averaging is that it can take into account the “weights” from multiple plausible models. “Weights” is the “Akaike weights” that can be considered as the possibility that a specific model has the least information loss (see the above explanation). Put differently, you can treat weights as model probabilities. Note that this dataset is not a good example to illustrate model average because the weight in the first best model is almost 1 and that in other best models are 0 (look at the last column of Table S6). Unfortunately, we have very limited datasets from our survey due to the low data sharing rate in the environmental sciences. Anyway, the logic and syntax are the same:

summary(model.avg(mod.candidate))

Call:
model.avg(object = mod.candidate)

Component model call: 
rma.mv(yi = lnRR, V = lnRRV, mods = ~<5 unique rhs>, random = list(~1 | 
     Study_ID, ~1 | ES_ID), data = dat2_Midolo_2019, method = ML, test = t, 
     sparse = TRUE)

Component models: 
       df  logLik   AICc  delta weight
123    16  -74.69 181.80   0.00      1
12     10  -92.66 205.50  23.70      0
2       9  -95.86 209.86  28.06      0
1       4 -135.48 279.00  97.20      0
(Null)  3 -138.73 283.48 101.68      0

Term codes: 
      elevation_log               trait elevation_log:trait 
                  1                   2                   3 

Model-averaged coefficients:  
(full average) 
                          Estimate Std. Error z value Pr(>|z|)    
intrcpt                   0.115566   0.141088   0.819 0.412727    
elevation_log            -0.019907   0.022636   0.879 0.379164    
traitLA                  -0.097898   0.204791   0.478 0.632624    
traitLMA                 -0.353229   0.188036   1.879 0.060310 .  
traitNarea               -0.468677   0.185908   2.521 0.011701 *  
traitNmass               -0.244480   0.165871   1.474 0.140505    
traitPmass               -0.888061   0.248237   3.577 0.000347 ***
traitSLA                  0.131673   0.198118   0.665 0.506294    
elevation_log:traitLA     0.003795   0.033306   0.114 0.909281    
elevation_log:traitLMA    0.066417   0.030266   2.194 0.028204 *  
elevation_log:traitNarea  0.088490   0.029835   2.966 0.003017 ** 
elevation_log:traitNmass  0.055879   0.026963   2.072 0.038222 *  
elevation_log:traitPmass  0.182724   0.041437   4.410 1.04e-05 ***
elevation_log:traitSLA   -0.025019   0.032394   0.772 0.439920    
 
(conditional average) 
                          Estimate Std. Error z value Pr(>|z|)    
intrcpt                   0.115566   0.141088   0.819 0.412727    
elevation_log            -0.019907   0.022636   0.879 0.379163    
traitLA                  -0.097898   0.204791   0.478 0.632624    
traitLMA                 -0.353229   0.188036   1.879 0.060310 .  
traitNarea               -0.468677   0.185908   2.521 0.011701 *  
traitNmass               -0.244480   0.165871   1.474 0.140505    
traitPmass               -0.888061   0.248237   3.577 0.000347 ***
traitSLA                  0.131673   0.198118   0.665 0.506294    
elevation_log:traitLA     0.003795   0.033306   0.114 0.909281    
elevation_log:traitLMA    0.066417   0.030266   2.194 0.028201 *  
elevation_log:traitNarea  0.088491   0.029834   2.966 0.003016 ** 
elevation_log:traitNmass  0.055879   0.026962   2.073 0.038218 *  
elevation_log:traitPmass  0.182726   0.041434   4.410 1.03e-05 ***
elevation_log:traitSLA   -0.025019   0.032394   0.772 0.439918    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The above output is showing all the estimates of model coefficients for moderators based on model average procedure. It is a typical format to show model estimates in R. The interpretation is quite easy. We also make a table (Table S8) showing results of model average based on glmulti for comparison. We see that the results are exactly the same.

Table S8 The estimates of model slopes for moderators based on the model average inference method.

Scale dependence is a widespread issue in the field of environmental sciences. Yet, it rarely has been accounted for in the current practice of environmental meta-analysis. A common type of scale dependence arises when the sampling unit varies in size, for example a quadrat or plot used to measure biodiversity can vary in size (e.g., 100 $cm^2$ or 1 $km^2$). Smaller plot sizes underestimate biodiversity differences, and effect size statistics that estimate them. However, studies that use smaller plots tend to have greater replication than studies using larger plots. The scale dependence of any meta-analytic inference arises due to a relationship among replicates (sample size), plot size and sampling variance $\nu_i$, which leads to covariation among the three. Basically, scale dependence will lead to biased estimates of sampling variance, which violates a basic assumption of meta-analysis; sampling variance is known because they (partially) serve as weights to average effect sizes across studies.

There are no established solutions on this issue. Incidentally, our team (Rebecca Spake) has a comprehensive simulation paper investigating the scale dependence in the context of meta-analytic model:

Spake R, Mori A S, Beckmann M, et al. Implications of scale dependence for cross‐study syntheses of biodiversity differences[J]. Ecology Letters, 2021, 24(2): 374-390.

Based on the results of this paper and our experience, we give three practical and easy-to-implement suggestions to mitigate the issue of scale dependence:

• 1. use lnRR rather than SMD as the effect size statistic to increase the accuracy of effect size and sampling variance estimates in terms of measurement of biodiversity, richness rather or species density. Visually assess relationships among effect size estimates, plot size, replication and variance.
• 2. conduct sensitivity analysis via, for example, fitting an unweighted model or using a different weighting scheme (e.g., weighting by an ordinal classification of study quality or non-parametric weighting) to get rid of the impact of inaccurate sampling variances; if using unweighted model, your dataset should be free of publication bias, otherwise you will get biased model coefficients
• 3. code the size of plot as a moderator and use meta-regression to control for its covariate effects on other moderators effects

Unfortunately, our survey did not find any dataset have coded plot size as a moderator (Plot.Size). Therefore, we are unable to provide a corresponding illustration. But remember the implementation is straightforward - you just need to use mods argument in rma.mv() function and supply Plot.Size variable to the mods using a formula: mods = ~ Plot.Size.