Systematic review and meta-analysis of anti-predator mechanism of eyespots: conspicuous pattern vs eye mimicry

1 Summary

This page contains the all code and data of Systematic review and meta-analysis of eyespot anti-predator mechanisms: conspicuous patterns vs. eye mimicry. Our meta-analytical study examined why butterfly eyespots intimidate birds and tested two hypotheses: the eye mimicry hypothesis and the conspicuousness hypothesis. Despite the existence of conflicting past empirical research, our results supported the idea of the conspicuousness hypothesis. We also found a key factor: larger eyespots were more effective in predator avoidance. Our findings suggest that birds avoid eyespots mainly due to their conspicuousness and emphasise the importance of such pattern characteristics in predator avoidance.

2 Setting ups

2.1 Load packages and custom function

Load packages, custom function, and dataset. Custom function is used for calculating the point estimate and sampling variance (var) of lnRR can be then calculated from raw data.

Code

pacman::p_load(ape,
               GoodmanKruskal,
               DT,
               dtplyr,
               ggokabeito,
               ggcorrplot,
               here, 
               knitr,
               tidyverse,
               patchwork,
               metafor,
               miWQS,
               multcomp,
               orchaRd,
               parallel,
               RColorBrewer)

# function for calculating effect size (lnRR) - the same as function.R
effect_lnRR <- function(dt) {
  # replace 0 in "C_mean", "T_sd", "C_sd", "C_proportion" with each minimum values
  # if proportion too extreme value, replace minimum value (only "T_proportion")

  dt <- dt %>%
  mutate(
    C_mean = ifelse(C_mean == 0, 0.04, C_mean),
    T_sd = ifelse(T_sd == 0, 0.01, T_sd),
    C_sd = ifelse(C_sd == 0, 0.05, C_sd),
    C_proportion = ifelse(C_proportion == 0, 0.01, C_proportion),
    T_proportion = ifelse(T_proportion < 0.01, 0.01, T_proportion),
    T_proportion = ifelse(T_proportion == 1, 0.9, T_proportion)
  )
    
  # copy dataset for adding effect size and its variation (lnRR /lnRR_var) column
  dt1 <- dt %>%
    mutate(
      lnRR = NA,
      lnRR_var = NA
    )

  for (i in seq_len(nrow(dt1))) {
    Tn <- dt1$Tn[i]
    Cn <- dt1$Cn[i]
    T_mean <- dt1$T_mean[i]
    C_mean <- dt1$C_mean[i]
    T_proportion <- dt1$T_proportion[i]
    C_proportion <- dt1$C_proportion[i]
    T_sd <- dt1$T_sd[i]
    C_sd <- dt1$C_sd[i]
    Response <- dt1$Response[i]
    Measure <- dt1$Measure[i]
    Study_design <- dt1$Study_design[i]
    Direction <- dt1$Direction[i]

    # continuous data - using escalc() (metafor package)
    if (Response == "continuous" & Study_design == "independent") {

      effect <- escalc(
        measure = "ROM",
        n1i = Tn, n2i = Cn,
        m1i = T_mean, m2 = C_mean,
        sd1i = T_sd, sd2i = C_sd,
        var.names = c("lnRR", "lnRR_var")
      )

      dt1$lnRR[i] <- effect$lnRR
      dt1$lnRR_var[i] <- effect$lnRR_var
    } else if (Response == "continuous" & Study_design == "dependent") {

      effect <- escalc(
        measure = "ROMC",
        ni = (Tn + Cn) / 2,
        m1i = T_mean, m2 = C_mean,
        sd1i = T_sd, sd2i = C_sd,
        ri = 0.5,
        var.names = c("lnRR", "lnRR_var")
      )

      dt1$lnRR[i] <- effect$lnRR
      dt1$lnRR_var[i] <- effect$lnRR_var
    }

    # proportion data (no sd values)
    else if (Response == "proportion1" & Study_design == "independent") {
      
    T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )

      # transform proportion value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro1 <- log(T_proportion / C_proportion)
      lnRR_var_pro1 <- (1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn) +
        1 / (C_proportion^2 * Cn))

      dt1$lnRR[i] <- lnRR_pro1
      dt1$lnRR_var[i] <- lnRR_var_pro1

    } else if (Response == "proportion1" & Study_design == "dependent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )

      # transform proportion value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro1 <- log(T_proportion / C_proportion)
      lnRR_var_pro1 <- (1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn) +
        1 / (C_proportion^2 * Cn)) -
        2 * 0.5 * sqrt((1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn))) *
          sqrt((1 / sqrt(8))^2 * (1 / (C_proportion^2 * Cn)))

      dt1$lnRR[i] <- lnRR_pro1
      dt1$lnRR_var[i] <- lnRR_var_pro1
    }

    # proportion data (exist sd values)
    else if (Response == "proportion2" & Study_design == "independent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )
      
      # transform proportion mean value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_SD <- sqrt(T_sd^2 / (4 * (T_proportion) * (1 - (T_proportion))))
      C_SD <- sqrt(C_sd^2 / (4 * (C_proportion) * (1 - (C_proportion))))

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro2 <- log(T_proportion / C_proportion)
      lnRR_var_pro2 <- (T_SD)^2 * (1 / (T_proportion^2 * Tn)) +
        (C_SD)^2 * (1 / (C_proportion^2 * Cn))

      dt1$lnRR[i] <- lnRR_pro2
      dt1$lnRR_var[i] <- lnRR_var_pro2

    } else if (Response == "proportion2" & Study_design == "dependent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )
      
      # transform proportion mean value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_SD <- sqrt(T_sd^2 / (4 * (T_proportion) * (1 - (T_proportion))))
      C_SD <- sqrt(C_sd^2 / (4 * (C_proportion) * (1 - (C_proportion))))

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro2 <- log(T_proportion / C_proportion)
      lnRR_var_pro2 <- (T_SD)^2 * (1 / (T_proportion^2 * Tn)) +
        (C_SD)^2 * (1 / (C_proportion^2 * Cn)) -
        2 * 0.5 * sqrt((T_SD)^2 * (1 / (T_proportion^2 * Tn))) *
          sqrt((C_SD)^2 * (1 / (C_proportion^2 * Cn)))

      dt1$lnRR[i] <- lnRR_pro2
      dt1$lnRR_var[i] <- lnRR_var_pro2
    }
  }

  return(dt1)
}

# setting tables
options(DT.options = list(dom = "Blfrtip",
                          scrollX = TRUE,
                          pageLength = 5,
                          columnDefs = list(list(targets = '_all', 
                                                 className = 'dt-center')),
                          buttons = c('copy', 'csv', 'excel', 'pdf')))

2.2 Dataset for analysis

Dataset

We calculated lnRR and lnRR variance and prepared the dataset for analysis. After that, we added effect size ID (obs_ID), calculated the lnRR and variance-covariance matrix (VCV), and log-transformed moderators.

Code

# read raw data()
dat_all <-  read.csv(here("data/all.csv"), header = T, fileEncoding = "CP932")

# calculate effect sizes and their variances
dt_all <- effect_lnRR(dat_all)

# add effect size ID (obs_ID)
dt_all$Obs_ID <- 1:nrow(dt_all)

# use vcalc to calculate variance-covariance matrix
VCV <- vcalc(vi = lnRR_var, 
             cluster = Cohort_ID,
             obs = Obs_ID,
             rho = 0.5,
             data = dt_all)

# convert pattern maximum diameter/length, area total area to logarithm
dt_all$area_pattern_T  <- dt_all$Area_pattern * dt_all$Number_pattern
dt_all$Log_diameter <- log(dt_all$Diameter_pattern)
dt_all$Log_area <- log(dt_all$Area_pattern)
dt_all$Log_background <- log(dt_all$Area_background)
dt_all$Log_area_pattern_T  <- log(dt_all$Area_pattern * dt_all$Number_pattern)

dt_all <- dt_all %>%
  mutate_if(is.character, as.factor)

Effect size distribution

These are the distributions of effect size and its variance.

Code

hist(dt_all$lnRR)
hist(dt_all$lnRR_var)

Figure 1— lnRR

Figure 2— lnRR variance

Description of dataset

3 Let’s meta-analysis and meta-regressions

3.1 Meta-analysis: overall effect size

First, we run meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins et al. 2003) and the partial I2 explained by each random factor.

Code

ma_all <- rma.mv(yi = lnRR,
                  V = VCV, 
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(ma_all)

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

   logLik   Deviance        AIC        BIC       AICc   
-251.0417   502.0833   512.0833   530.0569   512.3115   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0485  0.2203     33     no           Study_ID 
sigma^2.2  0.0000  0.0000     89     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    164     no          Cohort_ID 
sigma^2.4  0.2334  0.4831    270     no             Obs_ID 

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

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.1977  0.0581  3.4039  269  0.0008  0.0833  0.3120  *** 

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

Code

i2 <- round(i2_ml(ma_all), 4)
r2 <- round(r2_ml(ma_all), 4)

I2_Total	I2_Study_ID	I2_Shared_control_ID	I2_Cohort_ID	I2_Obs_ID
96.4974	16.6131	0	0	79.8844

R2_marginal	R2_conditional
0	0.1722

Which random effects are removed in our models?

We can remove Shared control ID and Cohort ID

Code

orchard_plot(ma_all,
              group = "Study_ID",
              xlab = "log response ratio (lnRR)", 
              angle = 45) + 
              scale_x_discrete(labels = c("Overall effect")) + 
              scale_colour_brewer(palette = "Accent") +
              scale_fill_brewer(palette = "Accent") 

Figure 3— Estimates of overall effect size and 95% confidence intervals

3.2 Meta-regressions: uni-moderator

Next, we performed uni-moderator meta-regression models with each of seven moderators: 1) treatment stimulus pattern types (eyespots vs. non-eyespots), 2) pattern area, 3) the number of pattern shapes, 4) prey material type, 5) maximum pattern diameter/length, 6) total pattern area, 7) total area of prey surface, and 8) prey shape type.

We used log-transformed data for pattern area, total pattern area, total area of prey surface, and pattern maximum diameter/length in our analysis to normalise these moderators.

1. Pattern type

Eyespot vs. non-eyespot patterns - is there a significant difference of effect size between two types of patterns?

3.2.0.1 Normal model

Code

mr_eyespot <- rma.mv(yi = lnRR,
                    V = VCV, 
                    mods = ~ Treatment_stimulus,
                    random = list(~1 | Study_ID,
                                  ~1 | Obs_ID),
                    test = "t",
                    method = "REML", 
                    sparse = TRUE,
                    data = dt_all)

summary(mr_eyespot)

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

   logLik   Deviance        AIC        BIC       AICc   
-250.1918   500.3835   508.3835   522.7475   508.5356   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0526  0.2293     33     no  Study_ID 
sigma^2.2  0.2335  0.4833    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.3291, p-val = 0.5666

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
intrcpt                          0.2181  0.0688   3.1722  268  0.0017   0.0827 
Treatment_stimulusnon_eyespot   -0.0603  0.1050  -0.5737  268  0.5666  -0.2671 
                                ci.ub     
intrcpt                        0.3535  ** 
Treatment_stimulusnon_eyespot  0.1465     

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

Code

r2_1 <- round(r2_ml(mr_eyespot), 4)

R2_marginal	R2_conditional
0.0027	0.186

3.2.0.2 Intercept-removed model

Code

# intercept-removed model 
mr_eyespot1 <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~ Treatment_stimulus - 1,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_all)

summary(mr_eyespot1)

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

   logLik   Deviance        AIC        BIC       AICc   
-250.1918   500.3835   508.3835   522.7475   508.5356   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0526  0.2293     33     no  Study_ID 
sigma^2.2  0.2335  0.4833    270     no    Obs_ID 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 268) = 5.7321, p-val = 0.0037

Model Results:

                               estimate      se    tval   df    pval    ci.lb 
Treatment_stimuluseyespot        0.2181  0.0688  3.1722  268  0.0017   0.0827 
Treatment_stimulusnon_eyespot    0.1579  0.0921  1.7132  268  0.0878  -0.0236 
                                ci.ub     
Treatment_stimuluseyespot      0.3535  ** 
Treatment_stimulusnon_eyespot  0.3393   . 

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

Code

orchard_plot(mr_eyespot,
            mod = "Treatment_stimulus",
            group = "Study_ID",
            xlab = "log response ratio (lnRR)",
            angle = 45) + 
            scale_colour_brewer(palette = "Set2") +
            scale_fill_brewer(palette = "Set2")

Figure 4— Eyespot vs. non-eyespot patterns

2. Pattern area (mm²)

Does a larger pattern area improve predator avoidance?

Code

mr_area <- rma.mv(yi = lnRR,
                  V = VCV, 
                  mods = ~ Log_area,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(mr_area)

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

   logLik   Deviance        AIC        BIC       AICc   
-246.5105   493.0210   501.0210   515.3849   501.1731   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0472  0.2172     33     no  Study_ID 
sigma^2.2  0.2265  0.4759    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 7.3561, p-val = 0.0071

Model Results:

          estimate      se     tval   df    pval    ci.lb   ci.ub     
intrcpt    -0.1972  0.1563  -1.2616  268  0.2082  -0.5049  0.1105     
Log_area    0.1119  0.0412   2.7122  268  0.0071   0.0307  0.1931  ** 

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

Code

r2_2 <- round(r2_ml(mr_area), 4)

R2_marginal	R2_conditional
0.0856	0.2433

Code

bubble_plot(mr_area,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Figure 5— Effect size and each pattern area

3. The number of pattern shapes

Does the number of pattern shapes affect predator avoidance?

Code

mr_num <- rma.mv(yi = lnRR,
                  V = VCV,
                  mods = ~ Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

summary(mr_num)

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

   logLik   Deviance        AIC        BIC       AICc   
-248.0452   496.0904   504.0904   518.4544   504.2425   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0442  0.2102     33     no  Study_ID 
sigma^2.2  0.2296  0.4792    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 5.2375, p-val = 0.0229

Model Results:

                estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt           0.3437  0.0854   4.0229  268  <.0001   0.1755   0.5119  *** 
Number_pattern   -0.0584  0.0255  -2.2886  268  0.0229  -0.1086  -0.0082    * 

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

Code

r2_3 <- round(r2_ml(mr_num), 4)

R2_marginal	R2_conditional
0.0246	0.182

Code

bubble_plot(mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") 

Figure 6— Effect size and number of patterns

4. Prey material type

Our dataset includes two material types of prey: real/imitation and artificial abstract butterflies. Is there a significant difference in effect size between the two stimuli types?

3.2.0.3 Normal model

Code

mr_prey_type <- rma.mv(yi = lnRR,
                  V = VCV, 
                  mods = ~ Type_prey,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(mr_prey_type)

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

   logLik   Deviance        AIC        BIC       AICc   
-250.2013   500.4026   508.4026   522.7666   508.5547   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0552  0.2351     33     no  Study_ID 
sigma^2.2  0.2329  0.4826    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.1246, p-val = 0.7244

Model Results:

               estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt          0.1829  0.0749  2.4401  268  0.0153   0.0353  0.3304  * 
Type_preyreal    0.0444  0.1259  0.3529  268  0.7244  -0.2035  0.2924    

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

Code

r2_4 <- round(r2_ml(mr_prey_type), 4)

R2_marginal	R2_conditional
0.0013	0.1928

Code

orchard_plot(mr_prey_type,
              mod = "Type_prey",
              group = "Study_ID",
              xlab = "Type of prey",
              angle = 45) +
              scale_colour_brewer(palette = "Set3") +
              scale_fill_brewer(palette = "Set3")

Figure 7— Effect size and prey type - real vs. abstract

3.2.0.4 Intercept-removed model

Code

mr_prey_type1 <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Type_prey -1,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                        test = "t",
                        method = "REML", 
                        sparse = TRUE,
                        data = dt_all)

summary(mr_prey_type1)

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

   logLik   Deviance        AIC        BIC       AICc   
-250.2013   500.4026   508.4026   522.7666   508.5547   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0552  0.2351     33     no  Study_ID 
sigma^2.2  0.2329  0.4826    270     no    Obs_ID 

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

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 268) = 5.4986, p-val = 0.0046

Model Results:

                     estimate      se    tval   df    pval   ci.lb   ci.ub    
Type_preyartificial    0.1829  0.0749  2.4401  268  0.0153  0.0353  0.3304  * 
Type_preyreal          0.2273  0.1012  2.2457  268  0.0255  0.0280  0.4266  * 

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

5. Diameter/length (mm)

Does a larger pattern diameter/length improve predator avoidance?

Code

mr_diameter <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~ Log_diameter,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_all)

summary(mr_diameter)

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

   logLik   Deviance        AIC        BIC       AICc   
-247.2115   494.4230   502.4230   516.7869   502.5751   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0440  0.2099     33     no  Study_ID 
sigma^2.2  0.2286  0.4781    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 6.0336, p-val = 0.0147

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt        -0.1757  0.1618  -1.0858  268  0.2785  -0.4942  0.1428    
Log_diameter    0.1931  0.0786   2.4563  268  0.0147   0.0383  0.3479  * 

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

Code

r2_5 <- round(r2_ml(mr_diameter), 4)

R2_marginal	R2_conditional
0.0662	0.2171

Code

bubble_plot(mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

Figure 8— Effect size and pattern diameter

6. Total pattern area (mm²)

Does total pattern area affect predator avoidance?

Code

mr_area_T  <- rma.mv(yi = lnRR,
                  V = VCV, 
                  mods = ~ Log_area_pattern_T ,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(mr_area_T)

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

   logLik   Deviance        AIC        BIC       AICc   
-248.0180   496.0359   504.0359   518.3999   504.1880   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0498  0.2231     33     no  Study_ID 
sigma^2.2  0.2290  0.4786    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 4.2763, p-val = 0.0396

Model Results:

                    estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt              -0.1783  0.1909  -0.9338  268  0.3512  -0.5542  0.1976    
Log_area_pattern_T    0.0872  0.0422   2.0679  268  0.0396   0.0042  0.1703  * 

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

Code

r2_6 <- round(r2_ml(mr_area_T), 4)

R2_marginal	R2_conditional
0.0518	0.2211

Code

bubble_plot(mr_area_T,
            mod = "Log_area_pattern_T",
            group = "Study_ID",
            xlab = "Log-transformed pattern total area")

{#fig-total area width=672}

7. Total area of prey surface (mm²)

Does the total prey surface area affect predator avoidance?

Code

mr_background <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Log_background,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                        test = "t",
                        method = "REML",
                        sparse = TRUE,
                        data = dt_all)

summary(mr_background)

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

   logLik   Deviance        AIC        BIC       AICc   
-249.8534   499.7068   507.7068   522.0707   507.8589   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0545  0.2334     33     no  Study_ID 
sigma^2.2  0.2331  0.4828    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.2291, p-val = 0.6326

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt           0.3997  0.4246   0.9413  268  0.3474  -0.4363  1.2356    
Log_background   -0.0289  0.0604  -0.4786  268  0.6326  -0.1479  0.0900    

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

Code

r2_7 <- round(r2_ml(mr_background), 4)

R2_marginal	R2_conditional
0.0042	0.1928

Code

bubble_plot(mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed total prey surface area")

Figure 9— Effect size and total prey surface area

8. Prey shape type

Our dataset includes four types of prey type: real butterfly, abstract artificial butterfly, abstract artificial caterpillar, and abstract artificial prey. Is there a significant difference in effect size between these prey?

Code

mr_prey_shape <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Shape_prey -1,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                                      test = "t",
                                      method = "REML", 
                                      sparse = TRUE,
                                      data = dt_all)

summary(mr_prey_shape)

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

   logLik   Deviance        AIC        BIC       AICc   
-246.3898   492.7797   504.7797   526.2806   505.1040   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0537  0.2317     33     no  Study_ID 
sigma^2.2  0.2305  0.4801    270     no    Obs_ID 

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

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 266) = 3.6583, p-val = 0.0064

Model Results:

                                estimate      se    tval   df    pval    ci.lb 
Shape_preyabstract_butterfly      0.3219  0.1057  3.0446  266  0.0026   0.1137 
Shape_preyabstract_caterpillar    0.0661  0.1259  0.5246  266  0.6003  -0.1819 
Shape_preyabstract_prey           0.0113  0.1842  0.0615  266  0.9510  -0.3514 
Shape_preybutterfly               0.2263  0.1004  2.2549  266  0.0250   0.0287 
                                 ci.ub     
Shape_preyabstract_butterfly    0.5301  ** 
Shape_preyabstract_caterpillar  0.3140     
Shape_preyabstract_prey         0.3741     
Shape_preybutterfly             0.4239   * 

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

Code

r2_8 <- round(r2_ml(mr_prey_shape), 4)

dat_all$Shape_prey <- factor(dat_all$Shape_prey)
levels(dat_all$Shape_prey)

[1] "abstract_butterfly"   "abstract_caterpillar" "abstract_prey"       
[4] "butterfly"           

Code

mat_ex <- cbind(contrMat(rep(1,length(mr_prey_shape$ci.ub)), type = "Tukey"))
sig_test <- summary(glht(mr_prey_shape, linfct= mat_ex), test = adjusted("none"))

sig_test

     Simultaneous Tests for General Linear Hypotheses

Fit: rma.mv(yi = lnRR, V = VCV, mods = ~Shape_prey - 1, random = list(~1 | 
    Study_ID, ~1 | Obs_ID), data = dt_all, method = "REML", test = "t", 
    sparse = TRUE)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0 -0.25583    0.16444  -1.556    0.120
3 - 1 == 0 -0.31057    0.21242  -1.462    0.144
4 - 1 == 0 -0.09561    0.14577  -0.656    0.512
3 - 2 == 0 -0.05474    0.22317  -0.245    0.806
4 - 2 == 0  0.16022    0.16104   0.995    0.320
4 - 3 == 0  0.21496    0.20980   1.025    0.306
(Adjusted p values reported -- none method)

R2_marginal	R2_conditional
0.0521	0.2312

Code

orchard_plot(mr_prey_shape,
              mod = "Shape_prey",
              group = "Study_ID",
              xlab = "Shape of prey",
              angle = 45) +
              scale_colour_brewer(palette = "Set3") +
              scale_fill_brewer(palette = "Set3") 

Figure 10— Effect size and prey shape

3.3 Correlation visualisation and choosing moderators

Before we run multi-moderator meta-regressions, we need to consider the correlation between moderators. Area, diameter and background seem to be correlated. Therefore, we visualised the correlation between these variables.

Visualisation of correlation

3.3.0.1 Correlation between continuous variables

Code

corr_cont <- round(cor(dt_all[, c("Diameter_pattern", "Area_pattern", "Number_pattern","Log_area_pattern_T", "Area_background")]),2)

p_cont <- ggcorrplot(corr_cont, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(a) Continuous variables")

corr_cont_log <- round(cor(dt_all[, c("Log_diameter", "Log_area", "Number_pattern","Log_area_pattern_T", "Log_background")]),2)

p_cont_log <- ggcorrplot(corr_cont_log, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(b) Log-transormed continuous variables")

p_cont + p_cont_log + plot_layout(guides = 'collect')

Figure 11— Correlation between coninuous moderators

All were correlated except the number of pattern shapes.

3.3.0.2 Correlation between continuous variables

Code

dat1 <- dt_all %>%
  dplyr::select("Treatment_stimulus", "Type_prey", "Shape_prey")

corr_cat <- GKtauDataframe(dat1)
plot(corr_cat)

Figure 12— Correlation between categorical variables

We should not include “Shape_prey (prey material type)” and “Type_prey (prey shape type)” in the same model as they are correlated.

Choose moderators

We used model R2 values to find better models and moderator VIF values to check multicollinearity between moderators. Higher R2 indicates a better model, and VIF > 2 indicates multicollinearity.

R2

Code

r2_area <- rma.mv(yi = lnRR,
                  V = VCV,
                  mods = ~ Log_area + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_all)

r2_area_T <- rma.mv(yi = lnRR,
                  V = VCV,
                  mods = ~ Log_area_pattern_T + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_all)

r2_diameter <- rma.mv(yi = lnRR,
                      V = VCV,
                      mods = ~ Log_diameter + Number_pattern,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

# check r2 values
r2_ml(r2_area)

   R2_marginal R2_conditional 
    0.08098982     0.22938197 

Code

r2_ml(r2_area_T)

   R2_marginal R2_conditional 
    0.08546244     0.23909238 

Code

r2_ml(r2_diameter)

   R2_marginal R2_conditional 
    0.06846176     0.20934212 

It seems pattern area is a slightly better predictor than diameter.

3.4 Meta-regressions: multi-moderator

Third, we also ran a multi-moderator meta-regression model, including treatment stimulus pattern types, pattern area, the number of pattern shapes, and prey material type variables due to moderator correlations.

Code

mr_all <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Treatment_stimulus + Log_area 
                + Number_pattern + Type_prey,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_all)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.6478   485.2957   499.2957   524.3538   499.7315   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0522  0.2285     33     no  Study_ID 
sigma^2.2  0.2244  0.4738    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 265) = 2.8780, p-val = 0.0233

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
intrcpt                         -0.0635  0.2073  -0.3065  265  0.7595  -0.4716 
Treatment_stimulusnon_eyespot   -0.0216  0.1065  -0.2027  265  0.8395  -0.2312 
Log_area                         0.0964  0.0446   2.1612  265  0.0316   0.0086 
Number_pattern                  -0.0527  0.0287  -1.8393  265  0.0670  -0.1092 
Type_preyreal                    0.1776  0.1358   1.3080  265  0.1920  -0.0898 
                                ci.ub    
intrcpt                        0.3446    
Treatment_stimulusnon_eyespot  0.1881    
Log_area                       0.1841  * 
Number_pattern                 0.0037  . 
Type_preyreal                  0.4450    

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

Code

r2_9 <- round(r2_ml(mr_all), 4)

R2_marginal	R2_conditional
0.0833	0.2563

3.5 Publication bias

Finally, we tested for publication bias using a funnel plot and Egger’s test. We also checked the decline effect by adding the year of publication as a moderator to the meta-analysis model.

Funnel plot

Code

# funnel plot - standard error
funnel(ma_all, yaxis = "sei",
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)" )

# funnel plot - inverse of standard error
funnel(ma_all, yaxis = "seinv",
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",  col = c(alpha("orange", 0.5)))

Figure 13— Effect size and its standard error

Figure 14— Effect size and its inverse standard error

Egger’s test

Code

df_bias <- dt_all %>% mutate(sqrt_inv_e_n = sqrt((Cn + Tn)/(Cn * Tn)))

bias_model <- rma.mv(yi = lnRR,
                      V = lnRR_var, 
                      mods = ~ 1 + sqrt_inv_e_n,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(bias_model)

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

   logLik   Deviance        AIC        BIC       AICc   
-248.4547   496.9094   504.9094   519.2734   505.0615   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0914  0.3024     33     no  Study_ID 
sigma^2.2  0.2074  0.4554    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.0568, p-val = 0.8118

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt         0.2491  0.1353   1.8404  268  0.0668  -0.0174  0.5155  . 
sqrt_inv_e_n   -0.0895  0.3757  -0.2383  268  0.8118  -0.8291  0.6501    

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

Code

bubble_plot(bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Figure 15— Egger’s test of publication bias

Decline effect

Code

year_model <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~  1 + Year,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(year_model)

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

   logLik   Deviance        AIC        BIC       AICc   
-250.1678   500.3355   508.3355   522.6995   508.4876   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0510  0.2259     33     no  Study_ID 
sigma^2.2  0.2339  0.4836    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.0147, p-val = 0.9035

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt    1.7561  12.8324   0.1368  268  0.8913  -23.5091  27.0212    
Year      -0.0008   0.0064  -0.1214  268  0.9035   -0.0134   0.0118    

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

Code

bubble_plot(year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

Figure 16— Decline effect

Multi-moderator

We also ran Egger’s test and checked the decline effect in the multi-moderator model.

Code

multi_bias <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~ Treatment_stimulus + Log_area 
                            + Number_pattern + Type_prey
                            + sqrt_inv_e_n + Year,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(multi_bias)

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

   logLik   Deviance        AIC        BIC       AICc   
-240.1219   480.2438   498.2438   530.3931   498.9552   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0419  0.2048     33     no  Study_ID 
sigma^2.2  0.2265  0.4759    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:7):
F(df1 = 6, df2 = 263) = 2.3105, p-val = 0.0343

Model Results:

                               estimate       se     tval   df    pval 
intrcpt                         -7.1975  14.4257  -0.4989  263  0.6182 
Treatment_stimulusnon_eyespot   -0.0470   0.1067  -0.4400  263  0.6603 
Log_area                         0.1128   0.0481   2.3465  263  0.0197 
Number_pattern                  -0.0467   0.0289  -1.6155  263  0.1074 
Type_preyreal                    0.1894   0.1316   1.4390  263  0.1513 
sqrt_inv_e_n                    -0.4846   0.3906  -1.2407  263  0.2158 
Year                             0.0036   0.0071   0.5038  263  0.6148 
                                  ci.lb    ci.ub    
intrcpt                        -35.6020  21.2071    
Treatment_stimulusnon_eyespot   -0.2572   0.1632    
Log_area                         0.0181   0.2074  * 
Number_pattern                  -0.1037   0.0102    
Type_preyreal                   -0.0698   0.4487    
sqrt_inv_e_n                    -1.2536   0.2845    
Year                            -0.0104   0.0176    

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

Code

bubble_plot(multi_bias,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Code

bubble_plot(multi_bias,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

3.6 Additional analyses

Dataset

We compared the effect size between two datasets (predator and prey) to see whether there was any difference in effect size between the two datasets.

Code

mr_dataset <- rma.mv(yi = lnRR,
                      V = VCV,
                      mods = ~ Dataset,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

summary(mr_dataset)

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

   logLik   Deviance        AIC        BIC       AICc   
-249.9719   499.9439   507.9439   522.3078   508.0959   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0519  0.2279     33     no  Study_ID 
sigma^2.2  0.2331  0.4828    270     no    Obs_ID 

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

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.4221, p-val = 0.5164

Model Results:

             estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt        0.1563  0.0875  1.7852  268  0.0754  -0.0161  0.3286  . 
Datasetprey    0.0756  0.1164  0.6497  268  0.5164  -0.1535  0.3047    

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

Code

orchard_plot(mr_dataset,
            mod = "Dataset",
            group = "Study_ID",
            xlab = "Dataset", 
            angle = 45) 

Figure 17— Effect size and dataset - predator and prey

Total pattern area

We also ran several multi-moderator models.

Code

mr_all_2 <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Treatment_stimulus + Log_area_pattern_T 
                + Number_pattern + Type_prey,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)
summary(mr_all_2)

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

   logLik   Deviance        AIC        BIC       AICc   
-242.3093   484.6186   498.6186   523.6767   499.0544   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0540  0.2323     33     no  Study_ID 
sigma^2.2  0.2227  0.4720    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 265) = 3.0623, p-val = 0.0172

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
intrcpt                         -0.0912  0.2069  -0.4406  265  0.6599  -0.4986 
Treatment_stimulusnon_eyespot   -0.0291  0.1067  -0.2723  265  0.7856  -0.2391 
Log_area_pattern_T               0.0994  0.0429   2.3196  265  0.0211   0.0150 
Number_pattern                  -0.0756  0.0276  -2.7383  265  0.0066  -0.1300 
Type_preyreal                    0.1747  0.1364   1.2806  265  0.2015  -0.0939 
                                 ci.ub     
intrcpt                         0.3163     
Treatment_stimulusnon_eyespot   0.1810     
Log_area_pattern_T              0.1838   * 
Number_pattern                 -0.0212  ** 
Type_preyreal                   0.4432     

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

Moderators associated with conspicuousness

Multi-moderator model; only include moderators related to conspicuousness (pattern area/total pattern area and number of patterns) in the model.

3.6.0.1 pattern area and number of patterns

Code

mr_cons <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area + Number_pattern,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons)

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

   logLik   Deviance        AIC        BIC       AICc   
-244.9501   489.9002   499.9002   517.8365   500.1301   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0436  0.2087     33     no  Study_ID 
sigma^2.2  0.2262  0.4756    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 267) = 4.8661, p-val = 0.0084

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          -0.0185  0.1916  -0.0966  267  0.9231  -0.3959  0.3588    
Log_area          0.0899  0.0427   2.1056  267  0.0362   0.0058  0.1740  * 
Number_pattern   -0.0406  0.0267  -1.5215  267  0.1293  -0.0932  0.0119    

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

Code

r2_ml(mr_cons)

   R2_marginal R2_conditional 
    0.08098982     0.22938197 

Code

bubble_plot(mr_cons,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Log-transformed area")

3.6.0.2 Total pattern area and number of patterns

Code

mr_cons1 <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area_pattern_T + Number_pattern,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons1)

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

   logLik   Deviance        AIC        BIC       AICc   
-244.6226   489.2452   499.2452   517.1814   499.4751   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0453  0.2128     33     no  Study_ID 
sigma^2.2  0.2243  0.4736    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 267) = 5.2145, p-val = 0.0060

Model Results:

                    estimate      se     tval   df    pval    ci.lb    ci.ub    
intrcpt              -0.0483  0.1929  -0.2506  267  0.8024  -0.4282   0.3315    
Log_area_pattern_T    0.0934  0.0412   2.2653  267  0.0243   0.0122   0.1745  * 
Number_pattern       -0.0625  0.0254  -2.4612  267  0.0145  -0.1125  -0.0125  * 

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

Code

r2_ml(mr_cons1)

   R2_marginal R2_conditional 
    0.08546244     0.23909238 

Code

bubble_plot(mr_cons1,
            mod = "Log_area_pattern_T",
            group = "Study_ID",
            xlab = "Log-transformed area")

3.6.0.3 Pattern area, number of patterns, and treatment stimulus

Code

mr_cons2 <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area + Number_pattern + Treatment_stimulus,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons2)

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

   logLik   Deviance        AIC        BIC       AICc   
-244.1838   488.3676   500.3676   521.8686   500.6920   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0480  0.2190     33     no  Study_ID 
sigma^2.2  0.2262  0.4756    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 266) = 3.2831, p-val = 0.0214

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
intrcpt                         -0.0034  0.1994  -0.0170  266  0.9865  -0.3960 
Log_area                         0.0907  0.0437   2.0761  266  0.0388   0.0047 
Number_pattern                  -0.0405  0.0269  -1.5039  266  0.1338  -0.0935 
Treatment_stimulusnon_eyespot   -0.0506  0.1025  -0.4939  266  0.6218  -0.2524 
                                ci.ub    
intrcpt                        0.3893    
Log_area                       0.1767  * 
Number_pattern                 0.0125    
Treatment_stimulusnon_eyespot  0.1512    

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

Code

r2_ml(mr_cons2)

   R2_marginal R2_conditional 
    0.08567312     0.24568295 

3.6.0.4 Total pattern area, number of patterns, and treatment stimulus

Code

mr_cons3 <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area_pattern_T + Number_pattern + Treatment_stimulus,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons3)

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

   logLik   Deviance        AIC        BIC       AICc   
-243.8159   487.6317   499.6317   521.1327   499.9560   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0499  0.2234     33     no  Study_ID 
sigma^2.2  0.2243  0.4736    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 266) = 3.5441, p-val = 0.0151

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
intrcpt                         -0.0349  0.1999  -0.1744  266  0.8617  -0.4284 
Log_area_pattern_T               0.0950  0.0421   2.2546  266  0.0250   0.0120 
Number_pattern                  -0.0625  0.0256  -2.4438  266  0.0152  -0.1129 
Treatment_stimulusnon_eyespot   -0.0572  0.1029  -0.5564  266  0.5784  -0.2598 
                                 ci.ub    
intrcpt                         0.3587    
Log_area_pattern_T              0.1779  * 
Number_pattern                 -0.0122  * 
Treatment_stimulusnon_eyespot   0.1453    

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

Code

r2_ml(mr_cons3)

   R2_marginal R2_conditional 
    0.09147406     0.25682132 

Moderators associated with prey

Code

mr_pre <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Treatment_stimulus + Type_prey,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_pre)

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

   logLik   Deviance        AIC        BIC       AICc   
-249.3570   498.7140   508.7140   526.6503   508.9439   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0589  0.2426     33     no  Study_ID 
sigma^2.2  0.2333  0.4830    270     no    Obs_ID 

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

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 267) = 0.2046, p-val = 0.8151

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
intrcpt                          0.2056  0.0889   2.3125  267  0.0215   0.0306 
Treatment_stimulusnon_eyespot   -0.0569  0.1093  -0.5204  267  0.6032  -0.2721 
Type_preyreal                    0.0338  0.1308   0.2580  267  0.7966  -0.2238 
                                ci.ub    
intrcpt                        0.3807  * 
Treatment_stimulusnon_eyespot  0.1584    
Type_preyreal                  0.2913    

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

Code

r2_ml(mr_pre)

   R2_marginal R2_conditional 
   0.003607887    0.204348070 

4 More information

Bird tree 🦜

Code

trees <- read.nexus(here("data/bird_phy.nex"))
plot(trees[1])

Figure 18— Phylogenetic tree of bird species

Check phylogenetic relatedness 🌲

We tested whether phylogeny should be included in our analyses.

Code

# data 
dat_predator <- filter(dat_all, Dataset == "predator")

# check if the tree is ultrametric
# is.ultrametric(trees[[1]])
# [1] TRUE

# example of variance-covariance matrix
vcv(trees[[1]], corr = T)

                    Gallus_gallus Cyanocitta_cristata Sturnus_vulgaris
Gallus_gallus                   1           0.0000000        0.0000000
Cyanocitta_cristata             0           1.0000000        0.4437522
Sturnus_vulgaris                0           0.4437522        1.0000000
Ficedula_hypoleuca              0           0.4437522        0.6390457
Emberiza_sulphurata             0           0.4437522        0.5328950
Parus_major                     0           0.4437522        0.5288567
Parus_caeruleus                 0           0.4437522        0.5288567
                    Ficedula_hypoleuca Emberiza_sulphurata Parus_major
Gallus_gallus                0.0000000           0.0000000   0.0000000
Cyanocitta_cristata          0.4437522           0.4437522   0.4437522
Sturnus_vulgaris             0.6390457           0.5328950   0.5288567
Ficedula_hypoleuca           1.0000000           0.5328950   0.5288567
Emberiza_sulphurata          0.5328950           1.0000000   0.5288567
Parus_major                  0.5288567           0.5288567   1.0000000
Parus_caeruleus              0.5288567           0.5288567   0.8176725
                    Parus_caeruleus
Gallus_gallus             0.0000000
Cyanocitta_cristata       0.4437522
Sturnus_vulgaris          0.5288567
Ficedula_hypoleuca        0.5288567
Emberiza_sulphurata       0.5288567
Parus_major               0.8176725
Parus_caeruleus           1.0000000

Code

# calculate lnRR and lnRR_var
dat_pred <- effect_lnRR(dat_predator)
dat_pred$Obs_ID <- 1:nrow(dat_pred)
V <- vcalc(lnRR_var, cluster = Cohort_ID, subgroup = Obs_ID, data = dat_pred)

# create function to run meta-analysis for 50 trees
phy_model <- function(cor_tree = vcv_tree){
            model <- rma.mv(yi = lnRR,
                            V = V,
                            random = list(~1 | Study_ID,
                                          ~1 | Shared_control_ID,
                                          ~1 | Cohort_ID,
                                          ~1 | Obs_ID,
                                          ~1 | Bird_species,
                                          ~1 | Phylogeny),
                            R = list(Phylogeny = cor_tree),
                            test = "t",
                            method = "REML",
                            sparse = TRUE,
                            data = dat_pred)
  model
}

# running 50 meta-analyses with 50 different trees
# We got 1000 phylogenetic trees from BirdTree.org (https://birdtree.org)
# Only 50 trees are used as in Nakagawa & Villemereuil (2019)(https://doi.org/10.1093/sysbio/syy089)

tree_50 <- trees[1:50]
vcv_tree_50 <- map(tree_50, ~vcv(.x, corr = TRUE))
ma_50 <- mclapply(vcv_tree_50, phy_model) # detectCores() = 12

# always best to save RDS files for analysis 
# which takes more than a min or so 
# saveRDS(ma_50, here("Rdata", "ma_50.RDS")) 

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

# combining the results
est_50 <- map_dbl(ma_50, ~ .x$b[[1]])
se_50 <-  map_dbl(ma_50, ~ .x$se)
df_50 <- c(rbind(est_50, se_50))

# creating an array required by pool.mi
my_array <- array(df_50, dim = c(1, 2, 50))

pool.mi(my_array)

#> Pooling estimates from 50 imputed analyses for 1 parameters. 
  pooled.mean pooled.total.se frac.miss.info   CI.1  CI.2 P.value
1       0.139           0.119              0 -0.094 0.373   0.242

  pooled.mean pooled.total.se se.within  se.between relative.inc.var
1   0.1393776       0.1192319 0.1192319 2.54026e-08     4.629898e-14
  frac.miss.info        CI.1      CI.2   p.value
1    2.00061e-06 -0.09431294 0.3730681 0.2424192

Code

# extract sigma^2 for averaging variance components
sigma2_mod <- do.call(rbind, lapply(ma_50, function(x) x$sigma2))
sigma2_mod <- data.frame(sigma2_mod)

colnames(sigma2_mod) <- c("sigma^2.1_Study_ID", "sigma^2.2_SharedControl_ID", 
                          "sigma^2.3_Cohort_ID", "sigma^2.4_Obs_ID", 
                          "sigma^2.5_BirdSpecies", "sigma^2.6_Phylogeny")

# easier to undersatnd if you round
round(colMeans(sigma2_mod), 2)

        sigma^2.1_Study_ID sigma^2.2_SharedControl_ID 
                      0.00                       0.09 
       sigma^2.3_Cohort_ID           sigma^2.4_Obs_ID 
                      0.10                       0.53 
     sigma^2.5_BirdSpecies        sigma^2.6_Phylogeny 
                      0.00                       0.00 

Data of Butterfly eyespots 🦋

Code

butterfly <- read.csv(here("data/butterfly_eyespots.csv"), header = T, fileEncoding = "CP932")
datatable(butterfly, extensions = "Buttons")

5 Figure gallery

These were the basis for the final figures presented in the paper.

Discrete moderators

Code

# extract the results of each model and combine them into one data frame for plotting
main  <- mod_results(ma_all, group = "Study_ID")
pattern <- mod_results(mr_eyespot, mod = "Treatment_stimulus", group = "Study_ID")
shape  <- mod_results(mr_prey_type, mod = "Type_prey", group = "Study_ID")

comb <- submerge(main, pattern, shape)

p1 <- orchard_plot(comb,
                    group = "Study_ID",
                    xlab = "log response ratio (lnRR)",
                    tree.order = c( "Artificial", "Real","Non_eyespot", "Eyespot", "Intrcpt"),
                    angle = 45,
                    branch.size = 2) +
                    theme(legend.direction = "horizontal",
                    legend.text = element_text(size = 6)) +
                    scale_colour_brewer(palette = "Pastel1") +
                    scale_fill_brewer(palette = "Pastel1") 
p1

p2 <- orchard_plot(mr_prey_shape,
                  mod = "Shape_prey",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6)) +
                  scale_colour_brewer(palette = "YlGnBu") +
                  scale_fill_brewer(palette = "YlGnBu")
p2

Figure 19— main figure - discrete moderators

Figure 20— Prey type

Continuous moderators

Code

# these figs were created by multi meta-regression model results using the log-transformed area

p3 <- bubble_plot(mr_area,
                  mod = "Log_area",
                  group = "Study_ID",
                  xlab = "Log-transformed area",
                  est.lwd = 1.2,
                  est.col = "deeppink3") +
                  scale_x_continuous(breaks = seq(0,7,1.5)) +
                  scale_y_continuous(breaks = seq(-3.0, 4.5, 1.0))

# number of patterns
p4 <- bubble_plot(mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns",
            est.lwd = 1.2,
            est.col = "turquoise4") +
      scale_x_continuous(breaks = seq(0,11,2)) +
      scale_y_continuous(breaks = seq(-3.0, 4.5, 1.0))

# combine
p3 / p4 + plot_annotation(tag_levels = "a")

# total pattern area
p5 <- bubble_plot(mr_area_T,
            mod = "Log_area_pattern_T",
            group = "Study_ID",
            xlab = "Log-transformed total pattern area")

# diameter
p6 <- bubble_plot(mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

# background
p7 <- bubble_plot(mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed total prey surface area")

p5 / p6 /p7 + plot_annotation(tag_levels = "a")

Figure 21— main figure - continuous moderators

Figure 22— pattern area

Publication bias

Code

# funnel plot

#pdf("funnel.pdf")

funnel(ma_all, yaxis = "seinv", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",
      # xlim = c(-4.0, 4.5), 
      # ylim = c(0.01, 60.0), 
      col = c(alpha("black", 0.5)),
      cex = 1.5) 

# pdf("funnel2.pdf", width = 15, height = 10)
# dev.off()

# Egger's test and decline effect
p8 <- bubble_plot(bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size") +
            scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5)) 

p9 <- bubble_plot(year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication") +  scale_y_continuous(breaks = seq(-2.5, 4.0, 1.5)) 

# combine plots created by bubble_plot()
pub <- p8 / p9
pub +  plot_annotation(tag_levels = 'a') 

Figure 23— Funnel plot

Figure 24— Egger’s test and Decline effect

Publication bias (multi-moderator)

Code

p10 <- bubble_plot(multi_bias,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

p11 <- bubble_plot(multi_bias,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

p10 / p11 +  plot_annotation(tag_levels = 'a')

Figure 25— Egger’s test and Decline effect (multi-moderator)

6 R session information

R version 4.3.1 (2023-06-16)

Platform: aarch64-apple-darwin20 (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: parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: RColorBrewer(v.1.1-3), orchaRd(v.2.0), multcomp(v.1.4-25), TH.data(v.1.1-2), MASS(v.7.3-60), survival(v.3.5-7), mvtnorm(v.1.2-3), miWQS(v.0.4.4), metafor(v.4.4-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.6-1.1), patchwork(v.1.1.3), 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), tidyverse(v.2.0.0), knitr(v.1.45), here(v.1.0.1), ggcorrplot(v.0.1.4.1), ggplot2(v.3.4.4), ggokabeito(v.0.1.0), dtplyr(v.1.3.1), DT(v.0.30), GoodmanKruskal(v.0.0.3) and ape(v.5.7-1)

loaded via a namespace (and not attached): gridExtra(v.2.3), glm2(v.1.2.1), sandwich(v.3.0-2), rlang(v.1.1.2), magrittr(v.2.0.3), compiler(v.4.3.1), mgcv(v.1.9-0), reshape2(v.1.4.4), vctrs(v.0.6.4), quantreg(v.5.97), pkgconfig(v.2.0.3), fastmap(v.1.1.1), ellipsis(v.0.3.2), backports(v.1.4.1), mcmc(v.0.9-7), labeling(v.0.4.3), pander(v.0.6.5), utf8(v.1.2.4), rmarkdown(v.2.25), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), MatrixModels(v.0.5-2), xfun(v.0.41), cachem(v.1.0.8), jsonlite(v.1.8.7), gmm(v.1.8), cluster(v.2.1.4), R6(v.2.5.1), bslib(v.0.6.0), stringi(v.1.8.1), Rsolnp(v.1.16), rlist(v.0.4.6.2), rpart(v.4.1.21), estimability(v.1.4.1), jquerylib(v.0.1.4), Rcpp(v.1.0.11), zoo(v.1.8-12), tmvmixnorm(v.1.1.1), base64enc(v.0.1-3), pacman(v.0.5.1), splines(v.4.3.1), nnet(v.7.3-19), timechange(v.0.2.0), tidyselect(v.1.2.0), rstudioapi(v.0.15.0), yaml(v.2.3.7), codetools(v.0.2-19), plyr(v.1.8.9), lattice(v.0.21-9), withr(v.2.5.2), coda(v.0.19-4), evaluate(v.0.23), tmvtnorm(v.1.5), foreign(v.0.8-85), pillar(v.1.9.0), corrplot(v.0.92), checkmate(v.2.3.0), stats4(v.4.3.1), generics(v.0.1.3), rprojroot(v.2.0.3), invgamma(v.1.1), mathjaxr(v.1.6-0), truncnorm(v.1.0-9), hms(v.1.1.3), munsell(v.0.5.0), scales(v.1.2.1), xtable(v.1.8-4), glue(v.1.6.2), emmeans(v.1.8.9), Hmisc(v.5.1-1), tools(v.4.3.1), data.table(v.1.14.8), SparseM(v.1.81), matrixNormal(v.0.1.1), grid(v.4.3.1), MCMCpack(v.1.6-3), crosstalk(v.1.2.0), colorspace(v.2.1-0), nlme(v.3.1-163), beeswarm(v.0.4.0), condMVNorm(v.2020.1), vipor(v.0.4.5), htmlTable(v.2.4.1), latex2exp(v.0.9.6), Formula(v.1.2-5), cli(v.3.6.1), fansi(v.1.0.5), gtable(v.0.3.4), sass(v.0.4.7), digest(v.0.6.33), farver(v.2.1.1), htmlwidgets(v.1.6.3), htmltools(v.0.5.7) and lifecycle(v.1.0.4)