Last updated: 2024-05-01
Checks: 7 0
Knit directory:
SSD_and_sexual_selection_2023/
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 7fd037e. 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
Untracked files:
Untracked: data/Data_SexSelSSD.csv
Untracked: data/Pyologeny_SexSelSSD.txt
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/index2.Rmd) and HTML
(docs/index2.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 |
|---|---|---|---|---|
| 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 | 9c4346c | LennartWinkler | 2023-05-03 | wflow_publish(republish = TRUE, all = T) |
| html | 9c4346c | LennartWinkler | 2023-05-03 | wflow_publish(republish = TRUE, all = T) |
| 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 | 57ca562 | LennartWinkler | 2023-04-30 | update |
| Rmd | 7a0f67c | LennartWinkler | 2023-04-30 | Set up 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 correct for phylogenetic non-independence (see Supporting Information). # 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
# Remove studies that did not correct for phylogenetic non-independence
MetaData=MetaData[MetaData$PhyloControlled=='Yes',]
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)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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-19.4865 38.9730 46.9730 57.4335 47.3897
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 59 no Study_ID no
sigma^2.2 0.0641 0.2531 102 no Index no
sigma^2.3 0.0211 0.1454 10 no Class yes
Test for Heterogeneity:
Q(df = 101) = 1268.8327, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2077 0.0853 2.4341 0.0149 0.0405 0.3750 *
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-22.6630 45.3260 51.3260 59.1714 51.5734
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1421 59 no Study_ID
sigma^2.2 0.0585 0.2419 102 no Index
Test for Heterogeneity:
Q(df = 101) = 1268.8327, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2283 0.0339 6.7344 <.0001 0.1619 0.2948 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next, we ran a series of models that test the effect of different moderators. Again we started with models including the phylogeny.
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_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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-5.0927 10.1854 22.1854 37.7561 23.0984
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0059 0.0767 59 no Study_ID no
sigma^2.2 0.0382 0.1954 102 no Index no
sigma^2.3 0.0213 0.1461 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 99) = 1029.8950, p-val < .0001
Test of Moderators (coefficients 2:3):
QM(df = 2) = 35.8112, p-val < .0001
Model Results:
estimate se zval pval
intrcpt 0.2312 0.0876 2.6399 0.0083
factor(SexSel_Episode)post-copulatory -0.3185 0.0774 -4.1136 <.0001
factor(SexSel_Episode)both 0.1564 0.0566 2.7619 0.0057
ci.lb ci.ub
intrcpt 0.0595 0.4028 **
factor(SexSel_Episode)post-copulatory -0.4703 -0.1668 ***
factor(SexSel_Episode)both 0.0454 0.2673 **
---
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.0083, 0.3897, .0001), method = 'fdr'),digit=3),row.names=cbind("Pre-copulatory","Post-copulatory","Both"))
colnames(tab1)<-cbind('P-value')
tab1 P-value
Pre-copulatory 0.012
Post-copulatory 0.390
Both 0.000Next 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
4.1746 -8.3492 11.6508 37.1895 14.2698
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0144 0.1198 59 no Study_ID no
sigma^2.2 0.0229 0.1513 102 no Index no
sigma^2.3 0.0166 0.1288 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 95) = 739.1552, p-val < .0001
Test of Moderators (coefficients 2:7):
QM(df = 6) = 69.8684, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.1255 0.0965 -1.3007 0.1933
factor(SexSel_Category)Density 0.3295 0.1226 2.6873 0.0072
factor(SexSel_Category)Mating system 0.4829 0.0827 5.8426 <.0001
factor(SexSel_Category)OSR 0.7280 0.1356 5.3697 <.0001
factor(SexSel_Category)Other 0.4671 0.0960 4.8656 <.0001
factor(SexSel_Category)Premating competition 0.4140 0.0837 4.9482 <.0001
factor(SexSel_Category)Trait-based 0.1316 0.0899 1.4643 0.1431
ci.lb ci.ub
intrcpt -0.3145 0.0636
factor(SexSel_Category)Density 0.0892 0.5698 **
factor(SexSel_Category)Mating system 0.3209 0.6449 ***
factor(SexSel_Category)OSR 0.4623 0.9937 ***
factor(SexSel_Category)Other 0.2789 0.6552 ***
factor(SexSel_Category)Premating competition 0.2500 0.5780 ***
factor(SexSel_Category)Trait-based -0.0445 0.3077
---
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.1933, 0.9494, 0.1121, 0.0013, 0.0001, .0001, 0.0004), 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.226
Trait-based 0.949
Density 0.157
Premating competition 0.002
Mating system 0.000
OSR 0.000
Other 0.001Next 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 = 100; method: REML)
logLik Deviance AIC BIC AICc
-18.7750 37.5500 47.5500 60.4748 48.2021
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 58 no Study_ID no
sigma^2.2 0.0670 0.2588 100 no Index no
sigma^2.3 0.0079 0.0891 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 98) = 1210.4699, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 4.4689, p-val = 0.0345
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3183 0.0809 3.9349 <.0001 0.1597 0.4768 ***
SSD_ProxyBody size -0.1471 0.0696 -2.1140 0.0345 -0.2835 -0.0107 *
---
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, 0.0082), method = 'fdr'),digit=3),row.names=cbind("Body mass","Body size"))
colnames(tab3)<-cbind('P-value')
tab3 P-value
Body mass 0.000
Body size 0.008Next 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-19.5742 39.1485 49.1485 62.1743 49.7868
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 59 no Study_ID no
sigma^2.2 0.0648 0.2545 102 no Index no
sigma^2.3 0.0203 0.1426 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 100) = 1267.1170, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 0.3170, p-val = 0.5734
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2179 0.0862 2.5265 0.0115 0.0489 0.3869 *
BodySizeControlledYes -0.0492 0.0874 -0.5631 0.5734 -0.2206 0.1221
---
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.0115, 0.1222), method = 'fdr'),digit=3),row.names=cbind("uncontrolled","controlled"))
colnames(tab4)<-cbind('P-value')
tab4 P-value
uncontrolled 0.023
controlled 0.122There 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 = 89; method: REML)
logLik Deviance AIC BIC AICc
-17.4392 34.8784 44.8784 57.2080 45.6192
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0068 0.0824 51 no Study_ID no
sigma^2.2 0.0603 0.2455 89 no Index no
sigma^2.3 0.0181 0.1344 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 87) = 1190.6640, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 1.4926, p-val = 0.2218
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3063 0.1099 2.7869 0.0053 0.0909 0.5217 **
SSD_SexBias_in_perc_F -0.0017 0.0014 -1.2217 0.2218 -0.0043 0.0010
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1To 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-119.9984 239.9969 249.9969 263.0227 250.6352
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 59 no Study_ID no
sigma^2.2 0.0000 0.0000 102 no Index no
sigma^2.3 0.0000 0.0003 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 100) = 17.9544, p-val = 1.0000
Test of Moderators (coefficient 2):
QM(df = 1) = 1.8587, p-val = 0.1728
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0727 0.2243 0.3243 0.7457 -0.3669 0.5124
SE_z 1.1696 0.8579 1.3633 0.1728 -0.5118 2.8510
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-19.6766 39.3532 49.3532 62.3791 49.9915
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 59 no Study_ID no
sigma^2.2 0.0647 0.2543 102 no Index no
sigma^2.3 0.0211 0.1452 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 100) = 1162.1943, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 0.1866, p-val = 0.6658
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 3.1624 6.8405 0.4623 0.6439 -10.2447 16.5695
Year -0.0015 0.0034 -0.4320 0.6658 -0.0081 0.0052
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here we ran all models without the phylogeny.
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-9.5144 19.0288 29.0288 42.0044 29.6740
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0198 0.1407 59 no Study_ID
sigma^2.2 0.0390 0.1974 102 no Index
Test for Residual Heterogeneity:
QE(df = 99) = 1029.8950, p-val < .0001
Test of Moderators (coefficients 2:3):
QM(df = 2) = 31.2049, p-val < .0001
Model Results:
estimate se zval pval
intrcpt 0.2016 0.0431 4.6795 <.0001
factor(SexSel_Episode)both 0.1959 0.0606 3.2334 0.0012
factor(SexSel_Episode)post-copulatory -0.2767 0.0841 -3.2922 0.0010
ci.lb ci.ub
intrcpt 0.1172 0.2861 ***
factor(SexSel_Episode)both 0.0771 0.3146 **
factor(SexSel_Episode)post-copulatory -0.4415 -0.1120 ***
---
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.3058, .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.306
Both 0.000Next 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
2.1308 -4.2615 13.7385 36.7234 15.8561
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0253 0.1592 59 no Study_ID
sigma^2.2 0.0218 0.1476 102 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 739.1552, p-val < .0001
Test of Moderators (coefficients 2:7):
QM(df = 6) = 71.1211, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.0745 0.0694 -1.0729 0.2833
factor(SexSel_Category)Other 0.4632 0.1004 4.6130 <.0001
factor(SexSel_Category)OSR 0.7080 0.1378 5.1368 <.0001
factor(SexSel_Category)Mating system 0.4649 0.0841 5.5272 <.0001
factor(SexSel_Category)Premating competition 0.3890 0.0884 4.4006 <.0001
factor(SexSel_Category)Density 0.3268 0.1266 2.5810 0.0099
factor(SexSel_Category)Trait-based 0.0731 0.0901 0.8116 0.4170
ci.lb ci.ub
intrcpt -0.2105 0.0616
factor(SexSel_Category)Other 0.2664 0.6600 ***
factor(SexSel_Category)OSR 0.4379 0.9781 ***
factor(SexSel_Category)Mating system 0.3001 0.6298 ***
factor(SexSel_Category)Premating competition 0.2158 0.5623 ***
factor(SexSel_Category)Density 0.0786 0.5749 **
factor(SexSel_Category)Trait-based -0.1035 0.2497
---
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.2833, 0.9816, 0.0176, .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.331
Trait-based 0.982
Density 0.025
Premating competition 0.000
Mating system 0.000
OSR 0.000
Other 0.000Next 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1906 0.0988 1.9290 0.0537 -0.0031 0.3843 .
ClassAmphibia -0.1204 0.1351 -0.8913 0.3728 -0.3852 0.1444
ClassAnimalia 0.0353 0.2173 0.1624 0.8710 -0.3906 0.4612
ClassAves -0.0396 0.1120 -0.3541 0.7233 -0.2591 0.1798
ClassInsecta 0.0435 0.1584 0.2745 0.7837 -0.2670 0.3540
ClassMammalia 0.1971 0.1103 1.7872 0.0739 -0.0191 0.4132 .
ClassNematoda 0.2315 0.2372 0.9759 0.3291 -0.2334 0.6965
ClassPisces 0.0677 0.2775 0.2440 0.8073 -0.4762 0.6116
ClassReptilia -0.1769 0.1373 -1.2881 0.1977 -0.4460 0.0923
ClassTrematoda -0.2027 0.2292 -0.8844 0.3765 -0.6518 0.2465
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0702 0.0922 0.7616 0.4463 -0.1104 0.2508
ClassActinopterygii 0.1204 0.1351 0.8913 0.3728 -0.1444 0.3852
ClassAnimalia 0.1557 0.2144 0.7264 0.4676 -0.2644 0.5759
ClassAves 0.0808 0.1061 0.7611 0.4466 -0.1273 0.2888
ClassInsecta 0.1639 0.1544 1.0618 0.2883 -0.1386 0.4665
ClassMammalia 0.3175 0.1043 3.0427 0.0023 0.1130 0.5220 **
ClassNematoda 0.3519 0.2345 1.5006 0.1335 -0.1077 0.8116
ClassPisces 0.1881 0.2752 0.6836 0.4942 -0.3512 0.7275
ClassReptilia -0.0564 0.1326 -0.4256 0.6704 -0.3163 0.2034
ClassTrematoda -0.0822 0.2264 -0.3633 0.7164 -0.5259 0.3614
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2259 0.1936 1.1671 0.2432 -0.1535 0.6053
ClassAmphibia -0.1557 0.2144 -0.7264 0.4676 -0.5759 0.2644
ClassActinopterygii -0.0353 0.2173 -0.1624 0.8710 -0.4612 0.3906
ClassAves -0.0749 0.2006 -0.3736 0.7087 -0.4681 0.3182
ClassInsecta 0.0082 0.2298 0.0356 0.9716 -0.4422 0.4586
ClassMammalia 0.1618 0.1996 0.8103 0.4178 -0.2295 0.5531
ClassNematoda 0.1962 0.2898 0.6771 0.4983 -0.3717 0.7642
ClassPisces 0.0324 0.3236 0.1001 0.9202 -0.6018 0.6666
ClassReptilia -0.2121 0.2158 -0.9832 0.3255 -0.6350 0.2107
ClassTrematoda -0.2380 0.2832 -0.8402 0.4008 -0.7931 0.3171
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1510 0.0527 2.8661 0.0042 0.0477 0.2542 **
ClassAnimalia 0.0749 0.2006 0.3736 0.7087 -0.3182 0.4681
ClassAmphibia -0.0808 0.1061 -0.7611 0.4466 -0.2888 0.1273
ClassActinopterygii 0.0396 0.1120 0.3541 0.7233 -0.1798 0.2591
ClassInsecta 0.0831 0.1346 0.6177 0.5368 -0.1806 0.3469
ClassMammalia 0.2367 0.0719 3.2921 0.0010 0.0958 0.3776 ***
ClassNematoda 0.2712 0.2220 1.2214 0.2219 -0.1640 0.7063
ClassPisces 0.1073 0.2646 0.4057 0.6850 -0.4113 0.6259
ClassReptilia -0.1372 0.1089 -1.2598 0.2078 -0.3507 0.0763
ClassTrematoda -0.1630 0.2134 -0.7641 0.4448 -0.5812 0.2552
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2341 0.1238 1.8903 0.0587 -0.0086 0.4768 .
ClassAves -0.0831 0.1346 -0.6177 0.5368 -0.3469 0.1806
ClassAnimalia -0.0082 0.2298 -0.0356 0.9716 -0.4586 0.4422
ClassAmphibia -0.1639 0.1544 -1.0618 0.2883 -0.4665 0.1386
ClassActinopterygii -0.0435 0.1584 -0.2745 0.7837 -0.3540 0.2670
ClassMammalia 0.1536 0.1332 1.1534 0.2488 -0.1074 0.4146
ClassNematoda 0.1880 0.2487 0.7561 0.4496 -0.2994 0.6754
ClassPisces 0.0242 0.2874 0.0843 0.9328 -0.5390 0.5874
ClassReptilia -0.2203 0.1563 -1.4098 0.1586 -0.5266 0.0860
ClassTrematoda -0.2461 0.2410 -1.0213 0.3071 -0.7185 0.2262
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3877 0.0489 7.9201 <.0001 0.2917 0.4836 ***
ClassInsecta -0.1536 0.1332 -1.1534 0.2488 -0.4146 0.1074
ClassAves -0.2367 0.0719 -3.2921 0.0010 -0.3776 -0.0958 ***
ClassAnimalia -0.1618 0.1996 -0.8103 0.4178 -0.5531 0.2295
ClassAmphibia -0.3175 0.1043 -3.0427 0.0023 -0.5220 -0.1130 **
ClassActinopterygii -0.1971 0.1103 -1.7872 0.0739 -0.4132 0.0191 .
ClassNematoda 0.0344 0.2211 0.1558 0.8762 -0.3990 0.4679
ClassPisces -0.1294 0.2639 -0.4903 0.6239 -0.6466 0.3878
ClassReptilia -0.3739 0.1072 -3.4892 0.0005 -0.5840 -0.1639 ***
ClassTrematoda -0.3997 0.2125 -1.8813 0.0599 -0.8162 0.0167 .
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.4221 0.2157 1.9573 0.0503 -0.0006 0.8448 .
ClassMammalia -0.0344 0.2211 -0.1558 0.8762 -0.4679 0.3990
ClassInsecta -0.1880 0.2487 -0.7561 0.4496 -0.6754 0.2994
ClassAves -0.2712 0.2220 -1.2214 0.2219 -0.7063 0.1640
ClassAnimalia -0.1962 0.2898 -0.6771 0.4983 -0.7642 0.3717
ClassAmphibia -0.3519 0.2345 -1.5006 0.1335 -0.8116 0.1077
ClassActinopterygii -0.2315 0.2372 -0.9759 0.3291 -0.6965 0.2334
ClassPisces -0.1638 0.3373 -0.4857 0.6272 -0.8248 0.4972
ClassReptilia -0.4084 0.2358 -1.7319 0.0833 -0.8705 0.0538 .
ClassTrematoda -0.4342 0.2988 -1.4532 0.1462 -1.0197 0.1514
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2583 0.2593 0.9962 0.3192 -0.2499 0.7665
ClassNematoda 0.1638 0.3373 0.4857 0.6272 -0.4972 0.8248
ClassMammalia 0.1294 0.2639 0.4903 0.6239 -0.3878 0.6466
ClassInsecta -0.0242 0.2874 -0.0843 0.9328 -0.5874 0.5390
ClassAves -0.1073 0.2646 -0.4057 0.6850 -0.6259 0.4113
ClassAnimalia -0.0324 0.3236 -0.1001 0.9202 -0.6666 0.6018
ClassAmphibia -0.1881 0.2752 -0.6836 0.4942 -0.7275 0.3512
ClassActinopterygii -0.0677 0.2775 -0.2440 0.8073 -0.6116 0.4762
ClassReptilia -0.2445 0.2763 -0.8852 0.3761 -0.7860 0.2969
ClassTrematoda -0.2704 0.3316 -0.8152 0.4149 -0.9204 0.3796
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0138 0.0953 0.1443 0.8853 -0.1731 0.2006
ClassPisces 0.2445 0.2763 0.8852 0.3761 -0.2969 0.7860
ClassNematoda 0.4084 0.2358 1.7319 0.0833 -0.0538 0.8705 .
ClassMammalia 0.3739 0.1072 3.4892 0.0005 0.1639 0.5840 ***
ClassInsecta 0.2203 0.1563 1.4098 0.1586 -0.0860 0.5266
ClassAves 0.1372 0.1089 1.2598 0.2078 -0.0763 0.3507
ClassAnimalia 0.2121 0.2158 0.9832 0.3255 -0.2107 0.6350
ClassAmphibia 0.0564 0.1326 0.4256 0.6704 -0.2034 0.3163
ClassActinopterygii 0.1769 0.1373 1.2881 0.1977 -0.0923 0.4460
ClassTrematoda -0.0258 0.2277 -0.1134 0.9097 -0.4721 0.4204
---
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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-15.3360 30.6720 54.6720 84.9335 58.6214
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0648 0.2546 102 no Index
Test for Residual Heterogeneity:
QE(df = 92) = 867.8593, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 23.0339, p-val = 0.0061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt -0.0121 0.2068 -0.0583 0.9535 -0.4173 0.3932
ClassReptilia 0.0258 0.2277 0.1134 0.9097 -0.4204 0.4721
ClassPisces 0.2704 0.3316 0.8152 0.4149 -0.3796 0.9204
ClassNematoda 0.4342 0.2988 1.4532 0.1462 -0.1514 1.0197
ClassMammalia 0.3997 0.2125 1.8813 0.0599 -0.0167 0.8162 .
ClassInsecta 0.2461 0.2410 1.0213 0.3071 -0.2262 0.7185
ClassAves 0.1630 0.2134 0.7641 0.4448 -0.2552 0.5812
ClassAnimalia 0.2380 0.2832 0.8402 0.4008 -0.3171 0.7931
ClassAmphibia 0.0822 0.2264 0.3633 0.7164 -0.3614 0.5259
ClassActinopterygii 0.2027 0.2292 0.8844 0.3765 -0.2465 0.6518
---
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(0.0042, 0.8853, .0001, 0.4463, 0.0537, 0.3192, 0.0587, 0.9535, 0.0503, 0.2432), method = 'fdr'),digit=3),row.names=cbind( "Animalia","Nematoda","Trematoda","Insecta","Pisces" ,"Actinopterygii","Amphibia","Mammalia","Reptilia","Aves"))
colnames(tab2)<-cbind('P-value')
tab2 P-value
Animalia 0.021
Nematoda 0.954
Trematoda 0.001
Insecta 0.558
Pisces 0.117
Actinopterygii 0.456
Amphibia 0.117
Mammalia 0.954
Reptilia 0.117
Aves 0.405Next 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 = 100; method: REML)
logLik Deviance AIC BIC AICc
-19.3573 38.7147 46.7147 57.0545 47.1448
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 58 no Study_ID
sigma^2.2 0.0708 0.2661 100 no Index
Test for Residual Heterogeneity:
QE(df = 98) = 1210.4699, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 9.6657, p-val = 0.0019
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3225 0.0434 7.4349 <.0001 0.2375 0.4076 ***
SSD_ProxyBody size -0.1844 0.0593 -3.1090 0.0019 -0.3006 -0.0681 **
---
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, 0.0006), method = 'fdr'),digit=3),row.names=cbind("Body mass","Body size"))
colnames(tab3)<-cbind('P-value')
tab3 P-value
Body mass 0.000
Body size 0.001Next 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-22.3916 44.7832 52.7832 63.2039 53.2043
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0200 0.1414 59 no Study_ID
sigma^2.2 0.0586 0.2420 102 no Index
Test for Residual Heterogeneity:
QE(df = 100) = 1267.1170, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 0.9102, p-val = 0.3401
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2422 0.0369 6.5671 <.0001 0.1699 0.3145 ***
BodySizeControlledYes -0.0886 0.0928 -0.9540 0.3401 -0.2705 0.0934
---
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.0712), method = 'fdr'),digit=3),row.names=cbind("uncontrolled","controlled"))
colnames(tab4)<-cbind('P-value')
tab4 P-value
uncontrolled 0.000
controlled 0.071Model_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 = 89; method: REML)
logLik Deviance AIC BIC AICc
-18.5867 37.1734 45.1734 55.0371 45.6612
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0221 0.1487 51 no Study_ID
sigma^2.2 0.0557 0.2359 89 no Index
Test for Residual Heterogeneity:
QE(df = 87) = 1190.6640, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 5.3018, p-val = 0.0213
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3546 0.0584 6.0739 <.0001 0.2402 0.4690 ***
SSD_SexBias_in_perc_F -0.0027 0.0012 -2.3026 0.0213 -0.0050 -0.0004 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model_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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-119.9984 239.9969 247.9969 258.4176 248.4179
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 59 no Study_ID
sigma^2.2 0.0000 0.0000 102 no Index
Test for Residual Heterogeneity:
QE(df = 100) = 17.9544, p-val = 1.0000
Test of Moderators (coefficient 2):
QM(df = 1) = 1.8587, p-val = 0.1728
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0727 0.2243 0.3243 0.7457 -0.3669 0.5124
SE_z 1.1696 0.8579 1.3633 0.1728 -0.5118 2.8510
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next 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 = 102; method: REML)
logLik Deviance AIC BIC AICc
-22.6104 45.2209 53.2209 63.6416 53.6419
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0233 0.1526 59 no Study_ID
sigma^2.2 0.0565 0.2376 102 no Index
Test for Residual Heterogeneity:
QE(df = 100) = 1162.1943, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 0.4138, p-val = 0.5200
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 5.2924 7.8708 0.6724 0.5013 -10.1341 20.7189
Year -0.0025 0.0039 -0.6433 0.5200 -0.0102 0.0051
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R version 4.3.3 (2024-02-29 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
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.35 coda_0.19-4.1
[10] RColorBrewer_1.1-3 ggplot2_3.5.0 data.table_1.15.2
[13] multcomp_1.4-25 TH.data_1.1-2 survival_3.5-8
[16] mvtnorm_1.2-4 psych_2.4.3 pwr_1.3-0
[19] MASS_7.3-60.0.1 metafor_4.4-0 numDeriv_2016.8-1.1
[22] metadat_1.2-0 Matrix_1.6-5 matrixcalc_1.0-6
[25] ape_5.7-1 workflowr_1.7.1
loaded via a namespace (and not attached):
[1] tensorA_0.36.2.1 gtable_0.3.4 xfun_0.42
[4] bslib_0.6.1 processx_3.8.4 lattice_0.22-5
[7] callr_3.7.5 mathjaxr_1.6-0 vctrs_0.6.5
[10] tools_4.3.3 ps_1.7.6 generics_0.1.3
[13] parallel_4.3.3 sandwich_3.1-0 tibble_3.2.1
[16] fansi_1.0.6 pkgconfig_2.0.3 distributional_0.4.0
[19] cubature_2.1.0 lifecycle_1.0.4 compiler_4.3.3
[22] stringr_1.5.1 git2r_0.33.0 munsell_0.5.0
[25] mnormt_2.1.1 getPass_0.2-4 codetools_0.2-19
[28] httpuv_1.6.14 htmltools_0.5.7 sass_0.4.9
[31] yaml_2.3.8 tidyr_1.3.1 later_1.3.2
[34] pillar_1.9.0 jquerylib_0.1.4 whisker_0.4.1
[37] cachem_1.0.8 nlme_3.1-164 tidyselect_1.2.1
[40] digest_0.6.35 stringi_1.8.3 purrr_1.0.2
[43] splines_4.3.3 rprojroot_2.0.4 fastmap_1.1.1
[46] grid_4.3.3 colorspace_2.1-0 cli_3.6.1
[49] magrittr_2.0.3 utf8_1.2.4 corpcor_1.6.10
[52] withr_3.0.0 scales_1.3.0 promises_1.2.1
[55] rmarkdown_2.26 httr_1.4.7 zoo_1.8-12
[58] evaluate_0.23 knitr_1.45 Rcpp_1.0.12
[61] glue_1.7.0 rstudioapi_0.16.0 jsonlite_1.8.8
[64] R6_2.5.1 fs_1.6.3