Pre-copulatory sexual selection predicts sexual size dimorphism: a meta-analysis of comparative studies
Main analyses
Lennart Winkler1, Robert P
Freckleton2, Tamas Szekely3, 4 &
Tim Janicke1,5
1Applied Zoology,
Technical University Dresden 2Department of
Zoology, University of Oxford, South Parks Road, Oxford OX1 3PS,
UK
3Milner Centre for Evolution, University of Bath,
Bath, UK84Department of Evolutionary Zoology and
Human Behaviour, University of Debrecen, Debrecen, Hungary
5Centre d’Écologie Fonctionnelle et Évolutive, UMR 5175,
CNRS, Université de Montpellier
Last updated: 2024-08-16
Checks: 7 0
Knit directory:
SSD_and_sexual_selection_2024/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20230430) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version cd21614. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rhistory
Unstaged changes:
Modified: analysis/_site.yml
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/index.Rmd) and HTML
(docs/index.html) files. If you’ve configured a remote Git
repository (see ?wflow_git_remote), click on the hyperlinks
in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | cd21614 | LennartWinkler | 2024-08-16 | Build site. |
| Rmd | 6eb630e | LennartWinkler | 2024-08-16 | wflow_publish(republish = TRUE, all = T) |
| html | 884cafd | LennartWinkler | 2024-05-11 | Build site. |
| Rmd | 84c3659 | LennartWinkler | 2024-05-11 | wflow_publish(all = T) |
| html | 84c3659 | LennartWinkler | 2024-05-11 | wflow_publish(all = T) |
| html | 66fef12 | LennartWinkler | 2024-05-01 | Build site. |
| Rmd | 7fd037e | LennartWinkler | 2024-05-01 | wflow_publish(all = T) |
| html | 7fd037e | LennartWinkler | 2024-05-01 | wflow_publish(all = T) |
| html | b028b7d | LennartWinkler | 2023-05-24 | Build site. |
| Rmd | c1d3311 | LennartWinkler | 2023-05-24 | wflow_publish(all = T) |
| html | 71d667c | LennartWinkler | 2023-05-03 | Build site. |
| Rmd | 378a58c | LennartWinkler | 2023-05-03 | wflow_publish(republish = TRUE, all = T) |
| html | 085f45f | LennartWinkler | 2023-05-03 | Build site. |
| Rmd | 85c143e | LennartWinkler | 2023-05-03 | update |
| html | 8edd942 | LennartWinkler | 2023-05-01 | Build site. |
| Rmd | d131cd9 | LennartWinkler | 2023-05-01 | wflow_publish(republish = TRUE, all = T) |
| Rmd | 57ca562 | LennartWinkler | 2023-04-30 | update |
| html | 054a12b | LennartWinkler | 2023-04-30 | Build site. |
| Rmd | 1a5561c | LennartWinkler | 2023-04-30 | wflow_publish(all = T) |
| Rmd | 4482c88 | LennartWinkler | 2023-04-30 | Start workflowr project. |
Supplementary material reporting R code for the manuscript ‘Pre-copulatory sexual selection predicts sexual size dimorphism: a meta-analysis of comparative studies’. Additional analyses excluding studies that did not control for phylogenetic non-independence (Supplement 2) can be found at: https://https://lennartwinkler.github.io/SSD_and_sexual_selection_2024/index2.html # Load and prepare data Before we started the analyses, we loaded all necessary packages and data.
rm(list = ls()) # Clear work environment
# Load R-packages ####
list_of_packages=cbind('ape','matrixcalc','metafor','Matrix','MASS','pwr','psych','multcomp','data.table','ggplot2','RColorBrewer','MCMCglmm','ggdist','cowplot','PupillometryR','dplyr','wesanderson','gridExtra')
lapply(list_of_packages, require, character.only = TRUE)
# Load data set ####
MetaData <- read.csv("./data/Data_SexSelSSD.csv", sep=";", header=TRUE) # Load data set
N_Studies <- length(summary(as.factor(MetaData$Study_ID))) # Number of included primary studies
Tree<- read.tree("./data/Pyologeny_SexSelSSD.txt") # Load phylogenetic tree
# Prune phylogenetic tree
MetaData_Class_Data <- unique(MetaData$Class)
Tree_Class<-drop.tip(Tree, Tree$tip.label[-na.omit(match(MetaData_Class_Data, Tree$tip.label))])
forcedC_Moderators <- as.matrix(forceSymmetric(vcv(Tree_Class, corr=TRUE)))
# Order moderator levels
MetaData$SexSel_Episode=as.factor(MetaData$SexSel_Episode)
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("post-copulatory"))
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("pre-copulatory"))
MetaData$Class=as.factor(MetaData$Class)
MetaData$z=as.numeric(MetaData$z)
# Set figure theme and colors
theme=theme(panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(0.2,0.5),
legend.title = element_blank(),
legend.text = element_text(colour="black", size=12),
axis.line.x = element_line(colour = "black", linewidth = 1),
axis.line.y = element_line(colour = "black", linewidth = 1),
axis.text.x = element_text(face="plain", color="black", size=12, angle=0),
axis.text.y = element_text(face="plain", color="black", size=12, angle=0),
axis.title.x = element_text(size=14,face="plain", margin = margin(r=0,10,0,0)),
axis.title.y = element_text(size=14,face="plain", margin = margin(r=10,0,0,0)),
axis.ticks = element_line(linewidth = 1),
axis.ticks.length = unit(.3, "cm"))
colpal=c('#003d6cff','#0c7accff','#5a9c97ff')
colpal2=brewer.pal(7, 'Dark2')
colpal3=brewer.pal(10, 'BrBG')
colpal4=c("grey50","grey65")
colpal5=c("#8C2D04","#CC4C02", "#EC7014", "#FE9929", "#FEC44F","#A8DDB5","#A8DDB5", "#7BCCC4", "#4EB3D3", "#2B8CBE", "#08589E")
colpal6=c("#8C2D04","#CC4C02", "#EC7014", "#FE9929", "#FEC44F","#A8DDB5", "#7BCCC4", "#4EB3D3", "#2B8CBE", "#08589E")
Meta_col=c('#5a9c97ff','#0c7accff','#003d6cff','black')Global models
We addressed the question if increasing sexual selection correlated with an increasingly male-biased SSD. For this we ran a global model including an observation-level index and the study identifier as random termson correlation coefficients that were positive if increasing sexual selection correlated with an increasingly male-biased SSD, but negative if increasing sexual selection correlated with an increasingly female-biased SSD. First, we ran a global model including the phylogeny:
Model_REML_Null = rma.mv(r ~ 1, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_Null)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-24.1841 48.3683 56.3683 67.7134 56.6989
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0215 0.1467 77 no Study_ID no
sigma^2.2 0.0469 0.2166 127 no Index no
sigma^2.3 0.0118 0.1084 11 no Class yes
Test for Heterogeneity:
Q(df = 126) = 1607.5311, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2584 0.0708 3.6509 0.0003 0.1197 0.3971 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Second, we ran a global model without the phylogeny:
Model_cREML_Null = rma.mv(r ~ 1, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_Null)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-25.6599 51.3199 57.3199 65.8287 57.5166
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0317 0.1781 77 no Study_ID
sigma^2.2 0.0451 0.2123 127 no Index
Test for Heterogeneity:
Q(df = 126) = 1607.5311, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2767 0.0312 8.8726 <.0001 0.2156 0.3378 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Moderator tests for phylogenetic models
Next, we ran a series of models that test the effect of different moderators. Again we started with models including the phylogeny.
Sexual selection episode
The first model explores the effect of the different episodes of sexual selection (i.e. pre-copulatory, post-copulatory or both):
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("pre-copulatory"))
Model_REML_by_SexSelEpisode = rma.mv(r ~ factor(SexSel_Episode), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelEpisode)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-7.9463 15.8925 27.8925 44.8142 28.6105
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0171 0.1308 77 no Study_ID no
sigma^2.2 0.0326 0.1807 127 no Index no
sigma^2.3 0.0115 0.1073 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 124) = 1378.0070, p-val < .0001
Test of Moderators (coefficients 2:3):
QM(df = 2) = 38.0922, p-val < .0001
Model Results:
estimate se zval pval
intrcpt 0.2582 0.0731 3.5332 0.0004
factor(SexSel_Episode)post-copulatory -0.3591 0.0805 -4.4613 <.0001
factor(SexSel_Episode)both 0.1412 0.0519 2.7198 0.0065
ci.lb ci.ub
intrcpt 0.1150 0.4014 ***
factor(SexSel_Episode)post-copulatory -0.5168 -0.2013 ***
factor(SexSel_Episode)both 0.0394 0.2429 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("post-copulatory"))
Model_REML_by_SexSelEpisode2 = rma.mv(r ~ factor(SexSel_Episode), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelEpisode2)
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("both"))
Model_REML_by_SexSelEpisode3 = rma.mv(r ~ factor(SexSel_Episode), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelEpisode3)Finally, we computed FDR corrected p-values:
tab1=as.data.frame(round(p.adjust(c(0.0005, 0.2849, .0001), method = 'fdr'),digit=3),row.names=cbind("Pre-copulatory","Post-copulatory","Both"))
colnames(tab1)<-cbind('P-value')
tab1 P-value
Pre-copulatory 0.001
Post-copulatory 0.285
Both 0.000Plot sexual selection mode (Figure 3)
Here we plot the sexual selection mode moderator:
MetaData$SexSel_Episode=factor(MetaData$SexSel_Episode, levels = c("both","post-copulatory" ,"pre-copulatory"))
ggplot(MetaData, aes(x=SexSel_Episode, y=r, fill = SexSel_Episode, colour = SexSel_Episode)) +
geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=1)+
geom_flat_violin(position = position_nudge(x = 0.25, y = 0),adjust =1, trim = F,alpha=0.6)+
geom_point(position = position_jitter(width = .1), size = 2.5,alpha=0.6,stroke=0,shape=19)+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelEpisode$ci.lb[1], x=3.25, xend= 3.25, yend= Model_REML_by_SexSelEpisode$ci.ub[1]), alpha=1,linewidth=1,color='black')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelEpisode2$ci.lb[1], x=2.25, xend= 2.25, yend= Model_REML_by_SexSelEpisode2$ci.ub[1]), alpha=1,linewidth=1,color='black')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelEpisode3$ci.lb[1], x=1.25, xend= 1.25, yend= Model_REML_by_SexSelEpisode3$ci.ub[1]), alpha=1,linewidth=1,color='black')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelEpisode$b[1,1], x=3.25), size = 3.5,alpha=1,stroke=1,shape=21,fill='grey50',color='black')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelEpisode2$b[1,1], x=2.25), size = 3.5,alpha=1,stroke=1,shape=21,fill='white',color='black')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelEpisode3$b[1,1], x=1.25), size = 3.5,alpha=1,stroke=1,shape=21,fill='black',color='black')+
ylab(expression(paste("Effect size (", italic("r"),')')))+xlab('Sexual selection episode')+coord_flip()+guides(fill = FALSE, colour = FALSE) +
scale_color_manual(values =colpal)+
scale_fill_manual(values =colpal)+
scale_x_discrete(labels=c("Pre- and \n post-copulatory","Post-copulatory" ,"Pre-copulatory"),expand=c(.1,0))+
annotate("text", x=1, y=1.2, label= "n = 52",size=4.5) +
annotate("text", x=2, y=1.2, label= "n = 19",size=4.5) +
annotate("text", x=3, y=1.2, label= "n = 56",size=4.5) + theme
Figure 3: Raincloud plot of correlation coefficients between SSD and the modes of sexual selection proxies (i.e. pre-copulatory, post-copulatory or both) including sample sizes and estimates with 95%CI from phylogenetic model.
Sexual selection category
Next we explored the effect of the sexual selection category (i.e. density, mating system, operational sex ratio (OSR), post-mating competition, pre-mating competition, trait-based, other):
MetaData$SexSel_Category=as.factor(MetaData$SexSel_Category)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Postmating competition"))
Model_REML_by_SexSelCat = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelCat)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
6.9109 -13.8217 6.1783 34.0532 8.1966
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0206 0.1436 77 no Study_ID no
sigma^2.2 0.0180 0.1343 127 no Index no
sigma^2.3 0.0103 0.1014 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 120) = 949.5416, p-val < .0001
Test of Moderators (coefficients 2:7):
QM(df = 6) = 86.9845, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.1169 0.0881 -1.3261 0.1848
factor(SexSel_Category)Density 0.3307 0.1204 2.7468 0.0060
factor(SexSel_Category)Mating system 0.4846 0.0786 6.1635 <.0001
factor(SexSel_Category)OSR 0.7187 0.1183 6.0759 <.0001
factor(SexSel_Category)Other 0.4501 0.0926 4.8600 <.0001
factor(SexSel_Category)Premating competition 0.4668 0.0805 5.7991 <.0001
factor(SexSel_Category)Trait-based 0.1204 0.0873 1.3800 0.1676
ci.lb ci.ub
intrcpt -0.2896 0.0559
factor(SexSel_Category)Density 0.0947 0.5667 **
factor(SexSel_Category)Mating system 0.3305 0.6387 ***
factor(SexSel_Category)OSR 0.4868 0.9505 ***
factor(SexSel_Category)Other 0.2686 0.6316 ***
factor(SexSel_Category)Premating competition 0.3091 0.6246 ***
factor(SexSel_Category)Trait-based -0.0506 0.2914
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Trait-based"))
Model_REML_by_SexSelCat2 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelCat2)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Density"))
Model_REML_by_SexSelCat3 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelCat3)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Premating competition"))
Model_REML_by_SexSelCat4 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelCat4)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Mating system"))
Model_REML_by_SexSelCat5 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelCat5)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("OSR"))
Model_REML_by_SexSelCat6 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelCat6)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Other"))
Model_REML_by_SexSelCat7 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SexSelCat7)Finally, we computed FDR corrected p-values:
tab2=as.data.frame(round(p.adjust(c(0.1888, 0.9711, 0.0722, .0001, .0001, .0001, .0001), method = 'fdr'),digit=3),row.names=cbind("Postmating competition","Trait-based","Density",'Premating competition',"Mating system","OSR","Other"))
colnames(tab2)<-cbind('P-value')
tab2 P-value
Postmating competition 0.220
Trait-based 0.971
Density 0.101
Premating competition 0.000
Mating system 0.000
OSR 0.000
Other 0.000Plot sexual selection category (Figure S1)
Here we plot the sexual selection category moderator:
MetaData$SexSel_Category=factor(MetaData$SexSel_Category, levels = rev(c( "Postmating competition","Trait-based" ,"Density","Premating competition" ,"Mating system" , "OSR", "Other")))
ggplot(MetaData, aes(x=SexSel_Category, y=r, fill = SexSel_Category, colour = SexSel_Category)) +
geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=1)+
geom_flat_violin(position = position_nudge(x = 0.25, y = 0),adjust =1, trim = F,alpha=0.6)+
geom_point(position = position_jitter(width = .1), size = 2.5,alpha=0.6,stroke=0,shape=19)+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat$b[1,1], x=7.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat2$b[1,1], x=6.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat3$b[1,1], x=5.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat4$b[1,1], x=4.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat5$b[1,1], x=3.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat6$b[1,1], x=2.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat7$b[1,1], x=1.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat$ci.lb[1], x=7.25, xend= 7.25, yend= Model_REML_by_SexSelCat$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat2$ci.lb[1], x=6.25, xend= 6.25, yend= Model_REML_by_SexSelCat2$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat3$ci.lb[1], x=5.25, xend= 5.25, yend= Model_REML_by_SexSelCat3$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat4$ci.lb[1], x=4.25, xend= 4.25, yend= Model_REML_by_SexSelCat4$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat5$ci.lb[1], x=3.25, xend= 3.25, yend= Model_REML_by_SexSelCat5$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat6$ci.lb[1], x=2.25, xend= 2.25, yend= Model_REML_by_SexSelCat6$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SexSelCat7$ci.lb[1], x=1.25, xend= 1.25, yend= Model_REML_by_SexSelCat7$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
ylab(expression(paste("Effect size (", italic("r"),')')))+xlab('Sexual selection category')+coord_flip()+guides(fill = FALSE, colour = FALSE) +
scale_color_manual(values =colpal2)+
scale_fill_manual(values =colpal2)+
scale_x_discrete(labels=rev(c( "Post-mating\ncompetition","Trait-based" ,"Population\ndensity","Pre-mating\ncompetition" ,"Mating\nsystem" , "Sex ratio", "Other")),expand=c(.1,0))+
annotate("text", x=7, y=1.2, label= "n = 19",size=4.5) +
annotate("text", x=6, y=1.2, label= "n = 18",size=4.5) +
annotate("text", x=5, y=1.2, label= "n = 5",size=4.5) +
annotate("text", x=4, y=1.2, label= "n = 25",size=4.5) +
annotate("text", x=3, y=1.2, label= "n = 36",size=4.5) +
annotate("text", x=2, y=1.2, label= "n = 7",size=4.5) +
annotate("text", x=1, y=1.2, label= "n = 17",size=4.5) + theme
Figure S1: Raincloud plot of correlation coefficients between SSD and different sexual selection categories including sample sizes and estimates with 95%CI from phylogenetic model.
Type of SSD measure
Next we explored the effect of the type of SSD measure (i.e. body mass or size):
MetaData$SSD_Proxy=as.factor(MetaData$SSD_Proxy)
MetaData$SSD_Proxy=relevel(MetaData$SSD_Proxy,c("Body mass"))
Model_REML_by_SSDMeasure = rma.mv(r ~ SSD_Proxy, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SSDMeasure)
Multivariate Meta-Analysis Model (k = 125; method: REML)
logLik Deviance AIC BIC AICc
-23.9352 47.8704 57.8704 71.9314 58.3833
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0244 0.1561 76 no Study_ID no
sigma^2.2 0.0470 0.2167 125 no Index no
sigma^2.3 0.0039 0.0625 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 123) = 1553.9370, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 2.7123, p-val = 0.0996
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3409 0.0678 5.0296 <.0001 0.2081 0.4737 ***
SSD_ProxyBody size -0.1055 0.0641 -1.6469 0.0996 -0.2311 0.0201 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$SSD_Proxy=relevel(MetaData$SSD_Proxy,c("Body size"))
Model_REML_by_SSDMeasure2 = rma.mv(r ~ SSD_Proxy, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_SSDMeasure2)Finally, we computed FDR corrected p-values:
tab3=as.data.frame(round(p.adjust(c(.0001, .0001), method = 'fdr'),digit=3),row.names=cbind("Body mass","Body size"))
colnames(tab3)<-cbind('P-value')
tab3 P-value
Body mass 0
Body size 0Plot Type of SSD measure (Figure S2A)
Here we plot the type of SSD measure moderator:
MetaData_NAProxy=MetaData[!is.na(MetaData$SSD_Proxy),]
MetaData_NAProxy$SSD_Proxy=as.factor(MetaData_NAProxy$SSD_Proxy)
MetaData_NAProxy$SSD_Proxy=relevel(MetaData_NAProxy$SSD_Proxy,c("Body size"))
ggplot(MetaData_NAProxy, aes(x=SSD_Proxy, y=r, fill = SSD_Proxy, colour = SSD_Proxy)) +
geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=1)+
geom_flat_violin(position = position_nudge(x = 0.25, y = 0),adjust =1, trim = F,alpha=0.6)+
geom_point(position = position_jitter(width = .1), size = 2.5,alpha=0.6,stroke=0,shape=19)+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SSDMeasure$b[1,1], x=2.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_SSDMeasure2$b[1,1], x=1.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SSDMeasure$ci.lb[1], x=2.25, xend= 2.25, yend= Model_REML_by_SSDMeasure$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_SSDMeasure2$ci.lb[1], x=1.25, xend= 1.25, yend= Model_REML_by_SSDMeasure2$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
ylab(expression(paste("Effect size (", italic("r"),')')))+xlab('SSD proxy')+coord_flip()+guides(fill = FALSE, colour = FALSE) +
scale_color_manual(values =colpal4)+
scale_fill_manual(values =colpal4)+
scale_x_discrete(labels=rev(c("Body \nmass ","Body \n size ")),expand=c(.15,0))+
annotate("text", x=2, y=1.1, label= "n = 58",size=3.5) +
annotate("text", x=1, y=1.1, label= "n = 67",size=3.5) + theme
Figure S2A: Raincloud plot of correlation coefficients for different types of SSD measures (i.e. body mass or size) including sample sizes and estimates with 95% CI from phylogenetic model.
SSD measure controlled for body size?
Next we explored the effect if the primary study controlled the SSD for body size (i.e. uncontrolled or controlled):
MetaData$BodySizeControlled=as.factor(MetaData$BodySizeControlled)
MetaData$BodySizeControlled=relevel(MetaData$BodySizeControlled,c("No"))
Model_REML_by_BodySizeCont = rma.mv(r ~ BodySizeControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_BodySizeCont)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-23.9170 47.8340 57.8340 71.9756 58.3382
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0215 0.1467 77 no Study_ID no
sigma^2.2 0.0470 0.2169 127 no Index no
sigma^2.3 0.0110 0.1050 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 125) = 1601.6900, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 0.9348, p-val = 0.3336
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2732 0.0709 3.8520 0.0001 0.1342 0.4121 ***
BodySizeControlledYes -0.0825 0.0853 -0.9668 0.3336 -0.2496 0.0847
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$BodySizeControlled=relevel(MetaData$BodySizeControlled,c("Yes"))
Model_REML_by_BodySizeCont2 = rma.mv(r ~ BodySizeControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_BodySizeCont2)Finally, we computed FDR corrected p-values:
tab4=as.data.frame(round(p.adjust(c(0.0002, 0.0566), method = 'fdr'),digit=3),row.names=cbind("uncontrolled","controlled"))
colnames(tab4)<-cbind('P-value')
tab4 P-value
uncontrolled 0.000
controlled 0.057Plot: SSD measure controlled for body size? (Figure S2B)
Here we plot effect sizes if type of SSD measure controlled for body size:
MetaData$BodySizeControlled=as.factor(MetaData$BodySizeControlled)
MetaData$BodySizeControlled=relevel(MetaData$BodySizeControlled,c("Yes"))
ggplot(MetaData, aes(x=BodySizeControlled, y=r, fill = BodySizeControlled, colour = BodySizeControlled)) +
geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=1)+
geom_flat_violin(position = position_nudge(x = 0.25, y = 0),adjust =1, trim = F,alpha=0.6)+
geom_point(position = position_jitter(width = .1), size = 2.5,alpha=0.6,stroke=0,shape=19)+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_BodySizeCont$b[1,1], x=2.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_REML_by_BodySizeCont2$b[1,1], x=1.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_BodySizeCont$ci.lb[1], x=2.25, xend= 2.25, yend= Model_REML_by_BodySizeCont$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_REML_by_BodySizeCont2$ci.lb[1], x=1.25, xend= 1.25, yend= Model_REML_by_BodySizeCont2$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
ylab(expression(paste("Effect size (", italic("r"),')')))+xlab('Controlling SSD \n for body size or mass')+coord_flip()+guides(fill = FALSE, colour = FALSE) +
scale_color_manual(values =colpal4)+
scale_fill_manual(values =colpal4)+
scale_x_discrete(labels=(c("Controlled","Uncontrolled")),expand=c(.15,0))+
annotate("text", x=2, y=1.1, label= "n = 110",size=3.5) +
annotate("text", x=1, y=1.1, label= "n = 17",size=3.5) + theme
Figure S2B: Raincloud plot of correlation coefficients for primary studies controlling SSD for body size or mass (uncontrolled or controlled) including sample sizes and estimates with 95% CI from phylogenetic model.
Phylogentic vs non-phylogenetic studies
Next we explored the effect of studies controlling for phylogeny vs non-phylogenetic studies:
MetaData$PhyloControlled=as.factor(MetaData$PhyloControlled)
MetaData$PhyloControlled=relevel(MetaData$PhyloControlled,c("No"))
Model_cREML_by_Phylo = rma.mv(r ~ PhyloControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_cREML_by_Phylo)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.7974 41.5949 51.5949 65.7364 52.0991
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0133 0.1154 77 no Study_ID no
sigma^2.2 0.0491 0.2216 127 no Index no
sigma^2.3 0.0120 0.1097 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 125) = 1428.7953, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 7.7357, p-val = 0.0054
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.4226 0.0926 4.5634 <.0001 0.2411 0.6040 ***
PhyloControlledYes -0.2052 0.0738 -2.7813 0.0054 -0.3498 -0.0606 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$PhyloControlled=relevel(MetaData$PhyloControlled,c("Yes"))
Model_cREML_by_Phylo2 = rma.mv(r ~ PhyloControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_cREML_by_Phylo2)Finally, we computed FDR corrected p-values:
tab5=as.data.frame(round(p.adjust(c(0.0001, 0.0028), method = 'fdr'),digit=3),row.names=cbind("Non-phylogenetic","With phylogeny"))
colnames(tab5)<-cbind('P-value')
tab5 P-value
Non-phylogenetic 0.000
With phylogeny 0.003Plot Phylogentic vs non-phylogenetic studies (Figure S2C)
Here we plot the type of SSD measure moderator:
MetaData$PhyloControlled=as.factor(MetaData$PhyloControlled)
MetaData$PhyloControlled=relevel(MetaData$PhyloControlled,c("No"))
ggplot(MetaData, aes(x=PhyloControlled, y=r, fill = PhyloControlled, colour = PhyloControlled)) +
geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=1)+
geom_flat_violin(position = position_nudge(x = 0.25, y = 0),adjust =1, trim = F,alpha=0.6)+
geom_point(position = position_jitter(width = .1), size = 2.5,alpha=0.6,stroke=0,shape=19)+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Phylo2$b[1,1], x=2.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Phylo$b[1,1], x=1.25), size = 3.5,alpha=1,stroke=0,shape=19,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Phylo2$ci.lb[1], x=2.25, xend= 2.25, yend= Model_cREML_by_Phylo2$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Phylo$ci.lb[1], x=1.25, xend= 1.25, yend= Model_cREML_by_Phylo$ci.ub[1]), alpha=1,linewidth=1,color='grey30')+
ylab(expression(paste("Effect size (", italic("r"),')')))+xlab('Controlling for phylogeny')+coord_flip()+guides(fill = FALSE, colour = FALSE) +
scale_color_manual(values =colpal4)+
scale_fill_manual(values =colpal4)+
scale_x_discrete(labels=(c("Phylogenetically \n uncorrected ","Phylogenetically \n corrected ")),expand=c(.15,0))+
annotate("text", x=1, y=1.1, label= "n = 23",size=3.5) +
annotate("text", x=2, y=1.1, label= "n = 104",size=3.5) + theme
Figure S2C: Raincloud plot of correlation coefficients for non-phylogenetic and phylogenetic analyses in primary studies including sample sizes and estimates with 95% CI from phylogenetic model (Table 4).
Percentage of species with female-biased SSD
There was variation in the primary studies regarding the typical SSD in the studied taxa (i.e. some studies focused on taxa with more male-biased SSD, while others on taxa with more female-biased SSD). Still, there was no significant relationship between the percentage of species with a female-biased SSD and effect sizes.
MetaData_SSDbias=MetaData
MetaData_SSDbias=MetaData_SSDbias[!is.na(MetaData_SSDbias$SSD_SexBias_in_perc_F),]
Model_REML_SSbias = rma.mv(r ~ SSD_SexBias_in_perc_F, V=Var_r, data = MetaData_SSDbias, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_SSbias)
Multivariate Meta-Analysis Model (k = 108; method: REML)
logLik Deviance AIC BIC AICc
-20.4032 40.8063 50.8063 64.1235 51.4063
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0273 0.1651 67 no Study_ID no
sigma^2.2 0.0453 0.2128 108 no Index no
sigma^2.3 0.0063 0.0794 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 106) = 1489.6620, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 3.6045, p-val = 0.0576
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3864 0.0832 4.6432 <.0001 0.2233 0.5495 ***
SSD_SexBias_in_perc_F -0.0022 0.0011 -1.8985 0.0576 -0.0044 0.0001 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Figure S3A
model_p1 <- predict(Model_REML_SSbias) # Extract model predictions
ggplot(MetaData_SSDbias, aes(x=as.numeric(SSD_SexBias_in_perc_F), y=r)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
theme(axis.text=element_text(size=13),
axis.title=element_text(size=14))+ theme(legend.position="none")+
geom_point(shape=16, size = 3,alpha=0.5)+xlab('% of species with \n female-biased SSD')+ylab(expression(paste("Effect size (", italic("r"),')')))+
geom_hline(yintercept=0, linetype="dashed", color = "black", linewidth=1)+
geom_line( aes(y = model_p1$pred), size = 1)+
geom_ribbon( aes(ymin = model_p1$ci.lb, ymax = model_p1$ci.ub, color = 'black',linetype=NA), alpha = .15,show.legend = F, outline.type = "both") +
theme
Figure S3A: Scatter plot of the sPercentage of species with a female-biased SSD. Dashed line marks a correlation coefficient of zero, black line represents predictions from REML model (controlling for phylogeny) and grey area represents 95% CI on model predictions.
Test for publication bias
To test for publication bias, we transformed r into z scores and ran multilevel mixed-effects models (restricted maximum likelihood) with z as the predictor and its standard error as the response with study ID and an observation level random effect. Models were weight by the mean standard error of z across all studies. While the variance in r depends on the effect size and the sample size, the variance in z is only dependent on the sample size. Hence, if z values correlate with the variance in z, this indicates that small studies were only published, if the effect was large, suggesting publication bias.
Model_REML_PublBias = rma.mv(z ~ SE_z, V=rep((mean(SE_z)*mean(SE_z))*N,length(SE_z)), data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML",control=list(rel.tol=1e-8))
summary(Model_REML_PublBias)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-166.7309 333.4618 343.4618 357.6033 343.9660
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 77 no Study_ID no
sigma^2.2 0.0000 0.0000 127 no Index no
sigma^2.3 0.0000 0.0002 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 125) = 14.7570, p-val = 1.0000
Test of Moderators (coefficient 2):
QM(df = 1) = 1.4805, p-val = 0.2237
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1965 0.1969 0.9978 0.3184 -0.1895 0.5825
SE_z 0.7067 0.5808 1.2168 0.2237 -0.4317 1.8452
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Figure S3B
# Extract model predictions
model_p2 <- predict(Model_REML_PublBias)
ggplot(MetaData, aes(x=SE_z, y=z)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
theme(axis.text=element_text(size=13),
axis.title=element_text(size=14))+ theme(legend.position="none")+ylab(expression(paste("Effect size (", italic("z"),')')))+
geom_point(shape=16, size = 3,alpha=0.5)+xlab('Standard error')+
geom_hline(yintercept=0, linetype="dashed", color = "black", linewidth=1)+
geom_line( aes(y = model_p2$pred), size = 1)+
geom_ribbon( aes(ymin = model_p2$ci.lb, ymax = model_p2$ci.ub, color = 'black',linetype=NA), alpha = .15,show.legend = F, outline.type = "both") +
theme
Figure S3B: Scatter plot of the standard error in z scores of each study against the z score of each primary study (transformed correlation coefficients). Dashed line marks a correlation coefficient of zero, black line represents predictions from REML model (controlling for phylogeny) and grey area represents 95% CI on model predictions.
Publication year
Next we explored the effect of the publication year of each study:
Model_REML_by_Year = rma.mv(r ~ Year, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
summary(Model_REML_by_Year)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.6385 41.2769 51.2769 65.4185 51.7811
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0165 0.1285 77 no Study_ID no
sigma^2.2 0.0462 0.2148 127 no Index no
sigma^2.3 0.0118 0.1087 11 no Class yes
Test for Residual Heterogeneity:
QE(df = 125) = 1344.0152, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 7.9645, p-val = 0.0048
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 12.9132 4.4856 2.8788 0.0040 4.1216 21.7048 **
Year -0.0063 0.0022 -2.8221 0.0048 -0.0107 -0.0019 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Plot publication year (Figure S3C)
Here we plot the publication year:
# Extract model predictions
Model_REML_by_Year = rma.mv(r ~ Year, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index, ~ 1 | Class), R = list(Class = forcedC_Moderators), method = "REML")
model_p4 <- predict(Model_REML_by_Year)
ggplot(MetaData, aes(x=as.numeric(Year), y=r)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
theme(axis.text=element_text(size=13),
axis.title=element_text(size=14))+ theme(legend.position="none")+
geom_point(shape=16, size = 3,alpha=0.5)+xlab('Publication year')+ylab(expression(paste("Effect size (", italic("r"),')')))+
geom_hline(yintercept=0, linetype="dashed", color = "black", linewidth=1)+
geom_line( aes(y = model_p4$pred), size = 1)+
geom_ribbon( aes(ymin = model_p4$ci.lb, ymax = model_p4$ci.ub, color = 'black',linetype=NA), alpha = .15,show.legend = F, outline.type = "both") +
theme
Figure S3C: Scatter plot of correlation coefficients against the publication year of each study. Dashed line marks a correlation coefficient of zero, black line represents predictions from REML model (controlling for phylogeny) and grey area represents 95% CI on model predictions.
Moderator tests for non-phylogenetic models
Here we ran all models without the phylogeny.
Sexual selection episode
The first model explores the effect of the sexual selection episode (i.e. pre-copulatory, post-copulatory or both):
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("pre-copulatory"))
Model_REML_by_cSexSelEpisode = rma.mv(r ~ factor(SexSel_Episode), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelEpisode)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-9.8630 19.7261 29.7261 43.8275 30.2346
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0241 0.1553 77 no Study_ID
sigma^2.2 0.0327 0.1809 127 no Index
Test for Residual Heterogeneity:
QE(df = 124) = 1378.0070, p-val < .0001
Test of Moderators (coefficients 2:3):
QM(df = 2) = 36.8836, p-val < .0001
Model Results:
estimate se zval pval
intrcpt 0.2470 0.0389 6.3556 <.0001
factor(SexSel_Episode)both 0.1694 0.0526 3.2225 0.0013
factor(SexSel_Episode)post-copulatory -0.3236 0.0820 -3.9455 <.0001
ci.lb ci.ub
intrcpt 0.1708 0.3231 ***
factor(SexSel_Episode)both 0.0664 0.2725 **
factor(SexSel_Episode)post-copulatory -0.4843 -0.1628 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("post-copulatory"))
Model_REML_by_cSexSelEpisode2 = rma.mv(r ~ factor(SexSel_Episode), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelEpisode2)
MetaData$SexSel_Episode=relevel(MetaData$SexSel_Episode,c("both"))
Model_REML_by_cSexSelEpisode3 = rma.mv(r ~ factor(SexSel_Episode), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelEpisode3)Finally, we computed FDR corrected p-values:
tab1=as.data.frame(round(p.adjust(c(0.0001, 0.2970, .0001), method = 'fdr'),digit=3),row.names=cbind("Pre-copulatory","Post-copulatory","Both"))
colnames(tab1)<-cbind('P-value')
tab1 P-value
Pre-copulatory 0.000
Post-copulatory 0.297
Both 0.000Sexual selection category
Next we explored the effect of the sexual selection category (i.e. density, mating system, operational sex ratio (OSR), post-mating competition, pre-mating competition, trait-based, other):
MetaData$SexSel_Category=as.factor(MetaData$SexSel_Category)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Postmating competition"))
Model_REML_by_cSexSelCat = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelCat)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
5.8687 -11.7375 6.2625 31.3500 7.8989
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0259 0.1608 77 no Study_ID
sigma^2.2 0.0180 0.1342 127 no Index
Test for Residual Heterogeneity:
QE(df = 120) = 949.5416, p-val < .0001
Test of Moderators (coefficients 2:7):
QM(df = 6) = 87.8526, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.0713 0.0679 -1.0500 0.2937
factor(SexSel_Category)Other 0.4440 0.0935 4.7493 <.0001
factor(SexSel_Category)OSR 0.7075 0.1191 5.9398 <.0001
factor(SexSel_Category)Mating system 0.4795 0.0785 6.1080 <.0001
factor(SexSel_Category)Premating competition 0.4355 0.0817 5.3291 <.0001
factor(SexSel_Category)Density 0.3297 0.1219 2.7042 0.0068
factor(SexSel_Category)Trait-based 0.0902 0.0863 1.0452 0.2959
ci.lb ci.ub
intrcpt -0.2045 0.0618
factor(SexSel_Category)Other 0.2608 0.6272 ***
factor(SexSel_Category)OSR 0.4740 0.9409 ***
factor(SexSel_Category)Mating system 0.3257 0.6334 ***
factor(SexSel_Category)Premating competition 0.2753 0.5956 ***
factor(SexSel_Category)Density 0.0907 0.5687 **
factor(SexSel_Category)Trait-based -0.0790 0.2594
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Trait-based"))
Model_REML_by_cSexSelCat2 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelCat2)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Density"))
Model_REML_by_cSexSelCat3 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelCat3)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Premating competition"))
Model_REML_by_cSexSelCat4 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelCat4)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Mating system"))
Model_REML_by_cSexSelCat5 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelCat5)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("OSR"))
Model_REML_by_cSexSelCat6 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelCat6)
MetaData$SexSel_Category=relevel(MetaData$SexSel_Category,c("Other"))
Model_REML_by_cSexSelCat7 = rma.mv(r ~ factor(SexSel_Category), V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSexSelCat7)Finally, we computed FDR corrected p-values:
tab2=as.data.frame(round(p.adjust(c(0.2937, 0.7259, 0.0110, .0001, .0001, .0001, .0001), method = 'fdr'),digit=3),row.names=cbind("Postmating competition","Trait-based","Density",'Premating competition',"Mating system","OSR","Other"))
colnames(tab2)<-cbind('P-value')
tab2 P-value
Postmating competition 0.343
Trait-based 0.726
Density 0.015
Premating competition 0.000
Mating system 0.000
OSR 0.000
Other 0.000Phylogenetic classes
Next we explored the effect of the phylogenetic classes:
MetaData$Class=relevel(MetaData$Class,c("Actinopterygii"))
Model_cREML_by_Class5 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class5)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2237 0.1104 2.0264 0.0427 0.0073 0.4400 *
ClassAmphibia -0.1024 0.1429 -0.7165 0.4737 -0.3825 0.1777
ClassAnimalia 0.0022 0.2490 0.0089 0.9929 -0.4859 0.4903
ClassAves -0.0041 0.1228 -0.0333 0.9734 -0.2448 0.2366
ClassInsecta 0.1610 0.1596 1.0088 0.3131 -0.1518 0.4737
ClassMalacostraca -0.2215 0.3536 -0.6263 0.5311 -0.9145 0.4716
ClassMammalia 0.1764 0.1197 1.4739 0.1405 -0.0582 0.4109
ClassNematoda 0.1983 0.2666 0.7438 0.4570 -0.3242 0.7207
ClassPisces 0.0346 0.2897 0.1196 0.9048 -0.5332 0.6025
ClassReptilia -0.0616 0.1385 -0.4446 0.6566 -0.3330 0.2099
ClassTrematoda -0.1544 0.2286 -0.6754 0.4994 -0.6025 0.2937
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Amphibia"))
Model_cREML_by_Class4 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class4)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1213 0.0908 1.3357 0.1816 -0.0567 0.2992
ClassActinopterygii 0.1024 0.1429 0.7165 0.4737 -0.1777 0.3825
ClassAnimalia 0.1046 0.2410 0.4341 0.6642 -0.3677 0.5769
ClassAves 0.0983 0.1063 0.9252 0.3549 -0.1099 0.3065
ClassInsecta 0.2634 0.1476 1.7838 0.0745 -0.0260 0.5527 .
ClassMalacostraca -0.1191 0.3480 -0.3422 0.7322 -0.8011 0.5629
ClassMammalia 0.2788 0.1031 2.7032 0.0069 0.0766 0.4809 **
ClassNematoda 0.3007 0.2591 1.1606 0.2458 -0.2071 0.8084
ClassPisces 0.1370 0.2829 0.4845 0.6280 -0.4174 0.6914
ClassReptilia 0.0408 0.1235 0.3305 0.7410 -0.2012 0.2829
ClassTrematoda -0.0520 0.2243 -0.2319 0.8166 -0.4917 0.3876
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Animalia"))
Model_cREML_by_Class11 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class11)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2259 0.2232 1.0119 0.3116 -0.2116 0.6634
ClassAmphibia -0.1046 0.2410 -0.4341 0.6642 -0.5769 0.3677
ClassActinopterygii -0.0022 0.2490 -0.0089 0.9929 -0.4903 0.4859
ClassAves -0.0063 0.2300 -0.0275 0.9781 -0.4570 0.4444
ClassInsecta 0.1587 0.2518 0.6305 0.5284 -0.3347 0.6522
ClassMalacostraca -0.2237 0.4033 -0.5546 0.5792 -1.0142 0.5668
ClassMammalia 0.1741 0.2285 0.7620 0.4460 -0.2738 0.6221
ClassNematoda 0.1961 0.3297 0.5946 0.5521 -0.4502 0.8423
ClassPisces 0.0324 0.3487 0.0930 0.9259 -0.6510 0.7159
ClassReptilia -0.0638 0.2384 -0.2676 0.7890 -0.5311 0.4035
ClassTrematoda -0.1567 0.3032 -0.5167 0.6054 -0.7508 0.4375
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Aves"))
Model_cREML_by_Class = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2196 0.0552 3.9772 <.0001 0.1114 0.3278 ***
ClassAnimalia 0.0063 0.2300 0.0275 0.9781 -0.4444 0.4570
ClassAmphibia -0.0983 0.1063 -0.9252 0.3549 -0.3065 0.1099
ClassActinopterygii 0.0041 0.1228 0.0333 0.9734 -0.2366 0.2448
ClassInsecta 0.1651 0.1286 1.2832 0.1994 -0.0871 0.4172
ClassMalacostraca -0.2174 0.3404 -0.6385 0.5232 -0.8846 0.4499
ClassMammalia 0.1805 0.0728 2.4783 0.0132 0.0377 0.3232 *
ClassNematoda 0.2024 0.2488 0.8133 0.4161 -0.2854 0.6901
ClassPisces 0.0387 0.2735 0.1416 0.8874 -0.4974 0.5748
ClassReptilia -0.0575 0.0996 -0.5773 0.5637 -0.2526 0.1377
ClassTrematoda -0.1503 0.2115 -0.7109 0.4771 -0.5648 0.2641
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Insecta"))
Model_cREML_by_Class8 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class8)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3846 0.1164 3.3033 0.0010 0.1564 0.6128 ***
ClassAves -0.1651 0.1286 -1.2832 0.1994 -0.4172 0.0871
ClassAnimalia -0.1587 0.2518 -0.6305 0.5284 -0.6522 0.3347
ClassAmphibia -0.2634 0.1476 -1.7838 0.0745 -0.5527 0.0260 .
ClassActinopterygii -0.1610 0.1596 -1.0088 0.3131 -0.4737 0.1518
ClassMalacostraca -0.3824 0.3555 -1.0757 0.2821 -1.0792 0.3144
ClassMammalia 0.0154 0.1259 0.1224 0.9026 -0.2314 0.2622
ClassNematoda 0.0373 0.2691 0.1386 0.8897 -0.4902 0.5648
ClassPisces -0.1263 0.2921 -0.4325 0.6654 -0.6988 0.4462
ClassReptilia -0.2225 0.1434 -1.5519 0.1207 -0.5036 0.0585
ClassTrematoda -0.3154 0.2343 -1.3463 0.1782 -0.7746 0.1438
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Malacostraca"))
Model_cREML_by_Class7 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class7)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0022 0.3359 0.0065 0.9948 -0.6562 0.6606
ClassInsecta 0.3824 0.3555 1.0757 0.2821 -0.3144 1.0792
ClassAves 0.2174 0.3404 0.6385 0.5232 -0.4499 0.8846
ClassAnimalia 0.2237 0.4033 0.5546 0.5792 -0.5668 1.0142
ClassAmphibia 0.1191 0.3480 0.3422 0.7322 -0.5629 0.8011
ClassActinopterygii 0.2215 0.3536 0.6263 0.5311 -0.4716 0.9145
ClassMammalia 0.3978 0.3395 1.1719 0.2412 -0.2675 1.0632
ClassNematoda 0.4197 0.4144 1.0129 0.3111 -0.3925 1.2319
ClassPisces 0.2561 0.4297 0.5960 0.5511 -0.5860 1.0982
ClassReptilia 0.1599 0.3462 0.4618 0.6442 -0.5187 0.8384
ClassTrematoda 0.0670 0.3936 0.1703 0.8648 -0.7044 0.8385
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Mammalia"))
Model_cREML_by_Class3 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class3)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.4000 0.0489 8.1767 <.0001 0.3041 0.4959 ***
ClassMalacostraca -0.3978 0.3395 -1.1719 0.2412 -1.0632 0.2675
ClassInsecta -0.0154 0.1259 -0.1224 0.9026 -0.2622 0.2314
ClassAves -0.1805 0.0728 -2.4783 0.0132 -0.3232 -0.0377 *
ClassAnimalia -0.1741 0.2285 -0.7620 0.4460 -0.6221 0.2738
ClassAmphibia -0.2788 0.1031 -2.7032 0.0069 -0.4809 -0.0766 **
ClassActinopterygii -0.1764 0.1197 -1.4739 0.1405 -0.4109 0.0582
ClassNematoda 0.0219 0.2475 0.0885 0.9295 -0.4632 0.5070
ClassPisces -0.1417 0.2723 -0.5204 0.6028 -0.6755 0.3920
ClassReptilia -0.2380 0.0964 -2.4678 0.0136 -0.4269 -0.0490 *
ClassTrematoda -0.3308 0.2092 -1.5810 0.1139 -0.7409 0.0793
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Nematoda"))
Model_cREML_by_Class10 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class10)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.4219 0.2426 1.7389 0.0821 -0.0536 0.8975 .
ClassMammalia -0.0219 0.2475 -0.0885 0.9295 -0.5070 0.4632
ClassMalacostraca -0.4197 0.4144 -1.0129 0.3111 -1.2319 0.3925
ClassInsecta -0.0373 0.2691 -0.1386 0.8897 -0.5648 0.4902
ClassAves -0.2024 0.2488 -0.8133 0.4161 -0.6901 0.2854
ClassAnimalia -0.1961 0.3297 -0.5946 0.5521 -0.8423 0.4502
ClassAmphibia -0.3007 0.2591 -1.1606 0.2458 -0.8084 0.2071
ClassActinopterygii -0.1983 0.2666 -0.7438 0.4570 -0.7207 0.3242
ClassPisces -0.1636 0.3614 -0.4527 0.6508 -0.8721 0.5448
ClassReptilia -0.2599 0.2567 -1.0124 0.3114 -0.7629 0.2432
ClassTrematoda -0.3527 0.3177 -1.1101 0.2670 -0.9754 0.2700
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Pisces"))
Model_cREML_by_Class6 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class6)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2583 0.2679 0.9642 0.3350 -0.2668 0.7834
ClassNematoda 0.1636 0.3614 0.4527 0.6508 -0.5448 0.8721
ClassMammalia 0.1417 0.2723 0.5204 0.6028 -0.3920 0.6755
ClassMalacostraca -0.2561 0.4297 -0.5960 0.5511 -1.0982 0.5860
ClassInsecta 0.1263 0.2921 0.4325 0.6654 -0.4462 0.6988
ClassAves -0.0387 0.2735 -0.1416 0.8874 -0.5748 0.4974
ClassAnimalia -0.0324 0.3487 -0.0930 0.9259 -0.7159 0.6510
ClassAmphibia -0.1370 0.2829 -0.4845 0.6280 -0.6914 0.4174
ClassActinopterygii -0.0346 0.2897 -0.1196 0.9048 -0.6025 0.5332
ClassReptilia -0.0962 0.2807 -0.3428 0.7317 -0.6463 0.4539
ClassTrematoda -0.1891 0.3374 -0.5604 0.5752 -0.8504 0.4722
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Reptilia"))
Model_cREML_by_Class2 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class2)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1621 0.0837 1.9358 0.0529 -0.0020 0.3262 .
ClassPisces 0.0962 0.2807 0.3428 0.7317 -0.4539 0.6463
ClassNematoda 0.2599 0.2567 1.0124 0.3114 -0.2432 0.7629
ClassMammalia 0.2380 0.0964 2.4678 0.0136 0.0490 0.4269 *
ClassMalacostraca -0.1599 0.3462 -0.4618 0.6442 -0.8384 0.5187
ClassInsecta 0.2225 0.1434 1.5519 0.1207 -0.0585 0.5036
ClassAves 0.0575 0.0996 0.5773 0.5637 -0.1377 0.2526
ClassAnimalia 0.0638 0.2384 0.2676 0.7890 -0.4035 0.5311
ClassAmphibia -0.0408 0.1235 -0.3305 0.7410 -0.2829 0.2012
ClassActinopterygii 0.0616 0.1385 0.4446 0.6566 -0.2099 0.3330
ClassTrematoda -0.0928 0.2215 -0.4192 0.6751 -0.5270 0.3413
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MetaData$Class=relevel(MetaData$Class,c("Trematoda"))
Model_cREML_by_Class9 = rma.mv(r ~ Class, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Class9)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-20.5793 41.1586 67.1586 102.9553 70.7272
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 77 no Study_ID
sigma^2.2 0.0492 0.2217 127 no Index
Test for Residual Heterogeneity:
QE(df = 116) = 1256.7618, p-val < .0001
Test of Moderators (coefficients 2:11):
QM(df = 10) = 15.7503, p-val = 0.1070
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0692 0.2051 0.3375 0.7357 -0.3328 0.4713
ClassReptilia 0.0928 0.2215 0.4192 0.6751 -0.3413 0.5270
ClassPisces 0.1891 0.3374 0.5604 0.5752 -0.4722 0.8504
ClassNematoda 0.3527 0.3177 1.1101 0.2670 -0.2700 0.9754
ClassMammalia 0.3308 0.2092 1.5810 0.1139 -0.0793 0.7409
ClassMalacostraca -0.0670 0.3936 -0.1703 0.8648 -0.8385 0.7044
ClassInsecta 0.3154 0.2343 1.3463 0.1782 -0.1438 0.7746
ClassAves 0.1503 0.2115 0.7109 0.4771 -0.2641 0.5648
ClassAnimalia 0.1567 0.3032 0.5167 0.6054 -0.4375 0.7508
ClassAmphibia 0.0520 0.2243 0.2319 0.8166 -0.3876 0.4917
ClassActinopterygii 0.1544 0.2286 0.6754 0.4994 -0.2937 0.6025
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Finally, we computed FDR corrected p-values:
tab2=as.data.frame(round(p.adjust(c(.0001, 0.0529, .0001, 0.1816, 0.0427, 0.3350, 0.9948, 0.0010, 0.7357, 0.0821, 0.3116), method = 'fdr'),digit=3),row.names=cbind( "Animalia","Nematoda","Trematoda","Insecta","Malacostraca","Pisces" ,"Actinopterygii","Amphibia","Mammalia","Reptilia","Aves"))
colnames(tab2)<-cbind('P-value')
tab2 P-value
Animalia 0.001
Nematoda 0.116
Trematoda 0.001
Insecta 0.285
Malacostraca 0.116
Pisces 0.409
Actinopterygii 0.995
Amphibia 0.004
Mammalia 0.809
Reptilia 0.151
Aves 0.409Plot phylogenetic classes (Figure 1)
Here we plot the phylogenetic classes:
MetaData$Class=as.factor(MetaData$Class)
MetaData$Class=factor(MetaData$Class, levels = (c("Animalia","Nematoda","Trematoda","Insecta","Malacostraca","Pisces" ,"Actinopterygii","Amphibia","Mammalia","Reptilia","Aves")))
MetaData_sorted=MetaData
MetaData_sorted$Class[MetaData_sorted$Class=='Pisces']='Actinopterygii'
ggplot() +
geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=1)+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class$ci.lb[1], x=11.35, xend= 11.35, yend= Model_cREML_by_Class$ci.ub[1]),linewidth=1,color=colpal5[11], alpha=1)+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class2$ci.lb[1], x=10.35, xend= 10.35, yend= Model_cREML_by_Class2$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[10])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class3$ci.lb[1], x=9.35, xend= 9.35, yend= Model_cREML_by_Class3$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[9])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class4$ci.lb[1], x=8.35, xend= 8.35, yend= Model_cREML_by_Class4$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[8])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class5$ci.lb[1], x=7.35, xend= 7.35, yend= Model_cREML_by_Class5$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[7])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class6$ci.lb[1], x=6.35, xend= 6.35, yend= Model_cREML_by_Class6$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[6])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class7$ci.lb[1], x=5.35, xend= 5.35, yend= Model_cREML_by_Class7$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[5])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class8$ci.lb[1], x=4.35, xend= 4.35, yend= Model_cREML_by_Class8$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[4])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class9$ci.lb[1], x=3.35, xend= 3.35, yend= Model_cREML_by_Class9$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[3])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class10$ci.lb[1], x=2.35, xend= 2.35, yend= Model_cREML_by_Class10$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[2])+
geom_segment(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class11$ci.lb[1], x=1.35, xend= 1.35, yend= Model_cREML_by_Class11$ci.ub[1]), alpha=1,linewidth=1,color=colpal5[1])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class$b[1,1], x=11.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class2$b[1,1], x=10.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class3$b[1,1], x=9.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class4$b[1,1], x=8.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class5$b[1,1], x=7.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class6$b[1,1], x=6.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class7$b[1,1], x=5.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class8$b[1,1], x=4.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class9$b[1,1], x=3.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class10$b[1,1], x=2.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class11$b[1,1], x=1.35), size = 3.5,alpha=1,stroke=0,shape=19,color='white')+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class$b[1,1], x=11.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[11])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class2$b[1,1], x=10.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[10])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class3$b[1,1], x=9.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[9])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class4$b[1,1], x=8.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[8])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class5$b[1,1], x=7.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[7])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class6$b[1,1], x=6.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[6])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class7$b[1,1], x=5.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[5])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class8$b[1,1], x=4.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[4])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class9$b[1,1], x=3.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[3])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class10$b[1,1], x=2.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[2])+
geom_point(inherit.aes = F,mapping = aes(y=Model_cREML_by_Class11$b[1,1], x=1.35), size = 3.5,alpha=1,stroke=0,shape=19,color=colpal5[1])+
ylab(expression(paste("Effect size (", italic("r"),')')))+xlab('Class')+coord_flip()+guides(fill = FALSE, colour = FALSE) +
scale_x_discrete(labels=(c("Animalia","Nematoda","Trematoda","Insecta","Malacostraca","Pisces" ,"Actinopterygii","Amphibia","Mammalia","Reptilia","Aves")),expand=c(.1,0))+
geom_point(data=MetaData,aes(x=Class, y=r, fill = Class, colour = Class),position = position_jitter(width = .1), size = 2.5,alpha=0.6,stroke=0,shape=19)+
scale_color_manual(values =colpal5)+
scale_fill_manual(values =colpal5)+
theme
Figure 1: Phylogeny based on classes including sample sizes for the number of studies and number of effect sizes, respectively, and estimates with 95%CI from non-phylogenetic model.
Type of SSD measure
Next we explored the effect of the type of SSD measure (i.e. body mass or size):
MetaData$SSD_Proxy=as.factor(MetaData$SSD_Proxy)
MetaData$SSD_Proxy=relevel(MetaData$SSD_Proxy,c("Body mass"))
Model_REML_by_cSSDMeasure = rma.mv(r ~ SSD_Proxy, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSSDMeasure)
Multivariate Meta-Analysis Model (k = 125; method: REML)
logLik Deviance AIC BIC AICc
-24.0901 48.1801 56.1801 67.4289 56.5191
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0268 0.1638 76 no Study_ID
sigma^2.2 0.0465 0.2156 125 no Index
Test for Residual Heterogeneity:
QE(df = 123) = 1553.9370, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 4.6542, p-val = 0.0310
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3484 0.0447 7.7937 <.0001 0.2608 0.4360 ***
SSD_ProxyBody size -0.1253 0.0581 -2.1574 0.0310 -0.2391 -0.0115 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$SSD_Proxy=relevel(MetaData$SSD_Proxy,c("Body size"))
Model_REML_by_cSSDMeasure2 = rma.mv(r ~ SSD_Proxy, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cSSDMeasure2)Finally, we computed FDR corrected p-values:
tab3=as.data.frame(round(p.adjust(c(.0001, .0001), method = 'fdr'),digit=3),row.names=cbind("Body mass","Body size"))
colnames(tab3)<-cbind('P-value')
tab3 P-value
Body mass 0
Body size 0SSD measure controlled for body size?
Next we explored the effect if the primary study controlled the SSD for body size (i.e. uncontrolled or controlled):
MetaData$BodySizeControlled=as.factor(MetaData$BodySizeControlled)
MetaData$BodySizeControlled=relevel(MetaData$BodySizeControlled,c("No"))
Model_REML_by_cBodySizeCont = rma.mv(r ~ BodySizeControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cBodySizeCont)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-25.2183 50.4365 58.4365 69.7498 58.7699
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0308 0.1755 77 no Study_ID
sigma^2.2 0.0455 0.2132 127 no Index
Test for Residual Heterogeneity:
QE(df = 125) = 1601.6900, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 1.2317, p-val = 0.2671
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2907 0.0336 8.6540 <.0001 0.2249 0.3566 ***
BodySizeControlledYes -0.0971 0.0875 -1.1098 0.2671 -0.2687 0.0744
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$BodySizeControlled=relevel(MetaData$BodySizeControlled,c("Yes"))
Model_REML_by_cBodySizeCont2 = rma.mv(r ~ BodySizeControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_REML_by_cBodySizeCont2)Finally, we computed FDR corrected p-values:
tab4=as.data.frame(round(p.adjust(c(.0001, 0.0166), method = 'fdr'),digit=3),row.names=cbind("uncontrolled","controlled"))
colnames(tab4)<-cbind('P-value')
tab4 P-value
uncontrolled 0.000
controlled 0.017Phylogentic vs non-phylogenetic studies
Next we explored the effect of studies controlling for phylogeny vs non-phylogenetic studies:
MetaData$PhyloControlled=as.factor(MetaData$PhyloControlled)
MetaData$PhyloControlled=relevel(MetaData$PhyloControlled,c("No"))
Model_cREML_by_cPhylo = rma.mv(r ~ PhyloControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_cPhylo)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-22.5511 45.1022 53.1022 64.4155 53.4355
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0234 0.1529 77 no Study_ID
sigma^2.2 0.0475 0.2180 127 no Index
Test for Residual Heterogeneity:
QE(df = 125) = 1428.7953, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 6.9580, p-val = 0.0083
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.4463 0.0716 6.2341 <.0001 0.3060 0.5866 ***
PhyloControlledYes -0.2073 0.0786 -2.6378 0.0083 -0.3614 -0.0533 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We then re-leveled the model for post-hoc comparisons:
MetaData$PhyloControlled=relevel(MetaData$PhyloControlled,c("Yes"))
Model_cREML_by_cPhylo2 = rma.mv(r ~ PhyloControlled, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_cPhylo2)Finally, we computed FDR corrected p-values:
tab5=as.data.frame(round(p.adjust(c(0.0001, 0.0001), method = 'fdr'),digit=3),row.names=cbind("Non-phylogenetic","With phylogeny"))
colnames(tab5)<-cbind('P-value')
tab5 P-value
Non-phylogenetic 0
With phylogeny 0Percentage of species with female-biased SSD
Model_cREML_SSbias = rma.mv(r ~ SSD_SexBias_in_perc_F, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_SSbias)
Multivariate Meta-Analysis Model (k = 108; method: REML)
logLik Deviance AIC BIC AICc
-20.7304 41.4609 49.4609 60.1146 49.8569
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0315 0.1774 67 no Study_ID
sigma^2.2 0.0447 0.2115 108 no Index
Test for Residual Heterogeneity:
QE(df = 106) = 1489.6620, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 5.7675, p-val = 0.0163
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3944 0.0532 7.4116 <.0001 0.2901 0.4988 ***
SSD_SexBias_in_perc_F -0.0025 0.0010 -2.4016 0.0163 -0.0045 -0.0005 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test for publication bias
Model_cREML_PublBias = rma.mv(z ~ SE_z, V=rep(((mean(SE_z)*mean(SE_z))*N),length(SE_z)), data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_PublBias)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-166.7309 333.4618 341.4618 352.7750 341.7951
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 77 no Study_ID
sigma^2.2 0.0000 0.0000 127 no Index
Test for Residual Heterogeneity:
QE(df = 125) = 14.7570, p-val = 1.0000
Test of Moderators (coefficient 2):
QM(df = 1) = 1.4805, p-val = 0.2237
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1965 0.1969 0.9978 0.3184 -0.1895 0.5825
SE_z 0.7067 0.5808 1.2168 0.2237 -0.4317 1.8452
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Publication year
Next we explored the effect of the publication year of each study:
Model_cREML_by_Year = rma.mv(r ~ Year, V=Var_r, data = MetaData, random = c(~ 1 | Study_ID,~ 1 | Index), method = "REML")
summary(Model_cREML_by_Year)
Multivariate Meta-Analysis Model (k = 127; method: REML)
logLik Deviance AIC BIC AICc
-22.2585 44.5171 52.5171 63.8303 52.8504
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0270 0.1643 77 no Study_ID
sigma^2.2 0.0441 0.2100 127 no Index
Test for Residual Heterogeneity:
QE(df = 125) = 1344.0152, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 7.4667, p-val = 0.0063
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 13.4012 4.8035 2.7899 0.0053 3.9865 22.8159 **
Year -0.0065 0.0024 -2.7325 0.0063 -0.0112 -0.0018 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Forest plot
# Sort data by Class
MetaData$Class=as.factor(MetaData$Class)
MetaData$Class=factor(MetaData$Class, levels = c("Animalia","Nematoda","Trematoda","Insecta","Malacostraca","Pisces" ,"Actinopterygii","Amphibia","Mammalia","Reptilia","Aves"))
MetaData_sorted <-MetaData[order(MetaData$r,decreasing = T),]
MetaData_sorted$Class[MetaData_sorted$Class=='Pisces']='Actinopterygii'
MetaData_sorted <-MetaData_sorted[order(MetaData_sorted$Class),]
# Add global effect size to data
forestData=MetaData_sorted[,c(3,5,10,15,18,19)]
forestData[,c(4:6)]=lapply(forestData[,c(4:6)],as.numeric)
forestData$SexSel_Episode=as.factor(forestData$SexSel_Episode)
forestData$SexSel_Episode=factor(forestData$SexSel_Episode, levels = c("pre-copulatory","post-copulatory" ,"both"))Figure 2
ggplot(data=forestData, aes(y=1:nrow(forestData), x=r, xmin=lCI, xmax=uCI,color=Class,shape=SexSel_Episode,fill=SexSel_Episode)) +
geom_vline(xintercept=Model_REML_Null$b[1,1], color='black', linetype='solid', alpha=.9,linewidth=0.75)+
annotate("rect", xmin = Model_REML_Null$ci.lb[1], xmax = Model_REML_Null$ci.ub[1], ymin = 1, ymax = nrow(forestData),alpha = .25)+
geom_vline(xintercept=0, color='black', linetype='dashed', alpha=.5,linewidth=0.8) +
geom_errorbarh(height=0,size = ifelse(1:nrow(forestData) == c(1), 3.3, 3.3), alpha=1,colour='white')+
geom_errorbarh(height=0,size = ifelse(1:nrow(forestData) == c(1), 3.3, 3.3), alpha=1)+
geom_point(size = ifelse(1:nrow(forestData) == c(1), 2, 2), alpha=1,stroke=1,colour='white')+
geom_point(size = ifelse(1:nrow(forestData) == c(1), 2, 2), alpha=1,stroke=1,colour='black')+
scale_y_continuous(breaks=seq(1, nrow(forestData), by=1), labels=forestData$AuthorsAndYear,limits=c(1,length(forestData$AuthorsAndYear)),expand=c(0.015,0.015)) +
labs(title='', x=expression(paste(italic("r "),'(95% CI)')), y = 'Study') +
theme_classic()+theme(axis.text.x=element_text(size=12),
axis.title=element_text(size=12))+ labs(fill = c("Sexual selection episode"),shape = c("Sexual selection mode"))+
scale_colour_manual(values=colpal6, labels=NULL,guide = "none")+
scale_shape_manual(values=c(21,21,21) ,labels=c("Pre-copulatory","Post-copulatory" ,"Both"))+
scale_fill_manual(values=c('grey50','white','black'),labels=c("Pre-copulatory","Post-copulatory" ,"Both"))
Figure 2: Forest plot including correlation coefficients (±95% CI), as well as the global effect size (in grey) with and without controlling for phylogeny. Different animal classes highlighted by background colors (see Figure 1 for class names) and sexual selection modes pre- (red), post-copulatory (blue) and both (green) by color of effect size.
R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] gridExtra_2.3 wesanderson_0.3.7 PupillometryR_0.0.5
[4] rlang_1.1.3 dplyr_1.1.4 cowplot_1.1.3
[7] ggdist_3.3.2 MCMCglmm_2.36 coda_0.19-4.1
[10] RColorBrewer_1.1-3 ggplot2_3.5.1 data.table_1.15.4
[13] multcomp_1.4-26 TH.data_1.1-2 survival_3.5-8
[16] mvtnorm_1.2-4 psych_2.4.6.26 pwr_1.3-0
[19] MASS_7.3-60.2 metafor_4.6-0 numDeriv_2016.8-1.1
[22] metadat_1.2-0 Matrix_1.7-0 matrixcalc_1.0-6
[25] ape_5.8 workflowr_1.7.1
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 farver_2.1.2 fastmap_1.1.1
[4] tensorA_0.36.2.1 mathjaxr_1.6-0 promises_1.3.0
[7] digest_0.6.35 lifecycle_1.0.4 processx_3.8.4
[10] magrittr_2.0.3 compiler_4.4.0 sass_0.4.9
[13] tools_4.4.0 utf8_1.2.4 yaml_2.3.8
[16] knitr_1.46 labeling_0.4.3 mnormt_2.1.1
[19] withr_3.0.0 purrr_1.0.2 grid_4.4.0
[22] fansi_1.0.6 git2r_0.33.0 colorspace_2.1-0
[25] scales_1.3.0 cli_3.6.2 rmarkdown_2.26
[28] generics_0.1.3 rstudioapi_0.16.0 httr_1.4.7
[31] cachem_1.0.8 stringr_1.5.1 splines_4.4.0
[34] parallel_4.4.0 vctrs_0.6.5 sandwich_3.1-0
[37] jsonlite_1.8.8 callr_3.7.6 jquerylib_0.1.4
[40] tidyr_1.3.1 glue_1.7.0 codetools_0.2-20
[43] ps_1.7.6 distributional_0.4.0 stringi_1.8.4
[46] cubature_2.1.0 gtable_0.3.5 later_1.3.2
[49] munsell_0.5.1 tibble_3.2.1 pillar_1.9.0
[52] htmltools_0.5.8.1 R6_2.5.1 rprojroot_2.0.4
[55] evaluate_0.23 lattice_0.22-6 highr_0.10
[58] corpcor_1.6.10 httpuv_1.6.15 bslib_0.7.0
[61] Rcpp_1.0.12 nlme_3.1-164 whisker_0.4.1
[64] xfun_0.43 fs_1.6.4 zoo_1.8-12
[67] getPass_0.2-4 pkgconfig_2.0.3