Does philopatry pay? Meta-analysis reveals little empirical support for dispersal being costly.

Supplmentary Material

Things to do for April

1. Rerun all the models
2. Put figures at the end
3. Add some text so the reader can understand
4. Go through code and add more annotations
5. Ask questions about things you do not understand

1 Setting up

1.1 Load 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
               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
               # more too add
)

eval(metafor:::.MuMIn)

1.2 Load data: paper and tree data

#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

# 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, 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, 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, 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

With all five random effects

# 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"))

# 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

# 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 

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.

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

Reduced multilevel model

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

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-280.0131   560.0262   568.0262   586.1263   568.0853   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3316    683     no         effectID 
sigma^2.2  0.0113  0.1061    204     no          paperID 
sigma^2.3  0.0274  0.1656    148     no  species_cleaned 

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

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0280  0.0300  -0.9335  682  0.3509  -0.0869  0.0309    

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

round(i2_ml(mod1),2)   

          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.69              72.27               7.40              18.02 

reduced_intercept<-orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

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.

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

tab_mod1

Fixed effect	estimate	lowerCL	upperCL	pval	lowerPR	upperPR
intercept	-0.028	-0.087	0.031	0.351	-0.787	0.731

3 Uni-moderator meta-regression models

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

Species class

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.6924   555.3847   573.3847   614.0438   573.6546   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1085  0.3293    683     no         effectID 
sigma^2.2  0.0081  0.0902    204     no          paperID 
sigma^2.3  0.0383  0.1957    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 677) = 0.5311, p-val = 0.7848

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
species_classActinopterygii    0.0018  0.1742   0.0105  677  0.9916  -0.3402 
species_classArachnida         0.0404  0.2879   0.1402  677  0.8885  -0.5249 
species_classAves             -0.0265  0.0372  -0.7129  677  0.4762  -0.0995 
species_classInsecta           0.0943  0.1500   0.6283  677  0.5300  -0.2003 
species_classMammalia         -0.0652  0.0459  -1.4191  677  0.1563  -0.1554 
species_classReptilia          0.0690  0.1441   0.4791  677  0.6320  -0.2139 
                              ci.ub    
species_classActinopterygii  0.3438    
species_classArachnida       0.6056    
species_classAves            0.0465    
species_classInsecta         0.3888    
species_classMammalia        0.0250    
species_classReptilia        0.3520    

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

r2_ml(mod_class)

   R2_marginal R2_conditional 
   0.006049711    0.303965551 

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 (blue circles) scaled by precision (1/SE).

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).

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

Study type

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.6093   559.2187   569.2187   591.8365   569.3076   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1097  0.3312    683     no         effectID 
sigma^2.2  0.0104  0.1018    204     no          paperID 
sigma^2.3  0.0298  0.1727    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 681) = 0.5914, p-val = 0.5538

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
study_typenatural        -0.0262  0.0308  -0.8491  681  0.3961  -0.0867  0.0344 
study_typesemi-natural   -0.0600  0.0680  -0.8825  681  0.3778  -0.1936  0.0735 
                          
study_typenatural         
study_typesemi-natural    

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

r2_ml(mod_type)

   R2_marginal R2_conditional 
   0.000703997    0.268549933 

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 (blue circles) scaled by precision (1/SE).

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).

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

Study design

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.7171   559.4342   569.4342   592.0520   569.5231   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1097  0.3313    683     no         effectID 
sigma^2.2  0.0122  0.1105    204     no          paperID 
sigma^2.3  0.0270  0.1644    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 681) = 0.5154, p-val = 0.5975

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
study_designcontrast   -0.0240  0.0316  -0.7575  681  0.4490  -0.0861  0.0381 
study_designgradient   -0.0420  0.0458  -0.9173  681  0.3593  -0.1320  0.0480 
                        
study_designcontrast    
study_designgradient    

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

r2_ml(mod_design)

   R2_marginal R2_conditional 
  0.0004010831   0.2636753327 

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 (blue circles) scaled by precision (1/SE).

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).

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

Dispersal type

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.1831   558.3663   570.3663   597.4988   570.4911   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1101  0.3318    683     no         effectID 
sigma^2.2  0.0122  0.1107    204     no          paperID 
sigma^2.3  0.0262  0.1617    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 680) = 0.5919, p-val = 0.6205

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
dispersal_typeB          -0.0912  0.0754  -1.2101  680  0.2266  -0.2393  0.0568 
dispersal_typebreeding   -0.0115  0.0443  -0.2605  680  0.7946  -0.0985  0.0754 
dispersal_typenatal      -0.0248  0.0328  -0.7551  680  0.4505  -0.0891  0.0396 
                          
dispersal_typeB           
dispersal_typebreeding    
dispersal_typenatal       

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

r2_ml(mod_disp)

   R2_marginal R2_conditional 
   0.002077027    0.260143736 

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. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

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 with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

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

Dispersal phase

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.3746   558.7492   568.7492   591.3670   568.8380   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1097  0.3313    683     no         effectID 
sigma^2.2  0.0097  0.0987    204     no          paperID 
sigma^2.3  0.0301  0.1735    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 681) = 1.0002, p-val = 0.3683

Model Results:

                            estimate      se     tval   df    pval    ci.lb 
dispersal_phaseprospection   -0.0651  0.0461  -1.4135  681  0.1580  -0.1556 
dispersal_phasesettlement    -0.0197  0.0316  -0.6240  681  0.5329  -0.0818 
                             ci.ub    
dispersal_phaseprospection  0.0254    
dispersal_phasesettlement   0.0423    

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

r2_ml(mod_phase)

   R2_marginal R2_conditional 
   0.002087416    0.267890548 

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 (blue circles) scaled by precision (1/SE).

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).

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

Sex

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.1389   558.2779   570.2779   597.4104   570.4027   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1096  0.3310    683     no         effectID 
sigma^2.2  0.0122  0.1104    204     no          paperID 
sigma^2.3  0.0273  0.1652    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 680) = 0.9000, p-val = 0.4408

Model Results:

      estimate      se     tval   df    pval    ci.lb   ci.ub    
sexB   -0.0722  0.0446  -1.6201  680  0.1057  -0.1597  0.0153    
sexF   -0.0103  0.0350  -0.2948  680  0.7682  -0.0789  0.0583    
sexM   -0.0148  0.0362  -0.4096  680  0.6823  -0.0860  0.0563    

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

r2_ml(mod_sex)

   R2_marginal R2_conditional 
   0.003419061    0.267391593 

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 (blue circles) scaled by precision (1/SE).

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).

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

Life history stage

3.1 Higher life history grouping

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.0796   554.1592   566.1592   593.2918   566.2840   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3315    683     no         effectID 
sigma^2.2  0.0090  0.0946    204     no          paperID 
sigma^2.3  0.0290  0.1703    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 680) = 2.2474, p-val = 0.0816

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
age_class_cleanadult      -0.0256  0.0313  -0.8169  680  0.4143  -0.0872 
age_class_cleanjuvenile    0.0237  0.0505   0.4687  680  0.6395  -0.0755 
age_class_cleanmix        -0.1549  0.0652  -2.3761  680  0.0178  -0.2828 
                           ci.ub    
age_class_cleanadult      0.0359    
age_class_cleanjuvenile   0.1228    
age_class_cleanmix       -0.0269  * 

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

r2_ml(mod_age1)

   R2_marginal R2_conditional 
    0.00930012     0.26371219 

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 (blue circles) scaled by precision (1/SE).

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).

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

# 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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-275.5914   551.1829   571.1829   616.3448   571.5137   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1102  0.3320    683     no         effectID 
sigma^2.2  0.0093  0.0962    204     no          paperID 
sigma^2.3  0.0284  0.1685    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 676) = 1.2805, p-val = 0.2572

Model Results:

              estimate      se     tval   df    pval    ci.lb    ci.ub    
age_classA     -0.0187  0.0342  -0.5450  676  0.5860  -0.0859   0.0486    
age_classJ      0.0202  0.0510   0.3967  676  0.6917  -0.0799   0.1204    
age_classJA    -0.2970  0.1827  -1.6251  676  0.1046  -0.6557   0.0618    
age_classJY    -0.2084  0.0972  -2.1444  676  0.0324  -0.3993  -0.0176  * 
age_classJYA   -0.0763  0.0922  -0.8272  676  0.4084  -0.2574   0.1048    
age_classY     -0.0445  0.0609  -0.7312  676  0.4649  -0.1641   0.0750    
age_classYA    -0.0498  0.0490  -1.0164  676  0.3098  -0.1460   0.0464    

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

r2_ml(mod_age2)

   R2_marginal R2_conditional 
     0.0138851      0.2649849 

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 (blue circles) scaled by precision (1/SE).

# table not created for this (not meaningful for some levels)

Fitness metric

3.3 Higher level

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.1211   558.2422   570.2422   597.3747   570.3670   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3315    683     no         effectID 
sigma^2.2  0.0115  0.1074    204     no          paperID 
sigma^2.3  0.0271  0.1646    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 680) = 0.9335, p-val = 0.4240

Model Results:

                                      estimate      se     tval   df    pval 
fitness_higher_levelindirect fitness   -0.2352  0.1700  -1.3842  680  0.1668 
fitness_higher_levelreproduction       -0.0163  0.0331  -0.4938  680  0.6216 
fitness_higher_levelsurvival           -0.0358  0.0342  -1.0465  680  0.2957 
                                        ci.lb   ci.ub    
fitness_higher_levelindirect fitness  -0.5689  0.0984    
fitness_higher_levelreproduction      -0.0813  0.0486    
fitness_higher_levelsurvival          -0.1029  0.0313    

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

r2_ml(mod_fit1)

   R2_marginal R2_conditional 
    0.00278045     0.26199113 

orchard_plot(mod_fit1, mod = "fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

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 (blue circles) scaled by precision (1/SE).

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).

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

• need to think what to do with classes less than 5 studies

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.1359   556.2718   578.2718   627.9336   578.6699   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1100  0.3316    683     no         effectID 
sigma^2.2  0.0088  0.0940    204     no          paperID 
sigma^2.3  0.0322  0.1794    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 675) = 0.6761, p-val = 0.7128

Model Results:

                                                           estimate      se 
fitness_metric_cleanage at first reproduction               -0.0272  0.1009 
fitness_metric_cleanindirect fitness                        -0.2300  0.1709 
fitness_metric_cleanlifetime breeding success               -0.0045  0.0386 
fitness_metric_cleanlifetime reproductive success           -0.0314  0.0423 
fitness_metric_cleanoffspring reproduction                   0.0393  0.1067 
fitness_metric_cleanoffspring survival                       0.0034  0.0478 
fitness_metric_cleanreproductive lifespan and/or attempts   -0.0668  0.1021 
fitness_metric_cleansurvival                                -0.0559  0.0381 
                                                              tval   df    pval 
fitness_metric_cleanage at first reproduction              -0.2694  675  0.7877 
fitness_metric_cleanindirect fitness                       -1.3457  675  0.1789 
fitness_metric_cleanlifetime breeding success              -0.1171  675  0.9068 
fitness_metric_cleanlifetime reproductive success          -0.7422  675  0.4582 
fitness_metric_cleanoffspring reproduction                  0.3683  675  0.7128 
fitness_metric_cleanoffspring survival                      0.0701  675  0.9441 
fitness_metric_cleanreproductive lifespan and/or attempts  -0.6535  675  0.5136 
fitness_metric_cleansurvival                               -1.4686  675  0.1424 
                                                             ci.lb   ci.ub    
fitness_metric_cleanage at first reproduction              -0.2253  0.1709    
fitness_metric_cleanindirect fitness                       -0.5655  0.1056    
fitness_metric_cleanlifetime breeding success              -0.0803  0.0713    
fitness_metric_cleanlifetime reproductive success          -0.1146  0.0517    
fitness_metric_cleanoffspring reproduction                 -0.1701  0.2487    
fitness_metric_cleanoffspring survival                     -0.0906  0.0973    
fitness_metric_cleanreproductive lifespan and/or attempts  -0.2673  0.1338    
fitness_metric_cleansurvival                               -0.1307  0.0189    

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

r2_ml(mod_fit2)

   R2_marginal R2_conditional 
   0.006343287    0.276277472 

orchard_plot(mod_fit2, mod = "fitness_metric_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

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 (blue circles) scaled by precision (1/SE).

# table not created for this (not meaningful for some levels)

Multi-generation

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.5152   559.0305   569.0305   591.6483   569.1193   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1093  0.3307    683     no         effectID 
sigma^2.2  0.0090  0.0946    204     no          paperID 
sigma^2.3  0.0324  0.1801    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 681) = 1.0550, p-val = 0.3488

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
whose_fitnessdescendant    0.0064  0.0461   0.1378  681  0.8905  -0.0842 
whose_fitnessfocal        -0.0361  0.0310  -1.1640  681  0.2448  -0.0969 
                          ci.ub    
whose_fitnessdescendant  0.0969    
whose_fitnessfocal       0.0248    

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

r2_ml(mod_gen)

   R2_marginal R2_conditional 
   0.001623982    0.275812105 

orchard_plot(mod_gen, mod = "whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S12. Mean fitness effect of 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 (blue circles) scaled by precision (1/SE).

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).

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

Study duration

3.5 Normal analysis

# 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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.7325   559.4650   569.4650   592.0828   569.5539   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3315    683     no         effectID 
sigma^2.2  0.0116  0.1076    204     no          paperID 
sigma^2.3  0.0276  0.1663    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 681) = 0.0266, p-val = 0.8705

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0199  0.0590  -0.3370  681  0.7362  -0.1358  0.0960    
ln_study_duration   -0.0035  0.0214  -0.1630  681  0.8705  -0.0455  0.0386    

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

r2_ml(mod_dur)

   R2_marginal R2_conditional 
  7.123617e-05   2.631849e-01 

bubble_plot(mod_dur,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S13. Add text.

# 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 = 683; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-11.8954   23.7909   33.7909   56.4087   33.8798   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0480  0.2192    683     no         effectID 
sigma^2.2  0.0074  0.0860    204     no          paperID 
sigma^2.3  0.0125  0.1117    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 681) = 2.7222, p-val = 0.0994

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0145  0.0410  -0.3531  681  0.7242  -0.0951  0.0661    
ln_study_duration    0.0244  0.0148   1.6499  681  0.0994  -0.0046  0.0535  . 

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

r2_ml(mod_dur_abs)

   R2_marginal R2_conditional 
   0.007591265    0.298068402 

bubble_plot(mod_dur_abs,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S14. Add text.

Selection bias

3.6 Normal analysis

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.7453   559.4907   569.4907   592.1085   569.5796   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1097  0.3312    683     no         effectID 
sigma^2.2  0.0119  0.1090    204     no          paperID 
sigma^2.3  0.0277  0.1664    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 681) = 0.4967, p-val = 0.6088

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN   -0.0165  0.0460  -0.3592  681  0.7196  -0.1069  0.0738    
fitness_main_focusY   -0.0315  0.0316  -0.9959  681  0.3196  -0.0935  0.0306    

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

r2_ml(mod_focus)

   R2_marginal R2_conditional 
  0.0002735648   0.2652602138 

orchard_plot(mod_focus, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S17. Add text.

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. Add text.

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

# 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 = 683; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-280.6331  9668.0097   571.2661   593.8986   571.3548   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3315    683     no         effectID 
sigma^2.2  0.0117  0.1080    204     no          paperID 
sigma^2.3  0.0264  0.1625    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 681) = 0.4780, p-val = 0.6202

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN   -0.0163  0.0457  -0.3562  681  0.7218  -0.1060  0.0735    
fitness_main_focusY   -0.0306  0.0313  -0.9772  681  0.3288  -0.0921  0.0309    

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

r2_ml(mod_focus1)

   R2_marginal R2_conditional 
  0.0002536483   0.2575663019 

# 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 = 683; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-280.5161  9667.7758   573.0322   600.1912   573.1565   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0119  0.1089    204     no          paperID 
sigma^2.2  0.0256  0.1601    148     no  species_cleaned 

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

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.1043  0.3230    164     no      N 
tau^2.2    0.1118  0.3344    519     no      Y 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 681) = 0.1066, p-val = 0.7441

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt               -0.0157  0.0451  -0.3483  681  0.7277  -0.1043  0.0729    
fitness_main_focusY   -0.0142  0.0436  -0.3265  681  0.7441  -0.0998  0.0713    

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

#r2_ml(mod_focusA1)

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

        df      AIC      BIC     AICc    logLik    LRT   pval             QE 
Full     6 573.0322 600.1912 573.1565 -280.5161               387470478.4986 
Reduced  5 571.2661 593.8986 571.3548 -280.6331 0.2339 0.6286 387470478.4986 

# homo
orchard_plot(mod_focus1, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S18. Add text.

# hetero
orchard_plot(mod_focus2, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S19. Add text.

Confirmation bias

3.8 Normal analysis

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.1498   558.2997   572.2997   603.9440   572.4666   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3316    683     no         effectID 
sigma^2.2  0.0176  0.1328    204     no          paperID 
sigma^2.3  0.0202  0.1420    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 679) = 0.3253, p-val = 0.8610

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0126  0.0841   0.1502  679  0.8806  -0.1524 
confirmation_biasCost      -0.0389  0.0400  -0.9743  679  0.3303  -0.1174 
confirmation_biasMixed     -0.0273  0.0398  -0.6868  679  0.4924  -0.1054 
confirmation_biasNeutral   -0.0072  0.0470  -0.1539  679  0.8777  -0.0995 
                           ci.ub    
confirmation_biasBenefit  0.1777    
confirmation_biasCost     0.0395    
confirmation_biasMixed    0.0508    
confirmation_biasNeutral  0.0851    

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

r2_ml(mod_confirm)

   R2_marginal R2_conditional 
   0.001341982    0.256948854 

orchard_plot(mod_confirm, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S20. Add text.

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. Add text.

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

# 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 = 683; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-280.4455  9667.6347   574.8911   606.5766   575.0570   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1103  0.3321    683     no         effectID 
sigma^2.2  0.0163  0.1278    204     no          paperID 
sigma^2.3  0.0190  0.1377    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 679) = 0.3141, p-val = 0.8686

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0125  0.0830   0.1503  679  0.8806  -0.1505 
confirmation_biasCost      -0.0380  0.0394  -0.9630  679  0.3359  -0.1153 
confirmation_biasMixed     -0.0258  0.0391  -0.6599  679  0.5095  -0.1025 
confirmation_biasNeutral   -0.0065  0.0463  -0.1403  679  0.8885  -0.0974 
                           ci.ub    
confirmation_biasBenefit  0.1754    
confirmation_biasCost     0.0394    
confirmation_biasMixed    0.0509    
confirmation_biasNeutral  0.0845    

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

r2_ml(mod_confirm1)

   R2_marginal R2_conditional 
   0.001314924    0.243439258 

# 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 = 683; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-242.1490  9591.0417   504.2981   549.5630   504.6255   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0150  0.1225    204     no          paperID 
sigma^2.2  0.0129  0.1135    148     no  species_cleaned 

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

            estim    sqrt  k.lvl  fixed    level 
tau^2.1    0.0296  0.1720     32     no  Benefit 
tau^2.2    0.1326  0.3642    210     no     Cost 
tau^2.3    0.0633  0.2516    289     no    Mixed 
tau^2.4    0.2037  0.4513    152     no  Neutral 

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

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 679) = 0.4496, p-val = 0.7177

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
intrcpt                     0.0236  0.0593   0.3985  679  0.6904  -0.0928 
confirmation_biasCost      -0.0733  0.0647  -1.1335  679  0.2574  -0.2002 
confirmation_biasMixed     -0.0594  0.0616  -0.9649  679  0.3350  -0.1804 
confirmation_biasNeutral   -0.0461  0.0712  -0.6475  679  0.5175  -0.1859 
                           ci.ub    
intrcpt                   0.1400    
confirmation_biasCost     0.0537    
confirmation_biasMixed    0.0615    
confirmation_biasNeutral  0.0937    

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

# they are not significantly different
anova(mod_confirm1, mod_confirm2)

        df      AIC      BIC     AICc    logLik     LRT   pval             QE 
Full    10 504.2981 549.5630 504.6255 -242.1490                387289225.0508 
Reduced  7 574.8911 606.5766 575.0570 -280.4455 76.5930 <.0001 387289225.0508 

# homo
orchard_plot(mod_confirm1, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S21. Add text.

# hetero
orchard_plot(mod_confirm2, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S22. Add text.

4 Interaction meta-regression models (2 moderators)

Sex x Species class

# 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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-273.9193   547.8386   583.8386   664.9158   584.8925   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1097  0.3312    683     no         effectID 
sigma^2.2  0.0115  0.1074    204     no          paperID 
sigma^2.3  0.0295  0.1718    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:15):
F(df1 = 15, df2 = 668) = 0.7704, p-val = 0.7114

Model Results:

                                   estimate      se     tval   df    pval 
sex_species_classB_Aves              0.0003  0.0580   0.0043  668  0.9966 
sex_species_classB_Insecta          -0.1459  0.2746  -0.5314  668  0.5953 
sex_species_classB_Mammalia         -0.1773  0.0668  -2.6535  668  0.0082 
sex_species_classB_Reptilia          0.1592  0.2727   0.5837  668  0.5596 
sex_species_classF_Actinopterygii   -0.0680  0.2061  -0.3298  668  0.7417 
sex_species_classF_Arachnida         0.0412  0.2791   0.1474  668  0.8828 
sex_species_classF_Aves             -0.0345  0.0439  -0.7866  668  0.4318 
sex_species_classF_Insecta           0.2390  0.1859   1.2856  668  0.1990 
sex_species_classF_Mammalia         -0.0125  0.0524  -0.2381  668  0.8119 
sex_species_classF_Reptilia          0.0727  0.1553   0.4682  668  0.6398 
sex_species_classM_Actinopterygii    0.0688  0.2056   0.3348  668  0.7379 
sex_species_classM_Aves             -0.0335  0.0444  -0.7557  668  0.4501 
sex_species_classM_Insecta          -0.1467  0.3773  -0.3888  668  0.6976 
sex_species_classM_Mammalia         -0.0121  0.0568  -0.2123  668  0.8320 
sex_species_classM_Reptilia         -0.0223  0.2312  -0.0964  668  0.9232 
                                     ci.lb    ci.ub     
sex_species_classB_Aves            -0.1136   0.1141     
sex_species_classB_Insecta         -0.6851   0.3933     
sex_species_classB_Mammalia        -0.3086  -0.0461  ** 
sex_species_classB_Reptilia        -0.3762   0.6946     
sex_species_classF_Actinopterygii  -0.4727   0.3368     
sex_species_classF_Arachnida       -0.5069   0.5892     
sex_species_classF_Aves            -0.1207   0.0517     
sex_species_classF_Insecta         -0.1260   0.6041     
sex_species_classF_Mammalia        -0.1154   0.0904     
sex_species_classF_Reptilia        -0.2323   0.3777     
sex_species_classM_Actinopterygii  -0.3348   0.4725     
sex_species_classM_Aves            -0.1207   0.0536     
sex_species_classM_Insecta         -0.8876   0.5942     
sex_species_classM_Mammalia        -0.1236   0.0995     
sex_species_classM_Reptilia        -0.4763   0.4317     

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

r2_ml(mod_sex_species_class)

   R2_marginal R2_conditional 
    0.01839621     0.28577240 

orchard_plot(mod_sex_species_class, mod = "sex_species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S23. Add text.

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. Add text.

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

Sex x Whose fitness

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.9692   555.9384   573.9384   614.5974   574.2082   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1096  0.3311    683     no         effectID 
sigma^2.2  0.0092  0.0958    204     no          paperID 
sigma^2.3  0.0314  0.1771    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 677) = 0.9567, p-val = 0.4538

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
sex_whose_fitnessB_descendant    0.0219  0.1503   0.1458  677  0.8842  -0.2732 
sex_whose_fitnessB_focal        -0.0771  0.0453  -1.7022  677  0.0892  -0.1660 
sex_whose_fitnessF_descendant   -0.0200  0.0547  -0.3658  677  0.7147  -0.1273 
sex_whose_fitnessF_focal        -0.0081  0.0366  -0.2217  677  0.8246  -0.0801 
sex_whose_fitnessM_descendant    0.0620  0.0628   0.9874  677  0.3238  -0.0613 
sex_whose_fitnessM_focal        -0.0318  0.0378  -0.8407  677  0.4008  -0.1060 
                                ci.ub    
sex_whose_fitnessB_descendant  0.3170    
sex_whose_fitnessB_focal       0.0118  . 
sex_whose_fitnessF_descendant  0.0873    
sex_whose_fitnessF_focal       0.0638    
sex_whose_fitnessM_descendant  0.1854    
sex_whose_fitnessM_focal       0.0425    

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

r2_ml(mod_sex_whose_fitness)

   R2_marginal R2_conditional 
   0.006892305    0.275099257 

orchard_plot(mod_sex_whose_fitness, mod = "sex_whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S23. Add text.

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. Add text.

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

Sex x Fitness (higher level)

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.9939   555.9878   577.9878   627.6497   578.3860   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1098  0.3314    683     no         effectID 
sigma^2.2  0.0124  0.1115    204     no          paperID 
sigma^2.3  0.0271  0.1647    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 675) = 0.7235, p-val = 0.6708

Model Results:

                                            estimate      se     tval   df 
sex_fitness_higher_levelB_reproduction       -0.0218  0.0715  -0.3045  675 
sex_fitness_higher_levelB_survival           -0.0935  0.0512  -1.8257  675 
sex_fitness_higher_levelF_indirect fitness   -0.2293  0.2540  -0.9025  675 
sex_fitness_higher_levelF_reproduction       -0.0013  0.0385  -0.0325  675 
sex_fitness_higher_levelF_survival           -0.0226  0.0434  -0.5214  675 
sex_fitness_higher_levelM_indirect fitness   -0.2284  0.2071  -1.1026  675 
sex_fitness_higher_levelM_reproduction       -0.0262  0.0425  -0.6176  675 
sex_fitness_higher_levelM_survival            0.0012  0.0449   0.0265  675 
                                              pval    ci.lb   ci.ub    
sex_fitness_higher_levelB_reproduction      0.7609  -0.1621  0.1186    
sex_fitness_higher_levelB_survival          0.0683  -0.1941  0.0071  . 
sex_fitness_higher_levelF_indirect fitness  0.3671  -0.7280  0.2695    
sex_fitness_higher_levelF_reproduction      0.9741  -0.0769  0.0744    
sex_fitness_higher_levelF_survival          0.6023  -0.1079  0.0626    
sex_fitness_higher_levelM_indirect fitness  0.2706  -0.6351  0.1783    
sex_fitness_higher_levelM_reproduction      0.5371  -0.1097  0.0572    
sex_fitness_higher_levelM_survival          0.9788  -0.0869  0.0893    

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

r2_ml(mod_sex_fitness)

   R2_marginal R2_conditional 
   0.007422826    0.270245437 

orchard_plot(mod_sex_fitness, mod = "sex_fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S24. Add text.

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. Add text.

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

Sex x Dispersal type

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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.1389   556.2778   580.2778   634.4365   580.7498   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1100  0.3317    683     no         effectID 
sigma^2.2  0.0162  0.1272    204     no          paperID 
sigma^2.3  0.0229  0.1512    148     no  species_cleaned 

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

Test of Moderators (coefficients 1:9):
F(df1 = 9, df2 = 674) = 0.4628, p-val = 0.8996

Model Results:

                              estimate      se     tval   df    pval    ci.lb 
sex_dispersal_typeB_B          -0.0938  0.1143  -0.8209  674  0.4120  -0.3183 
sex_dispersal_typeB_breeding   -0.0031  0.0990  -0.0309  674  0.9753  -0.1974 
sex_dispersal_typeB_natal      -0.0822  0.0508  -1.6202  674  0.1057  -0.1819 
sex_dispersal_typeF_B          -0.0814  0.1111  -0.7329  674  0.4639  -0.2996 
sex_dispersal_typeF_breeding   -0.0022  0.0535  -0.0403  674  0.9679  -0.1072 
sex_dispersal_typeF_natal      -0.0043  0.0389  -0.1104  674  0.9121  -0.0807 
sex_dispersal_typeM_B          -0.1086  0.1225  -0.8863  674  0.3758  -0.3490 
sex_dispersal_typeM_breeding   -0.0204  0.0584  -0.3496  674  0.7267  -0.1351 
sex_dispersal_typeM_natal      -0.0018  0.0411  -0.0438  674  0.9650  -0.0825 
                               ci.ub    
sex_dispersal_typeB_B         0.1306    
sex_dispersal_typeB_breeding  0.1912    
sex_dispersal_typeB_natal     0.0174    
sex_dispersal_typeF_B         0.1368    
sex_dispersal_typeF_breeding  0.1029    
sex_dispersal_typeF_natal     0.0721    
sex_dispersal_typeM_B         0.1319    
sex_dispersal_typeM_breeding  0.0942    
sex_dispersal_typeM_natal     0.0789    

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

r2_ml(mod_sex_dispersal_type)

   R2_marginal R2_conditional 
    0.00672015     0.26697673 

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 S25. Add text.

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. Add text.

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 Multi-moderator-regression model

Full model

# 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 = 683; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-274.8126  9656.3688   577.6252   640.9961   578.2539   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1083  0.3290    683     no         effectID 
sigma^2.2  0.0079  0.0891    204     no          paperID 
sigma^2.3  0.0305  0.1745    148     no  species_cleaned 

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

Test of Moderators (coefficients 2:11):
F(df1 = 10, df2 = 672) = 1.1939, p-val = 0.2915

Model Results:

                                  estimate      se     tval   df    pval 
intrcpt                            -0.2619  0.1980  -1.3228  672  0.1864 
sexF                                0.0464  0.0499   0.9286  672  0.3534 
sexM                                0.0427  0.0507   0.8424  672  0.3999 
age_class_cleanjuvenile             0.1033  0.0574   1.7981  672  0.0726 
age_class_cleanmix                 -0.0878  0.0679  -1.2916  672  0.1969 
whose_fitnessfocal                 -0.0506  0.0481  -1.0521  672  0.2931 
fitness_higher_levelreproduction    0.2121  0.1699   1.2483  672  0.2124 
fitness_higher_levelsurvival        0.1814  0.1725   1.0512  672  0.2935 
dispersal_typebreeding              0.0214  0.0842   0.2544  672  0.7993 
dispersal_typenatal                 0.0062  0.0770   0.0803  672  0.9360 
dispersal_phasesettlement           0.0378  0.0535   0.7075  672  0.4795 
                                    ci.lb   ci.ub    
intrcpt                           -0.6507  0.1269    
sexF                              -0.0517  0.1444    
sexM                              -0.0568  0.1422    
age_class_cleanjuvenile           -0.0095  0.2161  . 
age_class_cleanmix                -0.2212  0.0457    
whose_fitnessfocal                -0.1449  0.0438    
fitness_higher_levelreproduction  -0.1215  0.5457    
fitness_higher_levelsurvival      -0.1574  0.5201    
dispersal_typebreeding            -0.1440  0.1868    
dispersal_typenatal               -0.1450  0.1573    
dispersal_phasesettlement         -0.0672  0.1429    

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

round(r2_ml(mod_full)*100, 2)

   R2_marginal R2_conditional 
          1.96          27.63 

Multi-model inference

# multi-model selection
candidates <- dredge(mod_full, trace = 2)

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

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) 

# 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

# 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          

6 Publication bias

Egger regression: uni-moderator

# 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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.8086   559.6171   569.6171   592.2350   569.7060   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1100  0.3317    683     no         effectID 
sigma^2.2  0.0120  0.1095    204     no          paperID 
sigma^2.3  0.0267  0.1633    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 681) = 0.0040, p-val = 0.9494

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -0.0294  0.0401  -0.7330  681  0.4638  -0.1082  0.0494    
sqeffectN    0.0003  0.0045   0.0635  681  0.9494  -0.0085  0.0091    

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

round(r2_ml(mod_egger)*100, 2)

   R2_marginal R2_conditional 
          0.00          26.01 

small <- bubble_plot(mod_egger,
                     mod = "sqeffectN",
                     group = "paperID",
                     xlab = "sqrt(Effective N)",
                     ylab = "Effect Size: Zr",
                     g = TRUE)

#| fig-cap: "**Figure S27.** Add text."

small

Decline effect: uni-moderator

# 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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-224.1660   448.3319   458.3319   480.9497   458.4208   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0631  0.2511    683     no         effectID 
sigma^2.2  0.0125  0.1120    204     no          paperID 
sigma^2.3  0.0263  0.1622    148     no  species_cleaned 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 681) = 0.6145, p-val = 0.4334

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt   -2.4143  3.0537  -0.7906  681  0.4294  -8.4101  3.5815    
year       0.0012  0.0015   0.7839  681  0.4334  -0.0018  0.0042    

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

round(r2_ml(mod_dec)*100, 2)

   R2_marginal R2_conditional 
          0.22          38.26 

decline <- bubble_plot(mod_dec,
                       mod = "year",
                       group = "paperID",
                       xlab = "Publication year",
                       ylab = "Effect Size: Zr",
                       g = TRUE)

decline 

Figure S28. Add text.

Combined analysis

# 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 = 683; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.3441   558.6881   570.6881   597.8207   570.8129   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1096  0.3310    683     no         effectID 
sigma^2.2  0.0112  0.1057    204     no          paperID 
sigma^2.3  0.0294  0.1713    148     no  species_cleaned 

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

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 680) = 0.1896, p-val = 0.8273

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -2.0823  3.3446  -0.6226  680  0.5338  -8.6493  4.4846    
year         0.0010  0.0017   0.6136  680  0.5397  -0.0023  0.0043    
sqeffectN   -0.0002  0.0045  -0.0461  680  0.9633  -0.0091  0.0087    

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

round(r2_ml(mod_comb)*100, 2)

   R2_marginal R2_conditional 
          0.11          27.08 

# small-study
bubble_plot(mod_egger,
            mod = "sqeffectN",
            group = "paperID",
            xlab = "sqrt(Effective N)",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S29. Add text.

# decline
bubble_plot(mod_comb,
            mod = "year",
            group = "paperID",
            xlab = "Publication year",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S30. Add text.

Leave-1study-out (sensitivity analysis)

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, 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"))

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))

leaveoneout

Figure S31. Add text.

7 Figure gallery

Figure 1

reduced_intercept

Test caption

Figure —

8 R Session Information

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

R version 4.2.0 (2022-04-22)

Platform: x86_64-apple-darwin17.0 (64-bit)

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: matrixcalc(v.1.0-6), orchaRd(v.2.0), here(v.1.0.1), patchwork(v.1.1.3), ape(v.5.7-1), metafor(v.4.4-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.5-4.1), pander(v.0.6.5), gridExtra(v.2.3), emmeans(v.1.8.9), MuMIn(v.1.47.5), kableExtra(v.1.3.4), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.3), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.4) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): httr(v.1.4.7), jsonlite(v.1.8.7), viridisLite(v.0.4.2), splines(v.4.2.0), highr(v.0.10), stats4(v.4.2.0), vipor(v.0.4.7), yaml(v.2.3.7), pillar(v.1.9.0), lattice(v.0.22-5), glue(v.1.6.2), digest(v.0.6.33), rvest(v.1.0.3), colorspace(v.2.1-0), sandwich(v.3.1-0), htmltools(v.0.5.6.1), pkgconfig(v.2.0.3), xtable(v.1.8-4), mvtnorm(v.1.2-3), scales(v.1.3.0), webshot(v.0.5.5), svglite(v.2.1.2), tzdb(v.0.4.0), timechange(v.0.2.0), mgcv(v.1.9-0), farver(v.2.1.1), generics(v.0.1.3), TH.data(v.1.1-2), pacman(v.0.5.1), withr(v.2.5.2), cli(v.3.6.1), survival(v.3.5-7), magrittr(v.2.0.3), estimability(v.1.4.1), evaluate(v.0.23), fansi(v.1.0.5), nlme(v.3.1-163), MASS(v.7.3-60), xml2(v.1.3.5), beeswarm(v.0.4.0), tools(v.4.2.0), hms(v.1.1.3), lifecycle(v.1.0.4), multcomp(v.1.4-25), munsell(v.0.5.0), compiler(v.4.2.0), systemfonts(v.1.0.5), rlang(v.1.1.1), grid(v.4.2.0), rstudioapi(v.0.15.0), htmlwidgets(v.1.6.1), labeling(v.0.4.3), rmarkdown(v.2.25), gtable(v.0.3.4), codetools(v.0.2-19), R6(v.2.5.1), zoo(v.1.8-12), knitr(v.1.45), fastmap(v.1.1.1), utf8(v.1.2.3), rprojroot(v.2.0.4), mathjaxr(v.1.6-0), latex2exp(v.0.9.6), ggbeeswarm(v.0.7.2), stringi(v.1.7.12), parallel(v.4.2.0), Rcpp(v.1.0.11), vctrs(v.0.6.4), tidyselect(v.1.2.0), xfun(v.0.40) and coda(v.0.19-4)