Animal dispersal costs are not universal

Supplmentary Material

Author

Martinig, A. R., S. L. P. Burk, Y. Yang, M. Lagisz, and S. Nakagawa

Published

June 11, 2025

1 Setting up

1.1 Load packages

Code
#force install pacman and orchaRd

#install.packages("pacman")
#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)
#devtools::install_github("daniel1noble/orchaRd", force = TRUE)

#load all other packages
pacman::p_load(tidyverse, # tidy family and related pacakges below
               kableExtra, # nice tables
               MuMIn,  # multi-model inference
               emmeans, # post-hoc tests
               gridExtra, # may not use this
               pander,   # nice tables
               metafor,  # package for meta-analysis
               ape,      # phylogenetic analysis
               ggtree,
               ggstance,
               MuMIn,  # multi-model inference
               patchwork,   # putting ggplots together - you need to install via devtool
               here,         # making reading files easy
               orchaRd, # plotting orchard plots
               matrixcalc, # matrix calculations
               cowplot #multipanel plots
)

eval(metafor:::.MuMIn)

1.2 Load data: paper and tree data

Code
#the dataset with effect sizes
dat <- read.csv(here("data", "clean_data.csv"))

#dim(dat)
#head(dat)

#the phylogenetic tree
tree <- read.tree(here("data", "species_tree.tre"))

1.3 Calculate effect sizes and sampling variances using custom functions

Code
# function to calculate effect sizes
# Zr - correlation
# there is always n

effect_size <- function(m1, m2, sd1, sd2, n1, n2, n, # 12 arguments
                        est , se, p_val, direction, method){

  if(method == "mean_method"){
  
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s_pool <- sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r    
    
  }else if(method == "count_method"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- sqrt(m1)
    s2 <- sqrt(m2)
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r     
    
  }else if(method == "percent_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "percent_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    sd1 <- sd1/100
    sd2 <- sd2/100
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1))
    m2 <- asin(sqrt(m2))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "t_method1"){
    

    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "t_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "F_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
  
  }else if(method == "F_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "correlation_method1"){
    
    r <- est
    
  }else if(method == "correlation_method2"){
    
    r <- 2*sin((pi/6)*est)
    
  }else if(method == "estimate_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  }else if(method == "estimate_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  } 
  
    if(r >= 1){
    # if over 1, we use 0.99
    Zr <- atanh(0.99)

    }else if(r <= -1){

    Zr <- atanh(-0.99) # r = correlation

    } else {

    Zr <- atanh(r) # r = correlation

    }
  
  VZr <- 1 /(n - 3)
  
  data.frame(ri = r, yi = Zr , vi = VZr)
  
}


##########

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name","\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[ ,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "normal") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  # creating a function to run models
  run_rma1 <- function(name) {
      
    # variance covarinace matrix for sampling errors
    VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
    #VCV1 <- nearPD(VCV1)$mat
      
    rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }

    run_rma2 <- function(name) {
    
            # variance covarinace matrix for sampling errors
            VCVa <- vcalc(vi = dat2$vi_abs, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
            #VCVa <- nearPD(VCVa)$mat
               
               rma.mv(abs_yi2, V = VCVa,
               mods = ~relevel(dat2[[mod]], ref = name),
                      random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
    }
    
    run_rma3 <- function(name) {
      
              # variance covarinace matrix for sampling errors
                VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
                #VCV1 <- nearPD(VCV1)$mat
               
                rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list( # hetero in relation to mod
                            formula(paste("~", mod, "| effectID")),
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     struct = "DIAG",
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }


# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
# use different specifications for aboslute and hetero models 
if (type == "normal"){

    model_all <- purrr::map(level_names[-length(level_names)], run_rma1)

  }else if(type == "absolute"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma2)
    
  }else if(type == "hetero"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma3)
    
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- purrr::map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- purrr::map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

# absolute value analyses

# folded mean
folded_mu <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_mu
} 

# folded variance
folded_v <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2)
  # adding se to make bigger mean
  fold_v <-fold_se^2
  fold_v
} 

1.4 Preparing final dataset

2 Intercept meta-analytic model

2.1 Phylogenetic multilevel models

Code
# takes a while to run
mod <- rma.mv(yi = yi, 
              V = VCV,
              mod = ~ 1,
              data = dat,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned,
                            ~ 1 | phylogeny),
              R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)

# saving the runs
saveRDS(mod, here("Rdata", "mod.RDS"))
Code
# getting the saved model
mod <- readRDS(here("Rdata", "mod.RDS"))

summary(mod)

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

   logLik   Deviance        AIC        BIC       AICc   
-279.0913   558.1826   568.1826   590.7488   568.2724   

Variance Components:

            estim    sqrt  nlvls  fixed           factor    R 
sigma^2.1  0.1105  0.3324    675     no         effectID   no 
sigma^2.2  0.0110  0.1050    202     no          paperID   no 
sigma^2.3  0.0287  0.1694    146     no  species_cleaned   no 
sigma^2.4  0.0000  0.0000    146     no        phylogeny  yes 

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

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0301  0.0303  -0.9947  674  0.3202  -0.0896  0.0294    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# Phylogeney accounts for nothing (0.00%), so we can take it out. 
round(i2_ml(mod),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.74              71.90               7.17              18.67 
      I2_phylogeny 
              0.00 
Code
pred_mod <- predict.rma(mod)

  estimate <- pred_mod$pred
  lowerCL <- pred_mod$ci.lb
  upperCL <- pred_mod$ci.ub 
  lowerPR <- pred_mod$cr.lb
  upperPR <- pred_mod$cr.ub 
  
  table_mod <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod <- table_mod %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod as RDS
saveRDS(tab_mod, here("Rdata", "tab_mod.RDS"))

Table S1. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for intercept meta-analytic model with all five random effects.

Code
tab_mod <-readRDS(here("Rdata", "tab_mod.RDS"))

tab_mod
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.03 -0.09 0.029 0.32 -0.793 0.733

We remove phylogeny as a random effect from our meta-analytic model because phylogeny accounts for nothing (0%)

Code
mod1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)

summary(mod1)

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

   logLik   Deviance        AIC        BIC       AICc   
-243.2890   486.5780   494.5780   512.7536   494.6359   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2732    696     no         effectID 
sigma^2.2  0.0212  0.1457    206     no          paperID 
sigma^2.3  0.0179  0.1337    148     no  species_cleaned 

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

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0240  0.0299  -0.8029  695  0.4223  -0.0827  0.0347    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(i2_ml(mod1),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.06              63.69              18.11              15.26 
Code
reduced_intercept<-orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)
reduced_intercept

Figure S1. Overall meta-analytic mean effect size (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE).
Code
pred_mod1 <- predict.rma(mod1)

  estimate <- pred_mod1$pred
  lowerCL <- pred_mod1$ci.lb
  upperCL <- pred_mod1$ci.ub 
  lowerPR <- pred_mod1$cr.lb
  upperPR <- pred_mod1$cr.ub 
  
  table_mod1 <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod1$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod1 <- table_mod1 %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod1 as RDS
saveRDS(tab_mod1, here("Rdata", "tab_mod1.RDS"))

Table S2. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for reduced intercept meta-analytic model without phylogeny as a random effect.

Code
tab_mod1 <-readRDS(here("Rdata", "tab_mod1.RDS"))

tab_mod1
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.024 -0.083 0.035 0.422 -0.689 0.641

3 Uni-moderator meta-regression models

For each of our moderators, we run a uni-moderator meta-regression model

Code
mod_class <- rma.mv(yi = yi, V = VCV,
               mod = ~ species_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_class)

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

   logLik   Deviance        AIC        BIC       AICc   
-240.9746   481.9491   499.9491   540.7793   500.2138   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0737  0.2715    696     no         effectID 
sigma^2.2  0.0154  0.1241    206     no          paperID 
sigma^2.3  0.0302  0.1738    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 690) = 0.5620, p-val = 0.7607

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
species_classActinopterygii    0.0239  0.1665   0.1434  690  0.8860  -0.3030 
species_classArachnida         0.0426  0.2711   0.1573  690  0.8751  -0.4897 
species_classAves             -0.0262  0.0361  -0.7254  690  0.4684  -0.0970 
species_classInsecta           0.0971  0.1395   0.6960  690  0.4866  -0.1768 
species_classMammalia         -0.0544  0.0451  -1.2067  690  0.2280  -0.1430 
species_classReptilia          0.1153  0.1440   0.8007  690  0.4236  -0.1674 
                              ci.ub    
species_classActinopterygii  0.3508    
species_classArachnida       0.5749    
species_classAves            0.0447    
species_classInsecta         0.3709    
species_classMammalia        0.0341    
species_classReptilia        0.3980    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_class)
   R2_marginal R2_conditional 
    0.00899031     0.38789761 
Code
orchard_plot(mod_class, mod = "species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S2. Mean fitness effect of different taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_class <- all_models(mod_class,  mod = "species_class", type = "normal")

# saving tab_class as RDS
saveRDS(tab_class, here("Rdata", "tab_class.RDS"))

Table S3. Mean estimates of different taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_class <-readRDS(here("Rdata", "tab_class.RDS"))

tab_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Actinopterygii 0.001 -0.344 0.346 0.996 -0.850 0.851
Arachnida 0.039 -0.531 0.610 0.892 -0.925 1.004
Aves -0.027 -0.101 0.046 0.464 -0.808 0.753
Insecta 0.094 -0.202 0.391 0.533 -0.738 0.926
Mammalia -0.071 -0.163 0.022 0.133 -0.853 0.712
Reptilia 0.067 -0.219 0.353 0.645 -0.761 0.895
Actinopterygii-Arachnida 0.038 -0.626 0.703 0.910 -0.984 1.061
Actinopterygii-Aves -0.028 -0.378 0.321 0.873 -0.881 0.824
Actinopterygii-Insecta 0.093 -0.358 0.545 0.685 -0.806 0.993
Actinopterygii-Mammalia -0.071 -0.423 0.280 0.690 -0.925 0.782
Actinopterygii-Reptilia 0.066 -0.376 0.508 0.769 -0.828 0.961
Arachnida-Aves -0.067 -0.640 0.507 0.819 -1.033 0.899
Arachnida-Insecta 0.055 -0.586 0.696 0.866 -0.953 1.063
Arachnida-Mammalia -0.110 -0.685 0.465 0.708 -1.077 0.857
Arachnida-Reptilia 0.028 -0.607 0.663 0.931 -0.976 1.032
Aves-Insecta 0.122 -0.180 0.424 0.429 -0.712 0.956
Aves-Mammalia -0.043 -0.147 0.061 0.417 -0.827 0.741
Aves-Reptilia 0.095 -0.194 0.384 0.520 -0.735 0.924
Insecta-Mammalia -0.165 -0.470 0.140 0.289 -1.000 0.670
Insecta-Reptilia -0.027 -0.434 0.380 0.896 -0.905 0.850
Mammalia-Reptilia 0.138 -0.152 0.428 0.352 -0.692 0.967
Code
mod_type <- rma.mv(yi = yi,  
               V = VCV,
               mod = ~ study_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_type)

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

   logLik   Deviance        AIC        BIC       AICc   
-243.0019   486.0038   496.0038   518.7161   496.0910   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0745  0.2729    696     no         effectID 
sigma^2.2  0.0205  0.1430    206     no          paperID 
sigma^2.3  0.0198  0.1406    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.3754, p-val = 0.6872

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
study_typenatural        -0.0233  0.0306  -0.7623  694  0.4461  -0.0834  0.0368 
study_typesemi-natural   -0.0412  0.0668  -0.6171  694  0.5374  -0.1723  0.0899 
                          
study_typenatural         
study_typesemi-natural    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_type)
   R2_marginal R2_conditional 
  0.0002521658   0.3509052501 
Code
orchard_plot(mod_type, mod = "study_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S3. Mean fitness effect of natural and semi-natural study types. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_type <- all_models(mod_type,  mod = "study_type", type = "normal")

# saving tab_type as RDS
saveRDS(tab_type, here("Rdata", "tab_type.RDS"))

Table S4. Mean estimates of natural and semi-natural study types with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_type <-readRDS(here("Rdata", "tab_type.RDS"))

tab_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Natural -0.028 -0.089 0.033 0.364 -0.795 0.738
Semi-natural -0.062 -0.196 0.072 0.366 -0.838 0.714
Natural-Semi-natural -0.034 -0.165 0.098 0.616 -0.809 0.742
Code
mod_design <- rma.mv(yi = yi, V = VCV,
               mod = ~ study_design - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_design)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.9601   485.9202   495.9202   518.6325   496.0074   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2730    696     no         effectID 
sigma^2.2  0.0233  0.1526    206     no          paperID 
sigma^2.3  0.0158  0.1256    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.4424, p-val = 0.6427

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
study_designcontrast   -0.0189  0.0310  -0.6084  694  0.5431  -0.0798  0.0420 
study_designgradient   -0.0415  0.0459  -0.9043  694  0.3662  -0.1315  0.0486 
                        
study_designcontrast    
study_designgradient    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_design)
   R2_marginal R2_conditional 
  0.0008110263   0.3444046941 
Code
orchard_plot(mod_design, mod = "study_design", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S4. Mean fitness effect of gradient and contrast study designs. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_design <- all_models(mod_design,  mod = "study_design", type = "normal")

# saving tab_design as RDS
saveRDS(tab_design, here("Rdata", "tab_design.RDS"))

Table S5. Mean estimates of gradient and contrast study designs with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_design <-readRDS(here("Rdata", "tab_design.RDS"))

tab_design
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Contrast 0.007 -0.033 0.047 0.735 -0.575 0.589
Gradient -0.039 -0.106 0.028 0.257 -0.624 0.546
Contrast-Gradient -0.025 -0.107 0.057 0.553 -0.758 0.709
Code
mod_disp <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ dispersal_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.7712   485.5425   497.5425   524.7886   497.6649   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2732    696     no         effectID 
sigma^2.2  0.0230  0.1516    206     no          paperID 
sigma^2.3  0.0164  0.1281    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 0.4132, p-val = 0.7436

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
dispersal_typeBoth       -0.0810  0.0804  -1.0068  693  0.3144  -0.2389  0.0769 
dispersal_typebreeding   -0.0144  0.0424  -0.3392  693  0.7346  -0.0976  0.0689 
dispersal_typenatal      -0.0218  0.0325  -0.6695  693  0.5034  -0.0855  0.0420 
                          
dispersal_typeBoth        
dispersal_typebreeding    
dispersal_typenatal       

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_disp)
   R2_marginal R2_conditional 
   0.001361912    0.346299964 
Code
orchard_plot(mod_disp, mod = "dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S5. Mean fitness effect of natal and breeding dispersal (or both dispersal types grouped). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_disp <- all_models(mod_disp,  mod = "dispersal_type", type = "normal")

# saving tab_disp as RDS
saveRDS(tab_disp, here("Rdata", "tab_disp.RDS"))

Table S6. Mean estimates of natal and breeding dispersal (or both dispersal types combined) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_disp <-readRDS(here("Rdata", "tab_disp.RDS"))

tab_disp
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.092 -0.241 0.057 0.226 -0.867 0.683
Breeding -0.013 -0.100 0.075 0.772 -0.779 0.753
Natal -0.027 -0.092 0.038 0.410 -0.791 0.736
B-Breeding 0.079 -0.084 0.242 0.343 -0.699 0.857
B-Natal 0.065 -0.087 0.216 0.403 -0.711 0.840
Breeding-Natal -0.014 -0.100 0.071 0.742 -0.780 0.751
Code
mod_phase <- rma.mv(yi = yi, V = VCV,
               mod = ~ dispersal_phase - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_phase)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.4666   484.9331   494.9331   517.6455   495.0203   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0743  0.2726    696     no         effectID 
sigma^2.2  0.0216  0.1469    206     no          paperID 
sigma^2.3  0.0180  0.1342    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 1.0905, p-val = 0.3366

Model Results:

                            estimate      se     tval   df    pval    ci.lb 
dispersal_phaseprospection   -0.0652  0.0447  -1.4585  694  0.1452  -0.1529 
dispersal_phasesettlement    -0.0141  0.0310  -0.4538  694  0.6501  -0.0750 
                             ci.ub    
dispersal_phaseprospection  0.0226    
dispersal_phasesettlement   0.0469    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_phase)
   R2_marginal R2_conditional 
   0.003449693    0.349617988 
Code
orchard_plot(mod_phase, mod = "dispersal_phase", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S6. Mean fitness effect of prospection and settlement dispersal phases. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_phase <- all_models(mod_phase,  mod = "dispersal_phase", type = "normal")

# saving tab_phase as RDS
saveRDS(tab_phase, here("Rdata", "tab_phase.RDS"))

Table S7. Mean estimates of prospection and settlement dispersal phases with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_phase <-readRDS(here("Rdata", "tab_phase.RDS"))

tab_phase
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Prospection -0.067 -0.158 0.025 0.152 -0.836 0.703
Settlement -0.022 -0.085 0.041 0.493 -0.788 0.745
Prospection-Settlement 0.045 -0.042 0.132 0.314 -0.724 0.813
Code
mod_sex <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.2724   484.5447   496.5447   523.7909   496.6672   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0741  0.2722    696     no         effectID 
sigma^2.2  0.0239  0.1547    206     no          paperID 
sigma^2.3  0.0156  0.1249    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 1.0504, p-val = 0.3696

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
sexBoth     -0.0702  0.0420  -1.6722  693  0.0949  -0.1526  0.0122  . 
sexFemale   -0.0082  0.0338  -0.2433  693  0.8078  -0.0745  0.0581    
sexMale     -0.0037  0.0349  -0.1058  693  0.9158  -0.0723  0.0649    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex)
   R2_marginal R2_conditional 
   0.005223149    0.351422447 
Code
orchard_plot(mod_sex, mod = "sex", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S7. Mean fitness effect of females, males, or both sexes grouped. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex <- all_models(mod_sex,  mod = "sex", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))

Table S8. Mean estimates of females, males, or both sexes grouped with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex <-readRDS(here("Rdata", "tab_sex.RDS"))

tab_sex
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.073 -0.162 0.015 0.102 -0.841 0.694
F -0.013 -0.083 0.056 0.711 -0.779 0.753
M -0.017 -0.088 0.055 0.650 -0.782 0.749
B-F 0.060 -0.032 0.152 0.199 -0.708 0.828
B-M 0.057 -0.037 0.151 0.237 -0.711 0.825
F-M -0.003 -0.062 0.055 0.908 -0.768 0.761

3.1 Higher life history grouping

Code
mod_age1 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age1)

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

   logLik   Deviance        AIC        BIC       AICc   
-241.5610   483.1219   495.1219   522.3681   495.2444   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0750  0.2738    696     no         effectID 
sigma^2.2  0.0172  0.1312    206     no          paperID 
sigma^2.3  0.0219  0.1479    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 1.4496, p-val = 0.2272

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
age_class_cleanadult      -0.0190  0.0310  -0.6131  693  0.5400  -0.0799 
age_class_cleanjuvenile    0.0059  0.0513   0.1151  693  0.9084  -0.0948 
age_class_cleanmix        -0.1213  0.0607  -1.9999  693  0.0459  -0.2404 
                           ci.ub    
age_class_cleanadult      0.0419    
age_class_cleanjuvenile   0.1066    
age_class_cleanmix       -0.0022  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age1)
   R2_marginal R2_conditional 
   0.006433797    0.346790908 
Code
orchard_plot(mod_age1, mod = "age_class_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S8. Mean fitness effect of life history stage (adult, juvenile, or mixed). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_age1 <- all_models(mod_age1,  mod = "age_class_clean", type = "normal")

# saving tab_age1 as RDS
saveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))

Table S9. Mean estimates of life history stage (adult, juvenile, or mixed) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_age1 <-readRDS(here("Rdata", "tab_age1.RDS"))

tab_age1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Adult -0.004 -0.061 0.052 0.879 -0.735 0.726
Juvenile 0.012 -0.080 0.105 0.794 -0.722 0.747
Mix -0.101 -0.224 0.022 0.108 -0.839 0.638
Adult-Juvenile 0.017 -0.072 0.106 0.713 -0.717 0.751
Adult-Mix -0.096 -0.217 0.025 0.118 -0.835 0.642
Juvenile-Mix -0.113 -0.254 0.027 0.115 -0.855 0.629

3.2 Finer life stage grouping

Code
# age_class

mod_age2 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age2)

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

   logLik   Deviance        AIC        BIC       AICc   
-240.3580   480.7160   500.7160   546.0684   501.0405   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0754  0.2746    696     no         effectID 
sigma^2.2  0.0183  0.1354    206     no          paperID 
sigma^2.3  0.0198  0.1407    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 689) = 0.9165, p-val = 0.4929

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
age_classA     -0.0112  0.0338  -0.3309  689  0.7408  -0.0775  0.0552    
age_classJ      0.0042  0.0515   0.0811  689  0.9354  -0.0969  0.1053    
age_classJA    -0.3023  0.1796  -1.6835  689  0.0927  -0.6549  0.0503  . 
age_classJY    -0.1398  0.0882  -1.5850  689  0.1134  -0.3130  0.0334    
age_classJYA   -0.0662  0.0854  -0.7745  689  0.4389  -0.2340  0.1016    
age_classY     -0.0432  0.0574  -0.7528  689  0.4519  -0.1559  0.0695    
age_classYA    -0.0375  0.0467  -0.8028  689  0.4224  -0.1292  0.0542    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age2)
   R2_marginal R2_conditional 
    0.01178388     0.34371431 
Code
orchard_plot(mod_age2, mod = "age_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S9. Mean fitness effect of life history stage at the level of precision reported (A: adult, Y: yearling, J: juvenile). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)

3.3 Higher level

Code
mod_fit1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_higher_level - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit1)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.1277   484.2555   496.2555   523.5017   496.3779   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0744  0.2728    696     no         effectID 
sigma^2.2  0.0211  0.1453    206     no          paperID 
sigma^2.3  0.0180  0.1342    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 1.2241, p-val = 0.3000

Model Results:

                                      estimate      se     tval   df    pval 
fitness_higher_levelindirect fitness   -0.2389  0.1452  -1.6452  693  0.1004 
fitness_higher_levelreproduction       -0.0103  0.0324  -0.3193  693  0.7496 
fitness_higher_levelsurvival           -0.0338  0.0332  -1.0180  693  0.3090 
                                        ci.lb   ci.ub    
fitness_higher_levelindirect fitness  -0.5239  0.0462    
fitness_higher_levelreproduction      -0.0740  0.0533    
fitness_higher_levelsurvival          -0.0990  0.0314    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit1)
   R2_marginal R2_conditional 
   0.004159334    0.347294612 
Code
orchard_plot(mod_fit1, mod = "fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S10. Mean fitness effect of different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_fit1 <- all_models(mod_fit1,  mod = "fitness_higher_level", type = "normal")

# saving tab_fit1 as RDS
saveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))

Table S10. Mean estimates of different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_fit1 <-readRDS(here("Rdata", "tab_fit1.RDS"))

tab_fit1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Indirect fitness -0.237 -0.572 0.098 0.166 -1.068 0.595
Reproduction -0.019 -0.084 0.047 0.580 -0.782 0.745
Survival -0.038 -0.106 0.030 0.273 -0.802 0.726
Indirect fitness-Reproduction 0.218 -0.117 0.554 0.202 -0.613 1.050
Indirect fitness-Survival 0.199 -0.136 0.534 0.244 -0.632 1.030
Reproduction-Survival -0.019 -0.079 0.040 0.525 -0.783 0.744

3.4 Finer level

Code
mod_fit2 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_metric_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit2)

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

   logLik   Deviance        AIC        BIC       AICc   
-241.0975   482.1951   504.1951   554.0667   504.5856   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0740  0.2721    696     no         effectID 
sigma^2.2  0.0166  0.1287    206     no          paperID 
sigma^2.3  0.0257  0.1602    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 688) = 0.9608, p-val = 0.4656

Model Results:

                                                           estimate      se 
fitness_metric_cleanage at first reproduction               -0.0481  0.0922 
fitness_metric_cleanindirect fitness                        -0.2312  0.1459 
fitness_metric_cleanlifetime breeding success                0.0043  0.0372 
fitness_metric_cleanlifetime reproductive success           -0.0264  0.0400 
fitness_metric_cleanoffspring reproduction                   0.0341  0.0877 
fitness_metric_cleanoffspring survival                       0.0105  0.0445 
fitness_metric_cleanreproductive lifespan and/or attempts   -0.0438  0.0899 
fitness_metric_cleansurvival                                -0.0586  0.0370 
                                                              tval   df    pval 
fitness_metric_cleanage at first reproduction              -0.5218  688  0.6020 
fitness_metric_cleanindirect fitness                       -1.5847  688  0.1135 
fitness_metric_cleanlifetime breeding success               0.1166  688  0.9072 
fitness_metric_cleanlifetime reproductive success          -0.6598  688  0.5096 
fitness_metric_cleanoffspring reproduction                  0.3891  688  0.6974 
fitness_metric_cleanoffspring survival                      0.2359  688  0.8136 
fitness_metric_cleanreproductive lifespan and/or attempts  -0.4874  688  0.6262 
fitness_metric_cleansurvival                               -1.5850  688  0.1134 
                                                             ci.lb   ci.ub    
fitness_metric_cleanage at first reproduction              -0.2292  0.1329    
fitness_metric_cleanindirect fitness                       -0.5177  0.0553    
fitness_metric_cleanlifetime breeding success              -0.0688  0.0775    
fitness_metric_cleanlifetime reproductive success          -0.1050  0.0522    
fitness_metric_cleanoffspring reproduction                 -0.1381  0.2064    
fitness_metric_cleanoffspring survival                     -0.0768  0.0978    
fitness_metric_cleanreproductive lifespan and/or attempts  -0.2202  0.1326    
fitness_metric_cleansurvival                               -0.1312  0.0140    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit2)
   R2_marginal R2_conditional 
    0.01004347     0.36962120 
Code
orchard_plot(mod_fit2, mod = "fitness_metric_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 0)

Figure S11. Mean fitness effect of different fitness components at a finer resolution. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)
Code
mod_gen <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ whose_fitness - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_gen)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.6778   485.3556   495.3556   518.0679   495.4428   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0740  0.2721    696     no         effectID 
sigma^2.2  0.0179  0.1337    206     no          paperID 
sigma^2.3  0.0240  0.1548    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 1.1603, p-val = 0.3140

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
whose_fitnessdescendant    0.0118  0.0430   0.2753  694  0.7832  -0.0725 
whose_fitnessfocal        -0.0325  0.0310  -1.0486  694  0.2947  -0.0932 
                          ci.ub    
whose_fitnessdescendant  0.0962    
whose_fitnessfocal       0.0283    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_gen)
   R2_marginal R2_conditional 
   0.002332594    0.362620116 
Code
orchard_plot(mod_gen, mod = "whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S12. Mean fitness effect of focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_gen <- all_models(mod_gen,  mod = "whose_fitness", type = "normal")

# saving tab_gen as RDS
saveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))

Table S11. Mean estimates of non-dispersers and dispersers and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_gen <-readRDS(here("Rdata", "tab_gen.RDS"))

tab_gen
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Individual -0.038 -0.100 0.023 0.219 -0.808 0.731
Offspring 0.005 -0.086 0.096 0.907 -0.767 0.777
Individual-Offspring 0.044 -0.035 0.123 0.275 -0.727 0.815

3.5 Normal analysis

Code
# getting absolute values for effect sizes as we expect shorter years will have more fluctuations
dat$ln_study_duration <- log(dat$study_duration+1)

mod_dur <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ ln_study_duration,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur)

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

   logLik   Deviance        AIC        BIC       AICc   
-243.0283   486.0566   496.0566   518.7690   496.1438   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2731    696     no         effectID 
sigma^2.2  0.0215  0.1467    206     no          paperID 
sigma^2.3  0.0182  0.1348    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.0060, p-val = 0.9380

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0202  0.0599  -0.3368  694  0.7364  -0.1378  0.0974    
ln_study_duration   -0.0016  0.0211  -0.0778  694  0.9380  -0.0430  0.0398    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_dur)
   R2_marginal R2_conditional 
  2.035908e-05   3.473406e-01 
Code
bubble_plot(mod_dur,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S13. Mean fitness effect of study duration (log years). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# absolute value analysis
mod_dur_abs <- rma.mv(yi = yi_abs,
                V = VCVa,
               mod = ~ ln_study_duration,
               data = dat,
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur_abs)

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

  logLik  Deviance       AIC       BIC      AICc   
 28.7611  -57.5223  -47.5223  -24.8099  -47.4351   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0326  0.1806    696     no         effectID 
sigma^2.2  0.0081  0.0900    206     no          paperID 
sigma^2.3  0.0090  0.0948    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 3.7270, p-val = 0.0539

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0359  0.0396  -0.9060  694  0.3653  -0.1137  0.0419    
ln_study_duration    0.0268  0.0139   1.9305  694  0.0539  -0.0005  0.0540  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_dur_abs)
   R2_marginal R2_conditional 
    0.01232166     0.35185508 
Code
bubble_plot(mod_dur_abs,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S14. Mean fitness effect of study duration (log years) from absolute value analysis. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

3.6 Normal analysis

Code
mod_focus <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_focus)

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

   logLik   Deviance        AIC        BIC       AICc   
-243.0004   486.0008   496.0008   518.7132   496.0881   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0744  0.2728    696     no         effectID 
sigma^2.2  0.0213  0.1459    206     no          paperID 
sigma^2.3  0.0188  0.1371    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.4197, p-val = 0.6574

Model Results:

                                 estimate      se     tval   df    pval 
fitness_main_focusNon-selective   -0.0097  0.0462  -0.2102  694  0.8336 
fitness_main_focusSelective       -0.0281  0.0313  -0.8982  694  0.3694 
                                   ci.lb   ci.ub    
fitness_main_focusNon-selective  -0.1003  0.0809    
fitness_main_focusSelective      -0.0896  0.0334    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus)
   R2_marginal R2_conditional 
  0.0005343088   0.3503769299 
Code
orchard_plot(mod_focus, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S15. Mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_focus <- all_models(mod_focus,  mod = "fitness_main_focus", type = "normal")

# saving tab_focus as RDS
saveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))

Table S14. Mean estimates of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_focus <-readRDS(here("Rdata", "tab_focus.RDS"))

tab_focus
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
N -0.019 -0.110 0.073 0.690 -0.787 0.750
Y -0.034 -0.096 0.029 0.292 -0.799 0.732
N-Y -0.015 -0.103 0.073 0.734 -0.783 0.752

3.7 Comparing normal model to heteroscedastic model

Code
# homo
mod_focus1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_focus1)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-243.8844  2860.9222   497.7688   520.4956   497.8558   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2731    696     no         effectID 
sigma^2.2  0.0214  0.1464    206     no          paperID 
sigma^2.3  0.0173  0.1314    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.4016, p-val = 0.6694

Model Results:

                                 estimate      se     tval   df    pval 
fitness_main_focusNon-selective   -0.0095  0.0459  -0.2066  694  0.8364 
fitness_main_focusSelective       -0.0273  0.0310  -0.8792  694  0.3796 
                                   ci.lb   ci.ub    
fitness_main_focusNon-selective  -0.0995  0.0806    
fitness_main_focusSelective      -0.0882  0.0336    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus1)
   R2_marginal R2_conditional 
  0.0005038894   0.3419616519 
Code
# hetero
mod_focus2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ fitness_main_focus,
                 data = dat, 
                 random = list(
                   ~ fitness_main_focus | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_focus2)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-238.8875  2850.9283   489.7749   517.0470   489.8968   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0241  0.1553    206     no          paperID 
sigma^2.2  0.0163  0.1277    148     no  species_cleaned 

outer factor: effectID           (nlvls = 696)
inner factor: fitness_main_focus (nlvls = 2)

            estim    sqrt  k.lvl  fixed          level 
tau^2.1    0.0430  0.2073    164     no  Non-selective 
tau^2.2    0.0817  0.2859    532     no      Selective 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.2219, p-val = 0.6377

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
intrcpt                       -0.0057  0.0433  -0.1323  694  0.8948  -0.0908 
fitness_main_focusSelective   -0.0195  0.0415  -0.4711  694  0.6377  -0.1009 
                              ci.ub    
intrcpt                      0.0793    
fitness_main_focusSelective  0.0619    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
#r2_ml(mod_focusA1)

# they are not significantly different
anova(mod_focus1, mod_focus2)

        df      AIC      BIC     AICc    logLik    LRT   pval         QE 
Full     6 489.7749 517.0470 489.8968 -238.8875               15614.7766 
Reduced  5 497.7688 520.4956 497.8558 -243.8844 9.9939 0.0016 15614.7766 
Code
# homo
orchard_plot(mod_focus1, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S16. Maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# hetero
orchard_plot(mod_focus2, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S17. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

3.8 Normal analysis

Code
mod_confirm <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_confirm)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.2655   484.5309   498.5309   530.3080   498.6947   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0743  0.2726    696     no         effectID 
sigma^2.2  0.0275  0.1659    206     no          paperID 
sigma^2.3  0.0122  0.1104    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 692) = 0.3695, p-val = 0.8304

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0300  0.0829   0.3613  692  0.7180  -0.1328 
confirmation_biasCost      -0.0367  0.0389  -0.9428  692  0.3461  -0.1132 
confirmation_biasMixed     -0.0254  0.0394  -0.6460  692  0.5185  -0.1028 
confirmation_biasNeutral   -0.0008  0.0469  -0.0163  692  0.9870  -0.0929 
                           ci.ub    
confirmation_biasBenefit  0.1928    
confirmation_biasCost     0.0397    
confirmation_biasMixed    0.0519    
confirmation_biasNeutral  0.0914    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm)
   R2_marginal R2_conditional 
   0.002617906    0.349759459 
Code
orchard_plot(mod_confirm, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S18. Mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_confirm <- all_models(mod_confirm,  mod = "confirmation_bias", type = "normal")

# saving tab_confirm as RDS
saveRDS(tab_confirm, here("Rdata", "tab_confirm.RDS"))

Table S15. Mean estimates of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_confirm <-readRDS(here("Rdata", "tab_confirm.RDS"))

tab_confirm
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Benefit 0.012 -0.154 0.178 0.884 -0.764 0.788
Cost -0.042 -0.122 0.037 0.296 -0.804 0.720
Mixed -0.028 -0.107 0.050 0.483 -0.790 0.734
Neutral -0.009 -0.102 0.084 0.852 -0.772 0.755
Benefit-Cost -0.055 -0.227 0.118 0.536 -0.832 0.723
Benefit-Mixed -0.040 -0.212 0.132 0.645 -0.818 0.737
Benefit-Neutral -0.021 -0.198 0.156 0.815 -0.800 0.757
Cost-Mixed 0.014 -0.081 0.110 0.770 -0.750 0.778
Cost-Neutral 0.033 -0.071 0.138 0.531 -0.732 0.799
Mixed-Neutral 0.019 -0.085 0.124 0.719 -0.746 0.784

3.9 Comparing normal model to heteroscedastic model

Code
# homo
mod_confirm1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_confirm1)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-243.5314  2860.2161   501.0628   532.8802   501.2255   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2731    696     no         effectID 
sigma^2.2  0.0263  0.1623    206     no          paperID 
sigma^2.3  0.0110  0.1050    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 692) = 0.3600, p-val = 0.8371

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0297  0.0819   0.3623  692  0.7173  -0.1311 
confirmation_biasCost      -0.0358  0.0384  -0.9320  692  0.3516  -0.1112 
confirmation_biasMixed     -0.0241  0.0387  -0.6221  692  0.5341  -0.1001 
confirmation_biasNeutral   -0.0002  0.0463  -0.0046  692  0.9963  -0.0911 
                           ci.ub    
confirmation_biasBenefit  0.1905    
confirmation_biasCost     0.0396    
confirmation_biasMixed    0.0519    
confirmation_biasNeutral  0.0907    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm1)
   R2_marginal R2_conditional 
   0.002571009    0.335545088 
Code
# hetero
mod_confirm2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ confirmation_bias,
                 data = dat, 
                 random = list(
                   ~ confirmation_bias | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_confirm2)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-224.2388  2821.6311   468.4777   513.9312   468.7989   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0288  0.1698    206     no          paperID 
sigma^2.2  0.0052  0.0719    148     no  species_cleaned 

outer factor: effectID          (nlvls = 696)
inner factor: confirmation_bias (nlvls = 4)

            estim    sqrt  k.lvl  fixed    level 
tau^2.1    0.0129  0.1135     32     no  Benefit 
tau^2.2    0.0767  0.2769    210     no     Cost 
tau^2.3    0.0565  0.2377    289     no    Mixed 
tau^2.4    0.1380  0.3715    165     no  Neutral 

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

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 692) = 0.6362, p-val = 0.5919

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
intrcpt                     0.0503  0.0641   0.7840  692  0.4333  -0.0756 
confirmation_biasCost      -0.0921  0.0680  -1.3547  692  0.1760  -0.2255 
confirmation_biasMixed     -0.0779  0.0672  -1.1593  692  0.2467  -0.2099 
confirmation_biasNeutral   -0.0630  0.0736  -0.8552  692  0.3928  -0.2075 
                           ci.ub    
intrcpt                   0.1762    
confirmation_biasCost     0.0414    
confirmation_biasMixed    0.0540    
confirmation_biasNeutral  0.0816    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# they are not significantly different
anova(mod_confirm1, mod_confirm2)

        df      AIC      BIC     AICc    logLik     LRT   pval         QE 
Full    10 468.4777 513.9312 468.7989 -224.2388                15529.2933 
Reduced  7 501.0628 532.8802 501.2255 -243.5314 38.5851 <.0001 15529.2933 
Code
# homo
orchard_plot(mod_confirm1, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S19. Maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# hetero
orchard_plot(mod_confirm2, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S20. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

4 Interaction meta-regression models (2 moderators)

Code
# the interaction between 2 categorical variables are the same as creating a new variable with 2 categories combined so we use that approach 

dat$sex_species_class <- as.factor(paste(dat$sex, dat$species_class ,sep = "_"))

mod_sex_species_class <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_species_class - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_species_class)

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

   logLik   Deviance        AIC        BIC       AICc   
-236.6135   473.2270   509.2270   590.6511   510.2602   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0741  0.2722    696     no         effectID 
sigma^2.2  0.0207  0.1439    206     no          paperID 
sigma^2.3  0.0195  0.1397    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:15):
F(df1 = 15, df2 = 681) = 0.9762, p-val = 0.4787

Model Results:

                                        estimate      se     tval   df    pval 
sex_species_classBoth_Aves               -0.0077  0.0517  -0.1480  681  0.8824 
sex_species_classBoth_Insecta            -0.1445  0.2402  -0.6015  681  0.5477 
sex_species_classBoth_Mammalia           -0.1764  0.0646  -2.7291  681  0.0065 
sex_species_classBoth_Reptilia            0.1653  0.2740   0.6034  681  0.5465 
sex_species_classFemale_Actinopterygii   -0.0592  0.1944  -0.3047  681  0.7607 
sex_species_classFemale_Arachnida         0.0410  0.2612   0.1570  681  0.8753 
sex_species_classFemale_Aves             -0.0326  0.0413  -0.7909  681  0.4293 
sex_species_classFemale_Insecta           0.2482  0.1764   1.4071  681  0.1599 
sex_species_classFemale_Mammalia         -0.0130  0.0495  -0.2632  681  0.7925 
sex_species_classFemale_Reptilia          0.1290  0.1533   0.8415  681  0.4004 
sex_species_classMale_Actinopterygii      0.0948  0.1901   0.4989  681  0.6180 
sex_species_classMale_Aves               -0.0352  0.0416  -0.8466  681  0.3975 
sex_species_classMale_Insecta            -0.1101  0.3352  -0.3285  681  0.7426 
sex_species_classMale_Mammalia            0.0216  0.0540   0.3994  681  0.6897 
sex_species_classMale_Reptilia            0.0132  0.2204   0.0597  681  0.9524 
                                          ci.lb    ci.ub     
sex_species_classBoth_Aves              -0.1091   0.0938     
sex_species_classBoth_Insecta           -0.6160   0.3271     
sex_species_classBoth_Mammalia          -0.3033  -0.0495  ** 
sex_species_classBoth_Reptilia          -0.3727   0.7033     
sex_species_classFemale_Actinopterygii  -0.4409   0.3224     
sex_species_classFemale_Arachnida       -0.4719   0.5539     
sex_species_classFemale_Aves            -0.1136   0.0484     
sex_species_classFemale_Insecta         -0.0981   0.5946     
sex_species_classFemale_Mammalia        -0.1102   0.0841     
sex_species_classFemale_Reptilia        -0.1720   0.4299     
sex_species_classMale_Actinopterygii    -0.2784   0.4680     
sex_species_classMale_Aves              -0.1169   0.0464     
sex_species_classMale_Insecta           -0.7684   0.5481     
sex_species_classMale_Mammalia          -0.0845   0.1277     
sex_species_classMale_Reptilia          -0.4197   0.4460     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_species_class)
   R2_marginal R2_conditional 
     0.0283798      0.3704015 
Code
orchard_plot(mod_sex_species_class, mod = "sex_species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S21. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_species_class <- all_models(mod_sex_species_class,  mod = "sex_species_class", type = "normal")

# saving tab_sex_species_class as RDS
saveRDS(tab_sex_species_class, here("Rdata", "tab_sex_species_class.RDS"))

Table S16. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_species_class <-readRDS(here("Rdata", "tab_sex_species_class.RDS"))

tab_sex_species_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_Aves -0.001 -0.116 0.114 0.987 -0.777 0.775
B_Insecta -0.146 -0.689 0.397 0.597 -1.086 0.794
B_Mammalia -0.180 -0.312 -0.047 0.008 -0.959 0.599
B_Reptilia 0.158 -0.381 0.696 0.565 -0.780 1.095
F_Actinopterygii -0.069 -0.477 0.339 0.741 -0.938 0.800
F_Arachnida 0.040 -0.514 0.594 0.887 -0.906 0.987
F_Aves -0.036 -0.124 0.051 0.417 -0.809 0.736
F_Insecta 0.240 -0.128 0.607 0.201 -0.611 1.091
F_Mammalia -0.019 -0.124 0.087 0.731 -0.793 0.756
F_Reptilia 0.071 -0.238 0.380 0.651 -0.756 0.898
M_Actinopterygii 0.068 -0.339 0.475 0.742 -0.800 0.937
M_Aves -0.035 -0.123 0.054 0.441 -0.807 0.738
M_Insecta -0.143 -0.888 0.602 0.707 -1.212 0.927
M_Mammalia -0.016 -0.129 0.097 0.778 -0.792 0.760
M_Reptilia -0.024 -0.481 0.433 0.918 -0.917 0.869
B_Aves-B_Insecta -0.145 -0.700 0.409 0.607 -1.092 0.801
B_Aves-B_Mammalia -0.179 -0.346 -0.012 0.036 -0.964 0.607
B_Aves-B_Reptilia 0.159 -0.388 0.706 0.569 -0.784 1.101
B_Aves-F_Actinopterygii -0.068 -0.489 0.353 0.752 -0.943 0.807
B_Aves-F_Arachnida 0.041 -0.523 0.605 0.886 -0.912 0.994
B_Aves-F_Aves -0.035 -0.165 0.094 0.594 -0.814 0.743
B_Aves-F_Insecta 0.241 -0.141 0.623 0.217 -0.617 1.098
B_Aves-F_Mammalia -0.018 -0.165 0.130 0.815 -0.799 0.764
B_Aves-F_Reptilia 0.072 -0.253 0.396 0.663 -0.761 0.905
B_Aves-M_Actinopterygii 0.069 -0.351 0.490 0.746 -0.806 0.944
B_Aves-M_Aves -0.034 -0.163 0.096 0.611 -0.812 0.745
B_Aves-M_Insecta -0.142 -0.894 0.611 0.712 -1.217 0.933
B_Aves-M_Mammalia -0.015 -0.169 0.138 0.845 -0.798 0.767
B_Aves-M_Reptilia -0.023 -0.491 0.445 0.923 -0.922 0.876
B_Insecta-B_Mammalia -0.033 -0.591 0.524 0.906 -0.982 0.915
B_Insecta-B_Reptilia 0.304 -0.459 1.067 0.435 -0.779 1.386
B_Insecta-F_Actinopterygii 0.077 -0.601 0.756 0.823 -0.947 1.102
B_Insecta-F_Arachnida 0.186 -0.589 0.962 0.637 -0.904 1.277
B_Insecta-F_Aves 0.110 -0.439 0.659 0.694 -0.834 1.054
B_Insecta-F_Insecta 0.386 -0.269 1.041 0.248 -0.623 1.395
B_Insecta-F_Mammalia 0.128 -0.424 0.680 0.650 -0.818 1.073
B_Insecta-F_Reptilia 0.217 -0.406 0.841 0.494 -0.771 1.206
B_Insecta-M_Actinopterygii 0.215 -0.463 0.892 0.534 -0.809 1.239
B_Insecta-M_Aves 0.112 -0.438 0.661 0.690 -0.832 1.056
B_Insecta-M_Insecta 0.004 -0.918 0.925 0.994 -1.196 1.203
B_Insecta-M_Mammalia 0.130 -0.423 0.684 0.645 -0.816 1.076
B_Insecta-M_Reptilia 0.122 -0.587 0.831 0.735 -0.922 1.167
B_Mammalia-B_Reptilia 0.337 -0.210 0.885 0.227 -0.605 1.280
B_Mammalia-F_Actinopterygii 0.111 -0.312 0.534 0.607 -0.765 0.987
B_Mammalia-F_Arachnida 0.220 -0.347 0.786 0.446 -0.734 1.174
B_Mammalia-F_Aves 0.144 -0.004 0.291 0.056 -0.638 0.925
B_Mammalia-F_Insecta 0.419 0.035 0.804 0.033 -0.439 1.278
B_Mammalia-F_Mammalia 0.161 0.020 0.303 0.026 -0.619 0.942
B_Mammalia-F_Reptilia 0.251 -0.075 0.577 0.132 -0.583 1.085
B_Mammalia-M_Actinopterygii 0.248 -0.175 0.671 0.250 -0.628 1.124
B_Mammalia-M_Aves 0.145 -0.002 0.293 0.054 -0.636 0.927
B_Mammalia-M_Insecta 0.037 -0.717 0.791 0.923 -1.039 1.113
B_Mammalia-M_Mammalia 0.163 0.015 0.312 0.031 -0.618 0.945
B_Mammalia-M_Reptilia 0.156 -0.314 0.625 0.515 -0.744 1.055
B_Reptilia-F_Actinopterygii -0.226 -0.897 0.444 0.507 -1.245 0.793
B_Reptilia-F_Arachnida -0.117 -0.887 0.652 0.764 -1.204 0.969
B_Reptilia-F_Aves -0.194 -0.735 0.347 0.482 -1.133 0.745
B_Reptilia-F_Insecta 0.082 -0.565 0.729 0.803 -0.922 1.086
B_Reptilia-F_Mammalia -0.176 -0.718 0.366 0.524 -1.116 0.764
B_Reptilia-F_Reptilia -0.087 -0.653 0.479 0.764 -1.040 0.867
B_Reptilia-M_Actinopterygii -0.089 -0.760 0.581 0.794 -1.109 0.930
B_Reptilia-M_Aves -0.192 -0.733 0.349 0.486 -1.131 0.747
B_Reptilia-M_Insecta -0.300 -1.217 0.616 0.520 -1.496 0.895
B_Reptilia-M_Mammalia -0.174 -0.718 0.370 0.530 -1.115 0.767
B_Reptilia-M_Reptilia -0.182 -0.860 0.497 0.599 -1.206 0.843
F_Actinopterygii-F_Arachnida 0.109 -0.577 0.795 0.755 -0.920 1.138
F_Actinopterygii-F_Aves 0.033 -0.381 0.446 0.877 -0.839 0.904
F_Actinopterygii-F_Insecta 0.308 -0.237 0.854 0.267 -0.633 1.250
F_Actinopterygii-F_Mammalia 0.050 -0.366 0.466 0.812 -0.823 0.923
F_Actinopterygii-F_Reptilia 0.140 -0.366 0.646 0.588 -0.779 1.059
F_Actinopterygii-M_Actinopterygii 0.137 -0.324 0.598 0.559 -0.758 1.032
F_Actinopterygii-M_Aves 0.034 -0.379 0.448 0.871 -0.838 0.906
F_Actinopterygii-M_Insecta -0.074 -0.921 0.773 0.864 -1.217 1.069
F_Actinopterygii-M_Mammalia 0.053 -0.366 0.471 0.805 -0.821 0.927
F_Actinopterygii-M_Reptilia 0.045 -0.564 0.653 0.885 -0.935 1.024
F_Arachnida-F_Aves -0.076 -0.635 0.483 0.789 -1.026 0.873
F_Arachnida-F_Insecta 0.200 -0.463 0.862 0.555 -0.814 1.213
F_Arachnida-F_Mammalia -0.059 -0.620 0.502 0.837 -1.009 0.892
F_Arachnida-F_Reptilia 0.031 -0.600 0.662 0.924 -0.963 1.024
F_Arachnida-M_Actinopterygii 0.028 -0.657 0.714 0.936 -1.001 1.057
F_Arachnida-M_Aves -0.075 -0.634 0.484 0.793 -1.024 0.875
F_Arachnida-M_Insecta -0.183 -1.110 0.744 0.699 -1.386 1.021
F_Arachnida-M_Mammalia -0.056 -0.619 0.506 0.844 -1.008 0.895
F_Arachnida-M_Reptilia -0.064 -0.780 0.651 0.860 -1.114 0.985
F_Aves-F_Insecta 0.276 -0.098 0.650 0.148 -0.578 1.130
F_Aves-F_Mammalia 0.018 -0.108 0.143 0.782 -0.760 0.795
F_Aves-F_Reptilia 0.107 -0.208 0.422 0.504 -0.722 0.937
F_Aves-M_Actinopterygii 0.104 -0.309 0.518 0.620 -0.767 0.976
F_Aves-M_Aves 0.002 -0.074 0.077 0.967 -0.770 0.773
F_Aves-M_Insecta -0.107 -0.855 0.642 0.780 -1.179 0.965
F_Aves-M_Mammalia 0.020 -0.112 0.152 0.768 -0.759 0.799
F_Aves-M_Reptilia 0.012 -0.449 0.474 0.959 -0.884 0.908
F_Insecta-F_Mammalia -0.258 -0.635 0.118 0.179 -1.113 0.597
F_Insecta-F_Reptilia -0.169 -0.643 0.305 0.485 -1.071 0.733
F_Insecta-M_Actinopterygii -0.171 -0.717 0.374 0.537 -1.113 0.770
F_Insecta-M_Aves -0.274 -0.648 0.100 0.151 -1.128 0.580
F_Insecta-M_Insecta -0.382 -1.165 0.400 0.338 -1.479 0.714
F_Insecta-M_Mammalia -0.256 -0.635 0.123 0.185 -1.112 0.600
F_Insecta-M_Reptilia -0.264 -0.846 0.318 0.374 -1.227 0.700
F_Mammalia-F_Reptilia 0.090 -0.228 0.407 0.580 -0.741 0.920
F_Mammalia-M_Actinopterygii 0.087 -0.329 0.503 0.682 -0.786 0.960
F_Mammalia-M_Aves -0.016 -0.141 0.109 0.802 -0.794 0.762
F_Mammalia-M_Insecta -0.124 -0.874 0.626 0.745 -1.197 0.949
F_Mammalia-M_Mammalia 0.002 -0.099 0.103 0.965 -0.772 0.776
F_Mammalia-M_Reptilia -0.006 -0.469 0.458 0.981 -0.902 0.891
F_Reptilia-M_Actinopterygii -0.003 -0.509 0.503 0.992 -0.922 0.917
F_Reptilia-M_Aves -0.106 -0.420 0.209 0.511 -0.935 0.724
F_Reptilia-M_Insecta -0.214 -1.017 0.590 0.602 -1.325 0.898
F_Reptilia-M_Mammalia -0.087 -0.407 0.233 0.593 -0.919 0.744
F_Reptilia-M_Reptilia -0.095 -0.542 0.352 0.676 -0.983 0.793
M_Actinopterygii-M_Aves -0.103 -0.516 0.311 0.625 -0.975 0.769
M_Actinopterygii-M_Insecta -0.211 -1.058 0.636 0.625 -1.354 0.932
M_Actinopterygii-M_Mammalia -0.085 -0.503 0.334 0.691 -0.959 0.789
M_Actinopterygii-M_Reptilia -0.092 -0.701 0.516 0.766 -1.072 0.887
M_Aves-M_Insecta -0.108 -0.857 0.640 0.777 -1.180 0.964
M_Aves-M_Mammalia 0.018 -0.114 0.151 0.786 -0.761 0.797
M_Aves-M_Reptilia 0.010 -0.451 0.472 0.965 -0.885 0.906
M_Insecta-M_Mammalia 0.126 -0.625 0.878 0.741 -0.948 1.201
M_Insecta-M_Reptilia 0.119 -0.753 0.990 0.789 -1.043 1.280
M_Mammalia-M_Reptilia -0.008 -0.473 0.457 0.974 -0.905 0.890
Code
dat$sex_whose_fitness <- as.factor(paste(dat$sex, dat$whose_fitness ,sep = "_"))

mod_sex_whose_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_whose_fitness - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_whose_fitness)

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

   logLik   Deviance        AIC        BIC       AICc   
-241.1304   482.2607   500.2607   541.0909   500.5254   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0740  0.2721    696     no         effectID 
sigma^2.2  0.0205  0.1433    206     no          paperID 
sigma^2.3  0.0199  0.1410    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 690) = 1.1201, p-val = 0.3487

Model Results:

                                    estimate      se     tval   df    pval 
sex_whose_fitnessBoth_descendant      0.0341  0.1090   0.3130  690  0.7544 
sex_whose_fitnessBoth_focal          -0.0809  0.0433  -1.8666  690  0.0624 
sex_whose_fitnessFemale_descendant   -0.0108  0.0503  -0.2139  690  0.8307 
sex_whose_fitnessFemale_focal        -0.0067  0.0353  -0.1891  690  0.8501 
sex_whose_fitnessMale_descendant      0.0639  0.0568   1.1247  690  0.2611 
sex_whose_fitnessMale_focal          -0.0181  0.0364  -0.4958  690  0.6202 
                                      ci.lb   ci.ub    
sex_whose_fitnessBoth_descendant    -0.1799  0.2481    
sex_whose_fitnessBoth_focal         -0.1660  0.0042  . 
sex_whose_fitnessFemale_descendant  -0.1096  0.0881    
sex_whose_fitnessFemale_focal       -0.0759  0.0625    
sex_whose_fitnessMale_descendant    -0.0477  0.1755    
sex_whose_fitnessMale_focal         -0.0896  0.0535    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_whose_fitness)
   R2_marginal R2_conditional 
   0.009821673    0.359580965 
Code
orchard_plot(mod_sex_whose_fitness, mod = "sex_whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S22. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_whose_fitness <- all_models(mod_sex_whose_fitness,  mod = "sex_whose_fitness", type = "normal")

# saving tab_sex_whose_fitness as RDS
saveRDS(tab_sex_whose_fitness, here("Rdata", "tab_sex_whose_fitness.RDS"))

Table S17. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_whose_fitness <-readRDS(here("Rdata", "tab_sex_whose_fitness.RDS"))

tab_sex_whose_fitness
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_individual -0.078 -0.168 0.011 0.086 -0.849 0.693
B_offspring 0.020 -0.276 0.316 0.894 -0.801 0.841
F_individual -0.011 -0.084 0.062 0.760 -0.781 0.758
F_offspring -0.021 -0.129 0.087 0.704 -0.794 0.753
M_individual -0.034 -0.109 0.041 0.376 -0.803 0.736
M_offspring 0.061 -0.063 0.185 0.333 -0.715 0.837
B_individual-B_offspring 0.098 -0.198 0.395 0.515 -0.723 0.920
B_individual-F_individual 0.067 -0.028 0.162 0.168 -0.705 0.839
B_individual-F_offspring 0.057 -0.067 0.182 0.365 -0.719 0.833
B_individual-M_individual 0.045 -0.053 0.142 0.370 -0.728 0.817
B_individual-M_offspring 0.140 0.002 0.277 0.047 -0.639 0.918
B_offspring-F_individual -0.031 -0.329 0.266 0.836 -0.853 0.790
B_offspring-F_offspring -0.041 -0.349 0.268 0.795 -0.867 0.785
B_offspring-M_individual -0.054 -0.352 0.244 0.723 -0.876 0.768
B_offspring-M_offspring 0.041 -0.273 0.355 0.797 -0.787 0.869
F_individual-F_offspring -0.010 -0.112 0.093 0.855 -0.782 0.763
F_individual-M_individual -0.022 -0.087 0.043 0.499 -0.791 0.746
F_individual-M_offspring 0.073 -0.046 0.192 0.232 -0.703 0.848
F_offspring-M_individual -0.013 -0.118 0.092 0.810 -0.786 0.760
F_offspring-M_offspring 0.082 -0.051 0.215 0.227 -0.695 0.860
M_individual-M_offspring 0.095 -0.025 0.215 0.121 -0.680 0.870
Code
dat$sex_fitness_higher_level <- as.factor(paste(dat$sex, dat$fitness_higher_level ,sep = "_"))

mod_sex_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_fitness_higher_level - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_fitness)

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

   logLik   Deviance        AIC        BIC       AICc   
-240.7301   481.4602   503.4602   553.3319   503.8507   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0738  0.2717    696     no         effectID 
sigma^2.2  0.0246  0.1569    206     no          paperID 
sigma^2.3  0.0149  0.1220    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 688) = 1.0605, p-val = 0.3891

Model Results:

                                                 estimate      se     tval   df 
sex_fitness_higher_levelBoth_reproduction         -0.0212  0.0631  -0.3357  688 
sex_fitness_higher_levelBoth_survival             -0.0917  0.0480  -1.9089  688 
sex_fitness_higher_levelFemale_indirect fitness   -0.2306  0.2135  -1.0801  688 
sex_fitness_higher_levelFemale_reproduction        0.0071  0.0368   0.1927  688 
sex_fitness_higher_levelFemale_survival           -0.0291  0.0402  -0.7230  688 
sex_fitness_higher_levelMale_indirect fitness     -0.2278  0.1748  -1.3031  688 
sex_fitness_higher_levelMale_reproduction         -0.0171  0.0399  -0.4287  688 
sex_fitness_higher_levelMale_survival              0.0163  0.0418   0.3886  688 
                                                   pval    ci.lb   ci.ub    
sex_fitness_higher_levelBoth_reproduction        0.7372  -0.1452  0.1028    
sex_fitness_higher_levelBoth_survival            0.0567  -0.1860  0.0026  . 
sex_fitness_higher_levelFemale_indirect fitness  0.2805  -0.6498  0.1886    
sex_fitness_higher_levelFemale_reproduction      0.8473  -0.0651  0.0793    
sex_fitness_higher_levelFemale_survival          0.4700  -0.1080  0.0499    
sex_fitness_higher_levelMale_indirect fitness    0.1930  -0.5710  0.1154    
sex_fitness_higher_levelMale_reproduction        0.6682  -0.0954  0.0612    
sex_fitness_higher_levelMale_survival            0.6977  -0.0659  0.0984    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_fitness)
   R2_marginal R2_conditional 
    0.01181556     0.35629363 
Code
orchard_plot(mod_sex_fitness, mod = "sex_fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S23. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_fitness <- all_models(mod_sex_fitness,  mod = "sex_fitness_higher_level", type = "normal")

# saving tab_sex_fitness as RDS
saveRDS(tab_sex_fitness, here("Rdata", "tab_sex_fitness.RDS"))

Table S18. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_fitness <-readRDS(here("Rdata", "tab_sex_fitness.RDS"))

tab_sex_fitness
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_reproduction -0.023 -0.164 0.118 0.749 -0.800 0.754
B_survival -0.095 -0.196 0.006 0.066 -0.865 0.676
F_indirect fitness -0.231 -0.732 0.270 0.365 -1.144 0.682
F_reproduction -0.004 -0.081 0.073 0.922 -0.771 0.764
F_survival -0.026 -0.112 0.061 0.558 -0.794 0.743
M_indirect fitness -0.230 -0.638 0.178 0.269 -1.096 0.636
M_reproduction -0.028 -0.112 0.056 0.513 -0.796 0.740
M_survival -0.001 -0.089 0.088 0.988 -0.769 0.768
B_reproduction-B_survival -0.072 -0.231 0.088 0.378 -0.852 0.708
B_reproduction-F_indirect fitness -0.208 -0.726 0.310 0.431 -1.131 0.715
B_reproduction-F_reproduction 0.019 -0.132 0.170 0.803 -0.759 0.798
B_reproduction-F_survival -0.003 -0.159 0.153 0.972 -0.782 0.777
B_reproduction-M_indirect fitness -0.207 -0.636 0.223 0.345 -1.083 0.669
B_reproduction-M_reproduction -0.005 -0.159 0.149 0.949 -0.784 0.774
B_reproduction-M_survival 0.022 -0.135 0.180 0.781 -0.757 0.802
B_survival-F_indirect fitness -0.136 -0.644 0.371 0.598 -1.053 0.781
B_survival-F_reproduction 0.091 -0.017 0.199 0.098 -0.680 0.862
B_survival-F_survival 0.069 -0.045 0.183 0.233 -0.703 0.841
B_survival-M_indirect fitness -0.135 -0.552 0.282 0.525 -1.005 0.735
B_survival-M_reproduction 0.067 -0.046 0.180 0.247 -0.705 0.839
B_survival-M_survival 0.094 -0.023 0.211 0.114 -0.678 0.867
F_indirect fitness-F_reproduction 0.227 -0.275 0.729 0.375 -0.687 1.141
F_indirect fitness-F_survival 0.205 -0.297 0.708 0.423 -0.709 1.119
F_indirect fitness-M_indirect fitness 0.001 -0.602 0.604 0.997 -0.972 0.974
F_indirect fitness-M_reproduction 0.203 -0.300 0.706 0.428 -0.711 1.118
F_indirect fitness-M_survival 0.230 -0.272 0.733 0.369 -0.684 1.145
F_reproduction-F_survival -0.022 -0.104 0.060 0.600 -0.790 0.746
F_reproduction-M_indirect fitness -0.226 -0.636 0.184 0.279 -1.093 0.641
F_reproduction-M_reproduction -0.024 -0.102 0.054 0.542 -0.792 0.743
F_reproduction-M_survival 0.003 -0.085 0.091 0.944 -0.766 0.772
F_survival-M_indirect fitness -0.204 -0.615 0.207 0.330 -1.071 0.663
F_survival-M_reproduction -0.002 -0.093 0.089 0.962 -0.771 0.767
F_survival-M_survival 0.025 -0.064 0.114 0.579 -0.744 0.794
M_indirect fitness-M_reproduction 0.202 -0.208 0.612 0.334 -0.665 1.069
M_indirect fitness-M_survival 0.229 -0.182 0.640 0.274 -0.638 1.096
M_reproduction-M_survival 0.027 -0.067 0.122 0.570 -0.742 0.797
Code
dat$sex_dispersal_type <- as.factor(paste(dat$sex, dat$dispersal_type ,sep = "_"))

mod_sex_dispersal_type <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_dispersal_type - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_dispersal_type)

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

   logLik   Deviance        AIC        BIC       AICc   
-241.4964   482.9928   506.9928   561.3808   507.4557   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0745  0.2729    696     no         effectID 
sigma^2.2  0.0304  0.1744    206     no          paperID 
sigma^2.3  0.0087  0.0935    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:9):
F(df1 = 9, df2 = 687) = 0.5900, p-val = 0.8059

Model Results:

                                   estimate      se     tval   df    pval 
sex_dispersal_typeBoth_Both         -0.1189  0.1211  -0.9817  687  0.3266 
sex_dispersal_typeBoth_breeding      0.0131  0.0846   0.1544  687  0.8773 
sex_dispersal_typeBoth_natal        -0.0838  0.0475  -1.7628  687  0.0784 
sex_dispersal_typeFemale_Both       -0.0864  0.1219  -0.7088  687  0.4787 
sex_dispersal_typeFemale_breeding   -0.0071  0.0498  -0.1431  687  0.8863 
sex_dispersal_typeFemale_natal      -0.0019  0.0366  -0.0518  687  0.9587 
sex_dispersal_typeMale_Both         -0.0056  0.1386  -0.0401  687  0.9680 
sex_dispersal_typeMale_breeding     -0.0257  0.0533  -0.4817  687  0.6302 
sex_dispersal_typeMale_natal         0.0083  0.0388   0.2145  687  0.8302 
                                     ci.lb   ci.ub    
sex_dispersal_typeBoth_Both        -0.3566  0.1189    
sex_dispersal_typeBoth_breeding    -0.1530  0.1791    
sex_dispersal_typeBoth_natal       -0.1772  0.0095  . 
sex_dispersal_typeFemale_Both      -0.3257  0.1529    
sex_dispersal_typeFemale_breeding  -0.1049  0.0907    
sex_dispersal_typeFemale_natal     -0.0738  0.0700    
sex_dispersal_typeMale_Both        -0.2777  0.2666    
sex_dispersal_typeMale_breeding    -0.1303  0.0790    
sex_dispersal_typeMale_natal       -0.0679  0.0845    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_dispersal_type)
   R2_marginal R2_conditional 
   0.009425299    0.350676406 
Code
orchard_plot(mod_sex_dispersal_type, mod = "sex_dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, trunk.size = 0.5, angle = 45)

Figure S24. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types combined). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_dispersal_type <- all_models(mod_sex_dispersal_type,  mod = "sex_dispersal_type", type = "normal")

# saving tab_sex_dispersal_type as RDS
saveRDS(tab_sex_dispersal_type, here("Rdata", "tab_sex_dispersal_type.RDS"))

Table S19. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types grouped) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_dispersal_type <-readRDS(here("Rdata", "tab_sex_dispersal_type.RDS"))

tab_sex_dispersal_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_B -0.096 -0.321 0.129 0.403 -0.892 0.699
B_breeding -0.004 -0.200 0.192 0.967 -0.792 0.784
B_natal -0.083 -0.184 0.017 0.103 -0.853 0.686
F_B -0.081 -0.301 0.138 0.468 -0.875 0.713
F_breeding -0.004 -0.110 0.102 0.940 -0.774 0.766
F_natal -0.008 -0.086 0.070 0.840 -0.775 0.759
M_B -0.108 -0.350 0.133 0.379 -0.909 0.692
M_breeding -0.022 -0.137 0.094 0.712 -0.793 0.750
M_natal -0.004 -0.086 0.077 0.922 -0.771 0.763
B_B-B_breeding 0.092 -0.203 0.387 0.540 -0.726 0.910
B_B-B_natal 0.013 -0.223 0.249 0.916 -0.786 0.811
B_B-F_B 0.015 -0.288 0.317 0.923 -0.806 0.836
B_B-F_breeding 0.092 -0.150 0.334 0.456 -0.708 0.892
B_B-F_natal 0.088 -0.142 0.318 0.453 -0.709 0.885
B_B-M_B -0.012 -0.330 0.305 0.939 -0.839 0.814
B_B-M_breeding 0.074 -0.172 0.321 0.554 -0.727 0.876
B_B-M_natal 0.092 -0.139 0.323 0.435 -0.705 0.889
B_breeding-B_natal -0.079 -0.292 0.133 0.464 -0.871 0.713
B_breeding-F_B -0.077 -0.368 0.214 0.603 -0.893 0.739
B_breeding-F_breeding 0.000 -0.214 0.214 1.000 -0.792 0.792
B_breeding-F_natal -0.004 -0.208 0.200 0.970 -0.794 0.786
B_breeding-M_B -0.104 -0.413 0.204 0.506 -0.927 0.718
B_breeding-M_breeding -0.018 -0.238 0.202 0.875 -0.812 0.776
B_breeding-M_natal 0.000 -0.206 0.206 1.000 -0.790 0.790
B_natal-F_B 0.002 -0.232 0.237 0.985 -0.796 0.800
B_natal-F_breeding 0.079 -0.051 0.209 0.232 -0.695 0.853
B_natal-F_natal 0.075 -0.032 0.183 0.169 -0.695 0.846
B_natal-M_B -0.025 -0.281 0.231 0.848 -0.830 0.780
B_natal-M_breeding 0.062 -0.078 0.202 0.386 -0.714 0.837
B_natal-M_natal 0.079 -0.031 0.190 0.158 -0.692 0.850
F_B-F_breeding 0.077 -0.160 0.314 0.523 -0.722 0.876
F_B-F_natal 0.073 -0.153 0.300 0.526 -0.723 0.869
F_B-M_B -0.027 -0.288 0.234 0.837 -0.834 0.779
F_B-M_breeding 0.059 -0.183 0.302 0.630 -0.741 0.860
F_B-M_natal 0.077 -0.150 0.304 0.506 -0.719 0.873
F_breeding-F_natal -0.004 -0.113 0.105 0.944 -0.775 0.767
F_breeding-M_B -0.104 -0.362 0.154 0.427 -0.910 0.701
F_breeding-M_breeding -0.018 -0.130 0.095 0.758 -0.789 0.754
F_breeding-M_natal 0.000 -0.112 0.112 1.000 -0.771 0.771
F_natal-M_B -0.100 -0.349 0.148 0.427 -0.903 0.702
F_natal-M_breeding -0.014 -0.134 0.107 0.823 -0.786 0.759
F_natal-M_natal 0.004 -0.068 0.076 0.916 -0.762 0.770
M_B-M_breeding 0.087 -0.176 0.350 0.517 -0.720 0.894
M_B-M_natal 0.104 -0.145 0.354 0.412 -0.698 0.907
M_breeding-M_natal 0.018 -0.106 0.142 0.780 -0.755 0.791

5 Species-level models

There are two steps, which can be done using phylo_model().

5.1 Estimate species-level effect sizes.

Code
# function to estimate species-level effect sizes
phylo_model <- function(data = dat, tree, correlation = 0.5) {
  
  # accounting for within-species correlation
  Vm <- vcalc(vi = data$vi, cluster = data$shared_group, obs = data$effectID, rho = correlation)
  mod <- rma.mv(yi = yi, 
              V = Vm,
              mod = ~ 1,
              data = data,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | phylogeny),
              #R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)
  
  
  # get random effects
  # dat3 <- dat2[c("phylogeny", "yi", "vi")]
  REs <- ranef(mod)$phylogeny
  
  # calculate means and CIs
  dat <- data.frame(phylogeny = rownames(REs),
                    estimate = mod$beta[1] + REs$intrcpt,
                    se = sqrt((mod$se[1])^2 + (REs$se)^2)) %>%
    mutate(lb = estimate - se * qnorm(1 - 0.05/2), 
           ub = estimate + se * qnorm(1 - 0.05/2)) %>%
    arrange(phylogeny)
  # get species class data
  species_class <- dplyr::distinct(data, phylogeny, .keep_all = TRUE) %>%
    dplyr::select(species_class, phylogeny)
  
  # extract tip labels from the tree and merge with species class
  tip.label <- data.frame(phylogeny = tree$tip.label)  # extract tip label
  species_class2 <- left_join(tip.label, species_class, by = "phylogeny")
  
  return(list(data = dat, class = species_class2))  # return processed data and species class
}

# estimate species-level effect sizes
# you can specify your own correlation coefficient via `correlation`. For illustration purpose, I used 0.5.
result <- phylo_model(data = dat, tree = tree, correlation = 0.5)

# find the significant species
result$data <- result$data %>%
  mutate(
    p_value = 2 * (1 - pnorm(abs(estimate / se))),
    significance = ifelse(p_value < 0.05, "Significant", "NonSignificant")
  )
result$data <- result$data %>%
  mutate(significance_label = ifelse(p_value < 0.05, "*", NA))

5.2 Plot tree and add species-level effect size estimates

Code
# plot tree
p1 <- ggtree(tree, layout = "rectangular", cex = 0.3)

p2 <- p1 %<+% result$class + geom_tiplab(aes(color = species_class), size = 3, align = T, offset = 0.05) + geom_tippoint(aes(color = species_class)) + guides(color = "none")

p3 <- p2 + 
  geom_facet(
    panel = "Effect size", 
    data = result$data, 
    geom = ggstance::geom_pointrangeh,
    mapping = aes(x = estimate, xmin = lb, xmax = ub, color = species_class)
  ) +
  geom_facet(
    panel = "Effect size",
    data = result$data,
    geom = geom_text,
    mapping = aes(x = lb, label = significance_label), 
    hjust = 1.5, size = 3, color = "red"              # adjust position of star *
  ) +
  geom_facet(
    panel = "Effect size",
    data = result$data, 
    geom = geom_vline,
    mapping = aes(xintercept = 0),                    
    linetype = "dashed", color = "gray"
  ) +
  theme_tree2()
Code
# adjust widths of the plot facets as necessary to improve visualization 

facet_widths(p3, c(Tree = 0.6, `Effect size` = 0.4))

Figure S25. Phylogenetic tree (left panel) and species-level effect size estimates (right panel). Significant relationships are indicated with a red asterisk. Dispersal was positive for ruffed grouse (Bonasa umbellus; Zr = 0.42, 95% CIs: 0.16 to 0.68) and negative for brown jays (Cyanocorax morio; Zr = -0.30, 95% CIs: -0.58 to -0.01), Tasmanian nativehens (Gallinula mortierii; Zr = -0.34, 95% CIs: -0.60 to X), and cougars (Puma concolor; Zr = -0.61, 95% CIs: -0.90 to -0.32).

6 Multi-moderator-regression model

Code
# used "maximum likelihood" for model selection

#######################
# Multi-variable models
#######################

mod_full <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1 + sex +
                age_class_clean + whose_fitness + 
                fitness_higher_level +
                dispersal_type + dispersal_phase,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              method = "ML",
              test = "t",
              sparse = TRUE)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-238.0923  2849.3379   504.1845   567.8194   504.8013   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0728  0.2699    696     no         effectID 
sigma^2.2  0.0173  0.1316    206     no          paperID 
sigma^2.3  0.0225  0.1501    148     no  species_cleaned 

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

Test of Moderators (coefficients 2:11):
F(df1 = 10, df2 = 685) = 1.1945, p-val = 0.2910

Model Results:

                                  estimate      se     tval   df    pval 
intrcpt                            -0.2420  0.1744  -1.3872  685  0.1658 
sexFemale                           0.0461  0.0445   1.0352  685  0.3009 
sexMale                             0.0531  0.0452   1.1761  685  0.2400 
age_class_cleanjuvenile             0.0790  0.0558   1.4171  685  0.1569 
age_class_cleanmix                 -0.0691  0.0622  -1.1124  685  0.2664 
whose_fitnessfocal                 -0.0565  0.0415  -1.3603  685  0.1742 
fitness_higher_levelreproduction    0.2257  0.1445   1.5621  685  0.1187 
fitness_higher_levelsurvival        0.1878  0.1465   1.2814  685  0.2005 
dispersal_typebreeding             -0.0120  0.0894  -0.1348  685  0.8928 
dispersal_typenatal                -0.0084  0.0817  -0.1022  685  0.9186 
dispersal_phasesettlement           0.0374  0.0501   0.7461  685  0.4559 
                                    ci.lb   ci.ub    
intrcpt                           -0.5844  0.1005    
sexFemale                         -0.0413  0.1335    
sexMale                           -0.0356  0.1418    
age_class_cleanjuvenile           -0.0305  0.1886    
age_class_cleanmix                -0.1912  0.0529    
whose_fitnessfocal                -0.1381  0.0251    
fitness_higher_levelreproduction  -0.0580  0.5095    
fitness_higher_levelsurvival      -0.0999  0.4755    
dispersal_typebreeding            -0.1875  0.1634    
dispersal_typenatal               -0.1688  0.1521    
dispersal_phasesettlement         -0.0610  0.1358    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_full)*100, 2)
   R2_marginal R2_conditional 
          2.17          36.75 
Code
# multi-model selection
candidates <- dredge(mod_full, trace = 2)

# saving candidates as RDS
saveRDS(candidates, here("Rdata", "candidates.RDS"))
Code
candidates <-readRDS(here("Rdata", "candidates.RDS"))

candidates
Global model call: rma.mv(yi = yi, V = VCV, mods = ~1 + sex + age_class_clean + 
    whose_fitness + fitness_higher_level + dispersal_type + dispersal_phase, 
    random = list(~1 | effectID, ~1 | paperID, ~1 | species_cleaned), 
    data = dat, method = "ML", test = "t", sparse = TRUE)
---
Model selection table 
   (Int) age_cls_cln dsp_phs dsp_typ ftn_hgh_lvl sex whs_ftn df   logLik  AICc
2      +           +                                          6 -276.828 565.8
4      +           +       +                                  7 -275.933 566.0
34     +           +                                       +  7 -276.278 566.7
36     +           +       +                               +  8 -275.479 567.2
1      +                                                      4 -279.768 567.6
42     +           +                           +           +  9 -274.701 567.7
10     +           +                           +              8 -275.782 567.8
18     +           +                               +          8 -275.927 568.1
12     +           +       +                   +              9 -275.053 568.4
33     +                                                   +  5 -279.217 568.5
3      +                   +                                  5 -279.276 568.6
44     +           +       +                   +           + 10 -274.319 569.0
20     +           +       +                       +          9 -275.359 569.0
50     +           +                               +       +  9 -275.518 569.3
6      +           +               +                          8 -276.641 569.5
9      +                                       +              6 -278.813 569.8
35     +                   +                               +  6 -278.815 569.8
41     +                                       +           +  7 -277.874 569.9
8      +           +       +       +                          9 -275.830 569.9
17     +                                           +          6 -278.913 570.0
26     +           +                           +   +         10 -274.949 570.2
52     +           +       +                       +       + 10 -275.002 570.3
38     +           +               +                       +  9 -276.064 570.4
5      +                           +                          6 -279.306 570.7
58     +           +                           +   +       + 11 -274.187 570.8
40     +           +       +       +                       + 10 -275.349 571.0
11     +                   +                   +              7 -278.449 571.1
49     +                                           +       +  7 -278.515 571.2
28     +           +       +                   +   +         11 -274.474 571.3
46     +           +               +           +           + 11 -274.557 571.5
19     +                   +                       +          7 -278.735 571.6
14     +           +               +           +             10 -275.674 571.7
37     +                           +                       +  7 -278.774 571.7
43     +                   +                   +           +  8 -277.769 571.8
22     +           +               +               +         10 -275.770 571.9
7      +                   +       +                          7 -278.938 572.0
25     +                                       +   +          8 -278.038 572.3
60     +           +       +                   +   +       + 12 -273.929 572.3
16     +           +       +       +           +             11 -274.996 572.4
48     +           +       +       +           +           + 12 -274.220 572.9
24     +           +       +       +               +         11 -275.261 572.9
51     +                   +                       +       +  8 -278.367 572.9
54     +           +               +               +       + 11 -275.338 573.1
57     +                                       +   +       +  9 -277.439 573.1
13     +                           +           +              8 -278.487 573.2
39     +                   +       +                       +  8 -278.490 573.2
21     +                           +               +          8 -278.571 573.4
45     +                           +           +           +  9 -277.602 573.5
27     +                   +                   +   +          9 -277.898 574.1
30     +           +               +           +   +         12 -274.856 574.2
56     +           +       +       +               +       + 12 -274.881 574.2
53     +                           +               +       +  9 -278.186 574.6
15     +                   +       +           +              9 -278.197 574.7
62     +           +               +           +   +       + 13 -274.062 574.7
59     +                   +                   +   +       + 10 -277.406 575.1
23     +                   +       +               +          9 -278.444 575.2
32     +           +       +       +           +   +         13 -274.418 575.4
47     +                   +       +           +           + 10 -277.534 575.4
29     +                           +           +   +         10 -277.795 575.9
64     +           +       +       +           +   +       + 14 -273.837 576.3
55     +                   +       +               +       + 10 -278.088 576.5
61     +                           +           +   +       + 11 -277.225 576.8
31     +                   +       +           +   +         11 -277.687 577.8
63     +                   +       +           +   +       + 12 -277.207 578.9
   delta weight
2   0.00  0.122
4   0.25  0.107
34  0.94  0.076
36  1.39  0.061
1   1.81  0.049
42  1.89  0.047
10  2.00  0.045
18  2.29  0.039
12  2.59  0.033
33  2.74  0.031
3   2.86  0.029
44  3.19  0.025
20  3.21  0.024
50  3.52  0.021
6   3.72  0.019
9   3.97  0.017
35  3.97  0.017
41  4.13  0.015
8   4.15  0.015
17  4.17  0.015
26  4.45  0.013
52  4.55  0.012
38  4.62  0.012
5   4.96  0.010
58  4.99  0.010
40  5.25  0.009
11  5.28  0.009
49  5.42  0.008
28  5.56  0.008
46  5.73  0.007
19  5.86  0.007
14  5.90  0.006
37  5.93  0.006
43  5.97  0.006
22  6.09  0.006
7   6.26  0.005
25  6.51  0.005
60  6.55  0.005
16  6.61  0.004
48  7.13  0.003
24  7.14  0.003
51  7.17  0.003
54  7.29  0.003
57  7.37  0.003
13  7.41  0.003
39  7.42  0.003
21  7.58  0.003
45  7.69  0.003
27  8.28  0.002
30  8.40  0.002
56  8.45  0.002
53  8.86  0.001
15  8.88  0.001
62  8.89  0.001
59  9.36  0.001
23  9.38  0.001
32  9.60  0.001
47  9.62  0.001
29 10.14  0.001
64 10.53  0.001
55 10.72  0.001
61 11.07  0.000
31 11.99  0.000
63 13.10  0.000
Models ranked by AICc(x) 
Code
# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 2) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 2)) 

mr_averaged_aic2

Call:
model.avg(object = candidates, subset = delta < 2)

Component model call: 
rma.mv(yi = yi, V = VCV, mods = ~<7 unique rhs>, random = list(~1 | 
     effectID, ~1 | paperID, ~1 | species_cleaned), data = dat, method = ML, 
     test = t, sparse = TRUE)

Component models: 
       df  logLik   AICc delta weight
1       6 -276.83 565.78  0.00   0.24
12      7 -275.93 566.03  0.25   0.21
14      7 -276.28 566.72  0.94   0.15
124     8 -275.48 567.17  1.39   0.12
(Null)  4 -279.77 567.59  1.81   0.10
134     9 -274.70 567.67  1.89   0.09
13      8 -275.78 567.78  2.00   0.09

Term codes: 
     age_class_clean      dispersal_phase fitness_higher_level 
                   1                    2                    3 
       whose_fitness 
                   4 

Model-averaged coefficients:  
(full average) 
                                 Estimate Std. Error z value Pr(>|z|)
intrcpt                          -0.08622    0.11029   0.782    0.434
age_class_cleanjuvenile           0.06241    0.05321   1.173    0.241
age_class_cleanmix               -0.10467    0.07067   1.481    0.139
dispersal_phasesettlement         0.02130    0.04137   0.515    0.607
whose_fitnessoffspring            0.01788    0.03542   0.505    0.614
fitness_higher_levelreproduction  0.03936    0.11059   0.356    0.722
fitness_higher_levelsurvival      0.03285    0.10117   0.325    0.745
 
(conditional average) 
                                 Estimate Std. Error z value Pr(>|z|)  
intrcpt                          -0.08622    0.11029   0.782   0.4344  
age_class_cleanjuvenile           0.06911    0.05169   1.337   0.1812  
age_class_cleanmix               -0.11591    0.06502   1.783   0.0746 .
dispersal_phasesettlement         0.06429    0.04902   1.312   0.1896  
whose_fitnessoffspring            0.04929    0.04369   1.128   0.2592  
fitness_higher_levelreproduction  0.21658    0.17002   1.274   0.2027  
fitness_higher_levelsurvival      0.18075    0.17198   1.051   0.2933  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# relative importance of each predictor for all the models
importance <- sw(candidates)

importance
                     age_class_clean dispersal_phase whose_fitness
Sum of weights:      0.74            0.40            0.40         
N containing models:   32              32              32         
                     fitness_higher_level sex  dispersal_type
Sum of weights:      0.28                 0.20 0.14          
N containing models:   32                   32   32          

7 Publication bias

Code
# if each group n is available - assume n/2
dat$effectN <- ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n,  
                      (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))
  
dat$sqeffectN <- sqrt(dat$effectN)

mod_egger <- rma.mv(yi = yi, 
                    V = VCV,
                mod = ~ sqeffectN,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_egger)

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

   logLik   Deviance        AIC        BIC       AICc   
-243.1309   486.2619   496.2619   518.9742   496.3491   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0747  0.2734    696     no         effectID 
sigma^2.2  0.0225  0.1500    206     no          paperID 
sigma^2.3  0.0165  0.1283    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.0281, p-val = 0.8668

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -0.0282  0.0407  -0.6920  694  0.4891  -0.1082  0.0518    
sqeffectN    0.0007  0.0041   0.1677  694  0.8668  -0.0074  0.0088    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_egger)*100, 2)
   R2_marginal R2_conditional 
          0.01          34.28 
Code
small <- bubble_plot(mod_egger,
                     mod = "sqeffectN",
                     group = "paperID",
                     xlab = "sqrt(Effective N)",
                     ylab = "Effect Size: Zr",
                     g = TRUE)
Code
small

Figure S26. Relationship between sample size and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# creating publication year from "reference"
dat$reference <- as.character(dat$reference)
dat$year <- with(dat, substr(reference, nchar(reference)-4, nchar(reference)))
dat$year <- as.integer(ifelse(dat$year == "2017a" | dat$year == "2017b", 2017, dat$year))
# decline effect
mod_dec <- rma.mv(yi = yi, V = vi,
                     mod = ~ year,
                     data = dat, 
                     random = list(
                       ~ 1 | effectID,
                       ~ 1 | paperID,
                       ~ 1 | species_cleaned),
                     test = "t",
                     sparse = TRUE)
summary(mod_dec)

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

   logLik   Deviance        AIC        BIC       AICc   
-223.4412   446.8824   456.8824   479.5947   456.9696   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0615  0.2480    696     no         effectID 
sigma^2.2  0.0158  0.1258    206     no          paperID 
sigma^2.3  0.0222  0.1490    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.4765, p-val = 0.4903

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt   -2.1353  3.0623  -0.6973  694  0.4859  -8.1478  3.8773    
year       0.0011  0.0015   0.6903  694  0.4903  -0.0019  0.0040    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_dec)*100, 2)
   R2_marginal R2_conditional 
          0.17          38.32 
Code
decline <- bubble_plot(mod_dec,
                       mod = "year",
                       group = "paperID",
                       xlab = "Publication year",
                       ylab = "Effect Size: Zr",
                       g = TRUE)
Code
decline 

Figure S27. Relationship between publication year and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# All together

mod_comb <- rma.mv(yi = yi, 
                   V = VCV,
                   mod = ~ year + sqeffectN,
                   data = dat, 
                   random = list(
                     ~ 1 | effectID,
                     ~ 1 | paperID,
                     ~ 1 | species_cleaned),
                   # ~ 1 | phylogeny),
                   # R= list(phylogeny = cor_tree),
                   test = "t",
                   sparse = TRUE)
                   
summary(mod_comb)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.6991   485.3982   497.3982   524.6444   497.5206   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0744  0.2728    696     no         effectID 
sigma^2.2  0.0199  0.1411    206     no          paperID 
sigma^2.3  0.0210  0.1447    148     no  species_cleaned 

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

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 693) = 0.1902, p-val = 0.8268

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -1.9866  3.2641  -0.6086  693  0.5430  -8.3953  4.4220    
year         0.0010  0.0016   0.5997  693  0.5489  -0.0022  0.0042    
sqeffectN    0.0003  0.0042   0.0608  693  0.9515  -0.0079  0.0084    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_comb)*100, 2)
   R2_marginal R2_conditional 
          0.13          35.52 
Code
# combined sample size and year
bubble_plot(mod_comb,
            mod = "sqeffectN",
            group = "paperID",
            xlab = "sqrt(Effective N)",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S28. Combined effect of sample size and publication year on effect size (plotting sample size only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# combined sample size and year
bubble_plot(mod_comb,
            mod = "year",
            group = "paperID",
            xlab = "Publication year",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S29. Combined effect of sample size and publication year on effect size (plotting publication year only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
dat <- dat %>%
  mutate(leave_out = reference)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(model) {
  df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
  return(df)
}

# using dplyr to form data frame
MA_oneout <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout$left_out <- as.factor(MA_oneout$left_out)
MA_oneout$left_out <- factor(MA_oneout$left_out, levels = MA_oneout$left_out)

# saving the runs
saveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))
Code
MA_oneout <- readRDS(here("Rdata", "MA_oneout.RDS"))

# plotting
leaveoneout <- ggplot(MA_oneout) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Study left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
Code
leaveoneout

Figure S30. Leave one study out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals.
Code
dat <- dat %>%
  mutate(leave_out = species_cleaned)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# # writing function for extracting est, ci.lb, and ci.ub from all models
# est.func <- function(model) {
#   df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
#   return(df)
# }

# using dplyr to form data frame
MA_oneout1 <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout1$left_out <- as.factor(MA_oneout1$left_out)
MA_oneout1$left_out <- factor(MA_oneout1$left_out, levels = MA_oneout1$left_out)

# saving the runs
saveRDS(MA_oneout1, here("Rdata", "MA_oneout1.RDS"))
Code
MA_oneout1 <- readRDS(here("Rdata", "MA_oneout1.RDS"))

# plotting
leaveoneout1 <- ggplot(MA_oneout1) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Species left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
Code
leaveoneout1

Figure S31. Leave one species out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals.
Code
precision_plot<-funnel(mod1, 
       yaxis="seinv",
       # = "rstudent",
       xlab = "Standarized residuals",
       ylab = "Precision (inverse of SE)",
       ylim = c(0.001, 16),
       xlim = c(-10,15),
       digits=c(0,1))

Figure S32. Relationship between standardized residuals and precision (inverse of SE).

8 R Session Information

Code
# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.4.3 (2025-02-28)

Platform: x86_64-apple-darwin20

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: cowplot(v.1.1.3), matrixcalc(v.1.0-6), orchaRd(v.2.0), here(v.1.0.1), patchwork(v.1.3.0), ggstance(v.0.3.7), ggtree(v.3.14.0), ape(v.5.8-1), metafor(v.4.8-0), numDeriv(v.2016.8-1.1), metadat(v.1.4-0), Matrix(v.1.7-2), pander(v.0.6.5), gridExtra(v.2.3), emmeans(v.1.10.7), MuMIn(v.1.48.4), kableExtra(v.1.4.0), lubridate(v.1.9.4), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.4), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): tidyselect(v.1.2.1), viridisLite(v.0.4.2), vipor(v.0.4.7), farver(v.2.1.2), latex2exp(v.0.9.6), fastmap(v.1.2.0), lazyeval(v.0.2.2), TH.data(v.1.1-3), pacman(v.0.5.1), mathjaxr(v.1.6-0), digest(v.0.6.37), estimability(v.1.5.1), timechange(v.0.3.0), lifecycle(v.1.0.4), survival(v.3.8-3), tidytree(v.0.4.6), magrittr(v.2.0.3), compiler(v.4.4.3), rlang(v.1.1.5), tools(v.4.4.3), yaml(v.2.3.10), knitr(v.1.49), labeling(v.0.4.3), htmlwidgets(v.1.6.4), xml2(v.1.3.7), aplot(v.0.2.4), multcomp(v.1.4-28), withr(v.3.0.2), grid(v.4.4.3), stats4(v.4.4.3), xtable(v.1.8-4), colorspace(v.2.1-1), scales(v.1.3.0), MASS(v.7.3-64), cli(v.3.6.4), mvtnorm(v.1.3-3), rmarkdown(v.2.29), treeio(v.1.30.0), generics(v.0.1.3), rstudioapi(v.0.17.1), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), splines(v.4.4.3), parallel(v.4.4.3), ggplotify(v.0.1.2), vctrs(v.0.6.5), yulab.utils(v.0.1.9), sandwich(v.3.1-1), jsonlite(v.1.9.0), gridGraphics(v.0.5-1), hms(v.1.1.3), beeswarm(v.0.4.0), systemfonts(v.1.2.1), glue(v.1.8.0), codetools(v.0.2-20), stringi(v.1.8.4), gtable(v.0.3.6), munsell(v.0.5.1), pillar(v.1.10.1), htmltools(v.0.5.8.1), R6(v.2.6.1), rprojroot(v.2.0.4), evaluate(v.1.0.3), lattice(v.0.22-6), ggfun(v.0.1.8), Rcpp(v.1.0.14), svglite(v.2.1.3), coda(v.0.19-4.1), nlme(v.3.1-167), mgcv(v.1.9-1), xfun(v.0.51), fs(v.1.6.5), zoo(v.1.8-13) and pkgconfig(v.2.0.3)