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/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 |
|---|---|---|---|---|
| html | cd21614 | LennartWinkler | 2024-08-16 | Build site. |
| Rmd | 6eb630e | LennartWinkler | 2024-08-16 | wflow_publish(republish = TRUE, 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 | 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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-20.2553 40.5106 48.5106 59.0882 48.9146
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 62 no Study_ID no
sigma^2.2 0.0645 0.2539 105 no Index no
sigma^2.3 0.0194 0.1391 10 no Class yes
Test for Heterogeneity:
Q(df = 104) = 1299.5696, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2118 0.0828 2.5582 0.0105 0.0495 0.3741 *
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-23.4061 46.8121 52.8121 60.7453 53.0521
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0211 0.1452 62 no Study_ID
sigma^2.2 0.0577 0.2403 105 no Index
Test for Heterogeneity:
Q(df = 104) = 1299.5696, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2369 0.0334 7.0898 <.0001 0.1714 0.3024 ***
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-5.2629 10.5258 22.5258 38.2756 23.4100
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0065 0.0806 62 no Study_ID no
sigma^2.2 0.0375 0.1937 105 no Index no
sigma^2.3 0.0201 0.1419 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 102) = 1050.7775, p-val < .0001
Test of Moderators (coefficients 2:3):
QM(df = 2) = 37.2830, p-val < .0001
Model Results:
estimate se zval pval
intrcpt 0.2379 0.0858 2.7709 0.0056
factor(SexSel_Episode)post-copulatory -0.3301 0.0770 -4.2853 <.0001
factor(SexSel_Episode)both 0.1535 0.0558 2.7483 0.0060
ci.lb ci.ub
intrcpt 0.0696 0.4061 **
factor(SexSel_Episode)post-copulatory -0.4810 -0.1791 ***
factor(SexSel_Episode)both 0.0440 0.2629 **
---
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.0068, 0.3664, .0001), method = 'fdr'),digit=3),row.names=cbind("Pre-copulatory","Post-copulatory","Both"))
colnames(tab1)<-cbind('P-value')
tab1 P-value
Pre-copulatory 0.010
Post-copulatory 0.366
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
4.2234 -8.4469 11.5531 37.4028 14.0819
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0151 0.1229 62 no Study_ID no
sigma^2.2 0.0224 0.1497 105 no Index no
sigma^2.3 0.0158 0.1259 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 98) = 767.4092, p-val < .0001
Test of Moderators (coefficients 2:7):
QM(df = 6) = 71.8193, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.1275 0.0957 -1.3312 0.1831
factor(SexSel_Category)Density 0.3288 0.1226 2.6824 0.0073
factor(SexSel_Category)Mating system 0.4966 0.0819 6.0658 <.0001
factor(SexSel_Category)OSR 0.7321 0.1355 5.4018 <.0001
factor(SexSel_Category)Other 0.4692 0.0962 4.8773 <.0001
factor(SexSel_Category)Premating competition 0.4273 0.0833 5.1308 <.0001
factor(SexSel_Category)Trait-based 0.1486 0.0884 1.6812 0.0927
ci.lb ci.ub
intrcpt -0.3151 0.0602
factor(SexSel_Category)Density 0.0886 0.5690 **
factor(SexSel_Category)Mating system 0.3361 0.6570 ***
factor(SexSel_Category)OSR 0.4665 0.9977 ***
factor(SexSel_Category)Other 0.2806 0.6577 ***
factor(SexSel_Category)Premating competition 0.2640 0.5905 ***
factor(SexSel_Category)Trait-based -0.0246 0.3219 .
---
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.1884, 0.8292, 0.1185, 0.0009, 0.0001, .0001, 0.0005), 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.829
Density 0.166
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 = 103; method: REML)
logLik Deviance AIC BIC AICc
-19.5015 39.0029 49.0029 62.0785 49.6345
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 61 no Study_ID no
sigma^2.2 0.0671 0.2590 103 no Index no
sigma^2.3 0.0071 0.0845 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 101) = 1236.9117, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 4.6150, p-val = 0.0317
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3224 0.0783 4.1170 <.0001 0.1689 0.4758 ***
SSD_ProxyBody size -0.1467 0.0683 -2.1482 0.0317 -0.2806 -0.0129 *
---
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.0059), 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.006Next 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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-20.1921 40.3841 50.3841 63.5578 51.0027
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 62 no Study_ID no
sigma^2.2 0.0648 0.2546 105 no Index no
sigma^2.3 0.0188 0.1370 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 103) = 1297.2607, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 0.6526, p-val = 0.4192
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2256 0.0838 2.6934 0.0071 0.0614 0.3898 **
BodySizeControlledYes -0.0692 0.0856 -0.8078 0.4192 -0.2370 0.0987
---
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.0085, 0.1492), method = 'fdr'),digit=3),row.names=cbind("uncontrolled","controlled"))
colnames(tab4)<-cbind('P-value')
tab4 P-value
uncontrolled 0.017
controlled 0.149There 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 = 91; method: REML)
logLik Deviance AIC BIC AICc
-17.0719 34.1437 44.1437 56.5869 44.8666
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0072 0.0849 53 no Study_ID no
sigma^2.2 0.0585 0.2419 91 no Index no
sigma^2.3 0.0180 0.1343 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 89) = 1207.9887, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 1.8083, p-val = 0.1787
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3162 0.1089 2.9037 0.0037 0.1028 0.5296 **
SSD_SexBias_in_perc_F -0.0018 0.0013 -1.3447 0.1787 -0.0044 0.0008
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-131.2805 262.5610 272.5610 285.7347 273.1796
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 62 no Study_ID no
sigma^2.2 0.0000 0.0000 105 no Index no
sigma^2.3 0.0000 0.0003 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 103) = 15.0179, p-val = 1.0000
Test of Moderators (coefficient 2):
QM(df = 1) = 1.7899, p-val = 0.1809
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1410 0.1956 0.7207 0.4711 -0.2424 0.5243
SE_z 0.8522 0.6369 1.3379 0.1809 -0.3962 2.1006
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-20.0318 40.0635 50.0635 63.2372 50.6821
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0005 0.0229 62 no Study_ID no
sigma^2.2 0.0639 0.2527 105 no Index no
sigma^2.3 0.0189 0.1375 10 no Class yes
Test for Residual Heterogeneity:
QE(df = 103) = 1165.7023, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 1.0473, p-val = 0.3061
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 6.7004 6.3404 1.0568 0.2906 -5.7266 19.1275
Year -0.0032 0.0032 -1.0234 0.3061 -0.0094 0.0030
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-9.6012 19.2024 29.2024 42.3272 29.8274
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0202 0.1420 62 no Study_ID
sigma^2.2 0.0381 0.1951 105 no Index
Test for Residual Heterogeneity:
QE(df = 102) = 1050.7775, p-val < .0001
Test of Moderators (coefficients 2:3):
QM(df = 2) = 32.8599, p-val < .0001
Model Results:
estimate se zval pval
intrcpt 0.2096 0.0420 4.9950 <.0001
factor(SexSel_Episode)both 0.1966 0.0593 3.3171 0.0009
factor(SexSel_Episode)post-copulatory -0.2848 0.0834 -3.4143 0.0006
ci.lb ci.ub
intrcpt 0.1274 0.2919 ***
factor(SexSel_Episode)both 0.0804 0.3128 ***
factor(SexSel_Episode)post-copulatory -0.4482 -0.1213 ***
---
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.3049, .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.305
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
2.2014 -4.4028 13.5972 36.8619 15.6426
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0259 0.1609 62 no Study_ID
sigma^2.2 0.0212 0.1456 105 no Index
Test for Residual Heterogeneity:
QE(df = 98) = 767.4092, p-val < .0001
Test of Moderators (coefficients 2:7):
QM(df = 6) = 73.3544, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.0734 0.0695 -1.0561 0.2909
factor(SexSel_Category)Other 0.4643 0.1004 4.6254 <.0001
factor(SexSel_Category)OSR 0.7100 0.1377 5.1556 <.0001
factor(SexSel_Category)Mating system 0.4771 0.0835 5.7155 <.0001
factor(SexSel_Category)Premating competition 0.4007 0.0877 4.5670 <.0001
factor(SexSel_Category)Density 0.3292 0.1262 2.6083 0.0091
factor(SexSel_Category)Trait-based 0.0878 0.0889 0.9874 0.3234
ci.lb ci.ub
intrcpt -0.2095 0.0628
factor(SexSel_Category)Other 0.2675 0.6610 ***
factor(SexSel_Category)OSR 0.4401 0.9799 ***
factor(SexSel_Category)Mating system 0.3135 0.6406 ***
factor(SexSel_Category)Premating competition 0.2287 0.5726 ***
factor(SexSel_Category)Density 0.0818 0.5766 **
factor(SexSel_Category)Trait-based -0.0865 0.2621
---
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.2909, 0.7974, 0.0156, .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.339
Trait-based 0.797
Density 0.022
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1906 0.0991 1.9235 0.0544 -0.0036 0.3849 .
ClassAmphibia -0.1204 0.1355 -0.8887 0.3742 -0.3860 0.1452
ClassAnimalia 0.0353 0.2180 0.1618 0.8714 -0.3920 0.4625
ClassAves -0.0397 0.1123 -0.3532 0.7240 -0.2598 0.1804
ClassInsecta 0.0434 0.1589 0.2735 0.7845 -0.2679 0.3548
ClassMammalia 0.2046 0.1102 1.8560 0.0635 -0.0115 0.4207 .
ClassNematoda 0.2315 0.2378 0.9734 0.3304 -0.2346 0.6976
ClassPisces 0.0677 0.2784 0.2431 0.8079 -0.4780 0.6134
ClassReptilia -0.1148 0.1317 -0.8715 0.3835 -0.3729 0.1434
ClassTrematoda -0.2027 0.2298 -0.8820 0.3778 -0.6530 0.2477
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0702 0.0924 0.7593 0.4477 -0.1110 0.2514
ClassActinopterygii 0.1204 0.1355 0.8887 0.3742 -0.1452 0.3860
ClassAnimalia 0.1557 0.2150 0.7241 0.4690 -0.2658 0.5772
ClassAves 0.0808 0.1065 0.7587 0.4480 -0.1279 0.2894
ClassInsecta 0.1639 0.1548 1.0587 0.2897 -0.1395 0.4673
ClassMammalia 0.3251 0.1043 3.1166 0.0018 0.1206 0.5295 **
ClassNematoda 0.3519 0.2351 1.4968 0.1345 -0.1089 0.8128
ClassPisces 0.1881 0.2761 0.6813 0.4957 -0.3531 0.7293
ClassReptilia 0.0056 0.1268 0.0445 0.9645 -0.2428 0.2541
ClassTrematoda -0.0822 0.2270 -0.3623 0.7171 -0.5271 0.3627
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2259 0.1942 1.1635 0.2446 -0.1546 0.6064
ClassAmphibia -0.1557 0.2150 -0.7241 0.4690 -0.5772 0.2658
ClassActinopterygii -0.0353 0.2180 -0.1618 0.8714 -0.4625 0.3920
ClassAves -0.0749 0.2012 -0.3724 0.7096 -0.4693 0.3194
ClassInsecta 0.0082 0.2305 0.0355 0.9717 -0.4435 0.4599
ClassMammalia 0.1693 0.2001 0.8464 0.3973 -0.2228 0.5615
ClassNematoda 0.1962 0.2906 0.6753 0.4995 -0.3733 0.7657
ClassPisces 0.0324 0.3246 0.0998 0.9205 -0.6039 0.6687
ClassReptilia -0.1501 0.2127 -0.7057 0.4804 -0.5669 0.2667
ClassTrematoda -0.2379 0.2840 -0.8377 0.4022 -0.7946 0.3187
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1510 0.0528 2.8577 0.0043 0.0474 0.2545 **
ClassAnimalia 0.0749 0.2012 0.3724 0.7096 -0.3194 0.4693
ClassAmphibia -0.0808 0.1065 -0.7587 0.4480 -0.2894 0.1279
ClassActinopterygii 0.0397 0.1123 0.3532 0.7240 -0.1804 0.2598
ClassInsecta 0.0831 0.1349 0.6159 0.5380 -0.1814 0.3476
ClassMammalia 0.2443 0.0716 3.4126 0.0006 0.1040 0.3846 ***
ClassNematoda 0.2712 0.2226 1.2184 0.2231 -0.1650 0.7074
ClassPisces 0.1073 0.2655 0.4043 0.6860 -0.4130 0.6277
ClassReptilia -0.0751 0.1016 -0.7396 0.4595 -0.2742 0.1240
ClassTrematoda -0.1630 0.2139 -0.7620 0.4461 -0.5823 0.2563
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2341 0.1242 1.8851 0.0594 -0.0093 0.4774 .
ClassAves -0.0831 0.1349 -0.6159 0.5380 -0.3476 0.1814
ClassAnimalia -0.0082 0.2305 -0.0355 0.9717 -0.4599 0.4435
ClassAmphibia -0.1639 0.1548 -1.0587 0.2897 -0.4673 0.1395
ClassActinopterygii -0.0434 0.1589 -0.2735 0.7845 -0.3548 0.2679
ClassMammalia 0.1612 0.1332 1.2097 0.2264 -0.1000 0.4223
ClassNematoda 0.1880 0.2493 0.7543 0.4507 -0.3006 0.6767
ClassPisces 0.0242 0.2883 0.0841 0.9330 -0.5408 0.5893
ClassReptilia -0.1582 0.1515 -1.0446 0.2962 -0.4551 0.1387
ClassTrematoda -0.2461 0.2417 -1.0185 0.3085 -0.7198 0.2275
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3952 0.0483 8.1823 <.0001 0.3006 0.4899 ***
ClassInsecta -0.1612 0.1332 -1.2097 0.2264 -0.4223 0.1000
ClassAves -0.2443 0.0716 -3.4126 0.0006 -0.3846 -0.1040 ***
ClassAnimalia -0.1693 0.2001 -0.8464 0.3973 -0.5615 0.2228
ClassAmphibia -0.3251 0.1043 -3.1166 0.0018 -0.5295 -0.1206 **
ClassActinopterygii -0.2046 0.1102 -1.8560 0.0635 -0.4207 0.0115 .
ClassNematoda 0.0269 0.2215 0.1213 0.9034 -0.4073 0.4611
ClassPisces -0.1369 0.2646 -0.5175 0.6048 -0.6556 0.3817
ClassReptilia -0.3194 0.0993 -3.2164 0.0013 -0.5140 -0.1248 **
ClassTrematoda -0.4073 0.2129 -1.9133 0.0557 -0.8245 0.0099 .
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.4221 0.2162 1.9525 0.0509 -0.0016 0.8459 .
ClassMammalia -0.0269 0.2215 -0.1213 0.9034 -0.4611 0.4073
ClassInsecta -0.1880 0.2493 -0.7543 0.4507 -0.6767 0.3006
ClassAves -0.2712 0.2226 -1.2184 0.2231 -0.7074 0.1650
ClassAnimalia -0.1962 0.2906 -0.6753 0.4995 -0.7657 0.3733
ClassAmphibia -0.3519 0.2351 -1.4968 0.1345 -0.8128 0.1089
ClassActinopterygii -0.2315 0.2378 -0.9734 0.3304 -0.6976 0.2346
ClassPisces -0.1638 0.3383 -0.4843 0.6282 -0.8269 0.4992
ClassReptilia -0.3463 0.2330 -1.4865 0.1372 -0.8029 0.1103
ClassTrematoda -0.4342 0.2995 -1.4495 0.1472 -1.0213 0.1529
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2583 0.2602 0.9927 0.3208 -0.2517 0.7683
ClassNematoda 0.1638 0.3383 0.4843 0.6282 -0.4992 0.8269
ClassMammalia 0.1369 0.2646 0.5175 0.6048 -0.3817 0.6556
ClassInsecta -0.0242 0.2883 -0.0841 0.9330 -0.5893 0.5408
ClassAves -0.1073 0.2655 -0.4043 0.6860 -0.6277 0.4130
ClassAnimalia -0.0324 0.3246 -0.0998 0.9205 -0.6687 0.6039
ClassAmphibia -0.1881 0.2761 -0.6813 0.4957 -0.7293 0.3531
ClassActinopterygii -0.0677 0.2784 -0.2431 0.8079 -0.6134 0.4780
ClassReptilia -0.1825 0.2743 -0.6653 0.5059 -0.7200 0.3551
ClassTrematoda -0.2704 0.3327 -0.8126 0.4164 -0.9224 0.3817
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.0758 0.0868 0.8740 0.3821 -0.0942 0.2459
ClassPisces 0.1825 0.2743 0.6653 0.5059 -0.3551 0.7200
ClassNematoda 0.3463 0.2330 1.4865 0.1372 -0.1103 0.8029
ClassMammalia 0.3194 0.0993 3.2164 0.0013 0.1248 0.5140 **
ClassInsecta 0.1582 0.1515 1.0446 0.2962 -0.1387 0.4551
ClassAves 0.0751 0.1016 0.7396 0.4595 -0.1240 0.2742
ClassAnimalia 0.1501 0.2127 0.7057 0.4804 -0.2667 0.5669
ClassAmphibia -0.0056 0.1268 -0.0445 0.9645 -0.2541 0.2428
ClassActinopterygii 0.1148 0.1317 0.8715 0.3835 -0.1434 0.3729
ClassTrematoda -0.0879 0.2247 -0.3911 0.6958 -0.5284 0.3526
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-16.3766 32.7533 56.7533 87.3998 60.5581
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0653 0.2555 105 no Index
Test for Residual Heterogeneity:
QE(df = 95) = 881.3207, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 22.4089, p-val = 0.0077
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt -0.0121 0.2073 -0.0581 0.9536 -0.4184 0.3943
ClassReptilia 0.0879 0.2247 0.3911 0.6958 -0.3526 0.5284
ClassPisces 0.2704 0.3327 0.8126 0.4164 -0.3817 0.9224
ClassNematoda 0.4342 0.2995 1.4495 0.1472 -0.1529 1.0213
ClassMammalia 0.4073 0.2129 1.9133 0.0557 -0.0099 0.8245 .
ClassInsecta 0.2461 0.2417 1.0185 0.3085 -0.2275 0.7198
ClassAves 0.1630 0.2139 0.7620 0.4461 -0.2563 0.5823
ClassAnimalia 0.2379 0.2840 0.8377 0.4022 -0.3187 0.7946
ClassAmphibia 0.0822 0.2270 0.3623 0.7171 -0.3627 0.5271
ClassActinopterygii 0.2027 0.2298 0.8820 0.3778 -0.2477 0.6530
---
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.0043, 0.3821, .0001, 0.4477, 0.0544, 0.3208, 0.0594, 0.9536, 0.0509, 0.2446), 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.478
Trematoda 0.001
Insecta 0.497
Pisces 0.119
Actinopterygii 0.458
Amphibia 0.119
Mammalia 0.954
Reptilia 0.119
Aves 0.408Next 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 = 103; method: REML)
logLik Deviance AIC BIC AICc
-20.0433 40.0866 48.0866 58.5471 48.5033
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 61 no Study_ID
sigma^2.2 0.0706 0.2658 103 no Index
Test for Residual Heterogeneity:
QE(df = 101) = 1236.9117, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 9.8855, p-val = 0.0017
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3297 0.0428 7.6972 <.0001 0.2457 0.4136 ***
SSD_ProxyBody size -0.1836 0.0584 -3.1441 0.0017 -0.2980 -0.0691 **
---
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.0002), method = 'fdr'),digit=3),row.names=cbind("Body mass","Body size"))
colnames(tab3)<-cbind('P-value')
tab3 P-value
Body mass 0
Body size 0Next 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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-23.0332 46.0664 54.0664 64.6054 54.4746
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0206 0.1436 62 no Study_ID
sigma^2.2 0.0578 0.2404 105 no Index
Test for Residual Heterogeneity:
QE(df = 103) = 1297.2607, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 1.1120, p-val = 0.2916
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.2517 0.0362 6.9573 <.0001 0.1808 0.3226 ***
BodySizeControlledYes -0.0977 0.0927 -1.0545 0.2916 -0.2794 0.0839
---
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 = 91; method: REML)
logLik Deviance AIC BIC AICc
-18.4446 36.8892 44.8892 54.8438 45.3654
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0217 0.1474 53 no Study_ID
sigma^2.2 0.0548 0.2341 91 no Index
Test for Residual Heterogeneity:
QE(df = 89) = 1207.9887, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 6.1191, p-val = 0.0134
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3627 0.0564 6.4284 <.0001 0.2521 0.4733 ***
SSD_SexBias_in_perc_F -0.0029 0.0012 -2.4737 0.0134 -0.0051 -0.0006 *
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-131.2805 262.5610 270.5610 281.0999 270.9692
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 62 no Study_ID
sigma^2.2 0.0000 0.0000 105 no Index
Test for Residual Heterogeneity:
QE(df = 103) = 15.0179, p-val = 1.0000
Test of Moderators (coefficient 2):
QM(df = 1) = 1.7899, p-val = 0.1809
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.1410 0.1956 0.7207 0.4711 -0.2424 0.5243
SE_z 0.8522 0.6369 1.3379 0.1809 -0.3962 2.1006
---
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 = 105; method: REML)
logLik Deviance AIC BIC AICc
-22.7935 45.5871 53.5871 64.1260 53.9953
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0243 0.1560 62 no Study_ID
sigma^2.2 0.0544 0.2333 105 no Index
Test for Residual Heterogeneity:
QE(df = 103) = 1165.7023, p-val < .0001
Test of Moderators (coefficient 2):
QM(df = 1) = 1.5764, p-val = 0.2093
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 9.3614 7.2665 1.2883 0.1976 -4.8806 23.6034
Year -0.0045 0.0036 -1.2555 0.2093 -0.0116 0.0025
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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] tensorA_0.36.2.1 gtable_0.3.5 xfun_0.43
[4] bslib_0.7.0 processx_3.8.4 lattice_0.22-6
[7] callr_3.7.6 mathjaxr_1.6-0 vctrs_0.6.5
[10] tools_4.4.0 ps_1.7.6 generics_0.1.3
[13] parallel_4.4.0 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.4.0
[22] stringr_1.5.1 git2r_0.33.0 munsell_0.5.1
[25] mnormt_2.1.1 getPass_0.2-4 codetools_0.2-20
[28] httpuv_1.6.15 htmltools_0.5.8.1 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.4 purrr_1.0.2
[43] splines_4.4.0 rprojroot_2.0.4 fastmap_1.1.1
[46] grid_4.4.0 colorspace_2.1-0 cli_3.6.2
[49] magrittr_2.0.3 utf8_1.2.4 corpcor_1.6.10
[52] withr_3.0.0 scales_1.3.0 promises_1.3.0
[55] rmarkdown_2.26 httr_1.4.7 zoo_1.8-12
[58] evaluate_0.23 knitr_1.46 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.4