# ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # # # Supplementary R code for: # # THE COSTS AND BENEFITS OF PATERNAL CARE IN FISH: A META-ANALYSIS # # Authors: Rebecca L. Goldberg, Philip A. Downing, Ashleigh S. Griffin and Jonathan P. Green # Journal: Proceedings of the Royal Society B: Biological Sciences # DOI: 10.1098/rspb.2020.1759 # contact: beccalgoldberg@gmail.com # # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # # Packages library(ape) library(phytools) library(MCMCglmm) library(metafor) library(dplyr) library(ggplot2) ### DATA MANIPULATION ### ## All effect sizes. Read in 'GoldbergetalEffectSizeData.csv (used in PART 1-5, 7 & 8) ## # Create separate dataframe for each analysis # # ZrCond # condData <- zrData[!is.na(zrData$ZrCond),] # 43 effect sizes from 25 studies on 19 species condData$species <- factor(condData$species) condData$ID <- condData$species # need to use to control for repeated measures condData$obs <- paste(condData$species, 1:length(condData$species), sep = "") # need a unique identifier to estimate residual variance in metafor # ZrOff # offData <- zrData[!is.na(zrData$ZrOff),] # 36 effect sizes from 18 studies on 16 species offData$species <- factor(offData$species) offData$ID <- offData$species offData$obs <- paste(offData$species, 1:length(offData$species), sep = "") # ZrPref1 # prefData1 <- zrData[!is.na(zrData$ZrPref1),] # 57 effect sizes from 34 studies on 23 species prefData1$species <- factor(prefData1$species) prefData1$ID <- prefData1$species prefData1$obs <- paste(prefData1$species, 1:length(prefData1$species), sep = "") # ZrPref2 # prefData2 <- zrData[!is.na(zrData$ZrPref2),] # 30 effect sizes from 15 studies on 10 species prefData2$species <- factor(prefData2$species) prefData2$ID <- prefData2$species prefData2$obs <- paste(prefData2$species, 1:length(prefData2$species), sep = "") # ZrBrood # broodData <- zrData[!is.na(zrData$ZrBrood),] # 87 effect sizes from 50 studies on 25 species broodData$species <- factor(broodData$species) broodData$ID <- broodData$species broodData$obs <- paste(broodData$species, 1:length(broodData$species), sep = "") ## study average effect sizes. Read in GoldbergetalStudyLevelData.csv (used in PART 6 and PART 7) ## # separate dataframe for each analysis and collapse to species averages # # ZrCond # condStudy <- studyData[!is.na(studyData$ZrCond),] condSpecies <- group_by(condStudy, species) condSpecies <- summarize(condSpecies, common=common[1], ZrCond=weighted.mean(x=ZrCond, w=Ncond), Nstudy=length(species), Ncond=sum(Ncond)) condSpecies$VarCond <- 1/(condSpecies$Ncond-3) condSpecies$obs <- paste(condSpecies$species, 1:length(condSpecies$species), sep = "") # need a unique identifier to estimate residual variance in metafor condSpecies$lciCond <- condSpecies$ZrCond - (sqrt(condSpecies$VarCond)*1.96) condSpecies$uciCond <- condSpecies$ZrCond + (sqrt(condSpecies$VarCond)*1.96) # ZrOff # offStudy <- studyData[!is.na(studyData$ZrOff),] offSpecies<-group_by(offStudy, species) offSpecies<-summarize(offSpecies, common=common[1], ZrOff=weighted.mean(x=ZrOff,w=Noff), Noff=sum(Noff), Nstudy=length(species)) offSpecies$VarOff <- 1/(offSpecies$Noff-3) offSpecies$obs <- paste(offSpecies$species, 1:length(offSpecies$species), sep = "") offSpecies$lciOff <- offSpecies$ZrOff - (sqrt(offSpecies$VarOff)*1.96) offSpecies$uciOff <- offSpecies$ZrOff + (sqrt(offSpecies$VarOff)*1.96) # ZrPref1 # prefStudy1 <- studyData[!is.na(studyData$ZrPref1),] prefSpecies1<-group_by(prefStudy1, species) prefSpecies1<-summarize(prefSpecies1, common=common[1], ZrPref1=weighted.mean(x=ZrPref1,w=Npref1), Npref1=sum(Npref1), Nstudy=length(species)) prefSpecies1$VarPref1 <- 1/(prefSpecies1$Npref1-3) prefSpecies1$obs <- paste(prefSpecies1$species, 1:length(prefSpecies1$species), sep = "") prefSpecies1$lciPref1 <- prefSpecies1$ZrPref1 - (sqrt(prefSpecies1$VarPref1)*1.96) prefSpecies1$uciPref1 <- prefSpecies1$ZrPref1 + (sqrt(prefSpecies1$VarPref1)*1.96) # ZrPref2 # prefStudy2 <- studyData[!is.na(studyData$ZrPref2),] prefSpecies2<-group_by(prefStudy2, species) prefSpecies2<-summarize(prefSpecies2, common=common[1], ZrPref2=weighted.mean(x=ZrPref2,w=Npref2), Npref2=sum(Npref2), Nstudy=length(species)) prefSpecies2$VarPref2 <- 1/(prefSpecies2$Npref2-3) prefSpecies2$obs <- paste(prefSpecies2$species, 1:length(prefSpecies2$species), sep = "") prefSpecies2$lciPref2 <- prefSpecies2$ZrPref2 - (sqrt(prefSpecies2$VarPref2)*1.96) prefSpecies2$uciPref2 <- prefSpecies2$ZrPref2 + (sqrt(prefSpecies2$VarPref2)*1.96) # ZrBrood # broodStudy <- studyData[!is.na(studyData$ZrBrood),] broodSpecies<-group_by(broodStudy, species) broodSpecies<-summarize(broodSpecies, common=common[1], ZrBrood=weighted.mean(x=ZrBrood,w=Nbrood), Nbrood=sum(Nbrood), Nstudy=length(species)) broodSpecies$VarBrood <- 1/(broodSpecies$Nbrood-3) broodSpecies$obs <- paste(broodSpecies$species, 1:length(broodSpecies$species), sep = "") broodSpecies$lciBrood <- broodSpecies$ZrBrood - (sqrt(broodSpecies$VarBrood)*1.96) broodSpecies$uciBrood <- broodSpecies$ZrBrood + (sqrt(broodSpecies$VarBrood)*1.96) ## phylogenetic trees ## # Load in tree created in PART 8 # is.ultrametric(tree) zrData$species[which((zrData$species %in% tree$tip.label) == FALSE)] tree$tip.label[which((tree$tip.label %in% zrData$species) == FALSE)] length(tree$tip.label) # 48 species # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### WITHIN-STUDY VARIANCE ### ## need this in order to calculate I^2 values ## see Nakagawa & Santos (2012) Methodological issues and advances in biological meta-analysis. Evol. Ecol 26: 1253- 1274 ## ZrCond within-study variance wjCond <- 1/condData$VarCond # inverse measurement error / sampling error variance kCond <- length(condData$ZrCond) # number of species withinCond <- sum(wjCond * (kCond -1)) / (sum(wjCond)^2 - sum(wjCond^2)) # 0.026 ## ZrOff within-study variance wjOff <- 1/offData$VarOff # inverse measurement error / sampling error variance kOff <- length(offData$ZrOff) # number of species withinOff <- sum(wjOff * (kOff -1)) / (sum(wjOff)^2 - sum(wjOff^2)) # 0.022 ## ZrPref1 within-study variance wjPref1 <- 1/prefData1$VarPref1 # inverse measurement error / sampling error variance kPref1 <- length(prefData1$ZrPref1) # number of species withinPref1 <- sum(wjPref1 * (kPref1 -1)) / (sum(wjPref1)^2 - sum(wjPref1^2)) # 0.029 ## ZrPref2 within-study variance wjPref2 <- 1/prefData2$VarPref2 # inverse measurement error / sampling error variance kPref2 <- length(prefData2$ZrPref2) # number of species withinPref2 <- sum(wjPref2 * (kPref2 -1)) / (sum(wjPref2)^2 - sum(wjPref2^2)) # 0.043 ## ZrBrood within-study variance wjBrood <- 1/broodData$VarBrood # inverse measurement error / sampling error variance kBrood <- length(broodData$ZrBrood) # number of species withinBrood <- sum(wjBrood * (kBrood -1)) / (sum(wjBrood)^2 - sum(wjBrood^2)) # 0.019 ### OTHER VARIANCE PARAMETERS FROM MCMCglmm ### # between <- posterior.mode(model$VCV[,X]) i.e. the between-study variance (estimated from the data) # repeated <- posterior.mode(model$VCV[,X]) i.e. the repreated measures variance (estimated from the data) # phylo <- posterior.mode(model$VCV[,X]) i.e. phylogenetic variance (estimated from the data) # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 1.1: COSTS OF MALE CARE (Zr Condition) WITHOUT PHYLOGENY ### ## A. Mean effect size ## # metafor # summary(rma(ZrCond ~ 1, vi=VarCond, data=condData)) # no random effects: B = -0.24, lwr = -0.36, upr = -0.13 trimfill(rma(ZrCond ~ 1, vi=VarCond, data=condData)) # 9 studies missing from the right-hand side funnel(trimfill(rma(ZrCond ~ 1, vi=VarCond, data=condData))) # funnel plot with trim and fill condData$invSE <- 1/sqrt(condData$VarCond) summary(lm(scale(ZrCond) ~ invSE, data=condData)) # Egger's test: intercept = -0.40 (p = 0.35), slope = 0.07 (p = 0.32) meanModel1.1a <- rma.mv(ZrCond ~ 1, V=VarCond, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyCond), data = condData) summary(meanModel1.1a) # B = -0.26, lwr = -0.39, upr = -0.13 (Q = 214.4, p < 0.001) forest(meanModel1.1a) # within-study variance: withinCond = 0.026 I2.between <- 0.0526 / (0.0643 + 0.0000 + 0.0526 + withinCond) * 100 # 37 % I2.repeated <- 0.0000 / (0.0643 + 0.0000 + 0.0526 + withinCond) * 100 # 0 % # MCMCglmm # priorA <- list(R=list(V=1,nu=0.02), G=list(G1=list(V=1,nu=0.02), G2=list(V=1,nu=0.02))) MEV <- condData$VarCond meanModel1.1b <- MCMCglmm(ZrCond ~ 1, random = ~ID+StudyCond, data=condData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorA, verbose = FALSE) # run the above three times meanModel1.1b.1 <- meanModel1.1b meanModel1.1b.2 <- meanModel1.1b meanModel1.1b.3 <- meanModel1.1b # chain convergence hist(meanModel1.1b.1$Liab) plot(meanModel1.1b.1$VCV) # all parameters well mixed plot(meanModel1.1b.1$Sol) # intercept estimate well mixed autocorr(meanModel1.1b.1$VCV) # correlation between successive samples < 0.1 for all parameters autocorr(meanModel1.1b.1$Sol) # correlation between successive samples < 0.1 for intercept meanModel1bSols <- mcmc.list(list(meanModel1.1b.1$Sol, meanModel1.1b.2$Sol, meanModel1.1b.3$Sol)) plot(meanModel1bSols) gelman.diag(meanModel1bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel1.1b.1$VCV) # units passed halfwidth heidel.diag(meanModel1.1b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel1.1b.1) # effective sample = 1000 posterior.mode(meanModel1.1b.1$Sol) # Zr = -0.23 HPDinterval(meanModel1.1b.1$Sol) # lwr = -0.45, upr = -0.09 (exp(2*-0.23)-1) / (exp(2*-0.23)+1) # r = -0.23 (small effect size) between <- posterior.mode(meanModel1.1b.1$VCV[,2]) repeated <- posterior.mode(meanModel1.1b.1$VCV[,1]) residual <- posterior.mode(meanModel1.1b.1$VCV[,4]) I2.between <- between / (between + repeated + residual + withinCond) * 100 # 12 % I2.repeated <- repeated / (between + repeated + residual + withinCond) * 100 # 8 % ### PART 1.2: COSTS OF MALE CARE (Zr Condition) WITH PHYLOGENY ### ## A. Mean effect size ## # metafor #random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyCond) dropTip <- tree$tip.label[which(is.na(match(tree$tip.label, condData$species)))] condTree <- drop.tip(tree, dropTip, trim.internal=T) condTree <- makeNodeLabel(condTree, method = "number") condTree <- chronoMPL(condTree) condCor <- vcv.phylo(condTree, cor=TRUE) meanModel1.2a <- rma.mv(ZrCond ~ 1, V=VarCond, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyCond, ~ 1 | species), R = list(species=condCor), data=condData) summary(meanModel1.2a) # B = -0.26, lwr = -0.39, upr = -0.13, (Q = 216.19, p < 0.001) (exp(2*-0.26)-1) / (exp(2*-0.26)+1) # r = -0.25 (small effect size) I2.between <- 0.0526 / (0.0643 + 0.0000 + 0.0526 + 0.0000 + withinCond) * 100 # 37 % I2.repeated <- 0.0000 / (0.0643 + 0.0000 + 0.0526 + 0.0000 + withinCond) * 100 # 0 % I2.phylo <- 0.0000 / (0.0643 + 0.0000 + 0.0526 + 0.0000 + withinCond) * 100 # 0 % # MCMCglmm # condData$animal <- condData$species priorB <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1,nu=0.02), G2=list(V=1,nu=0.02), G3=list(V=1,nu=0.02))) MEV <- condData$VarCond INtree <- inverseA(condTree, nodes="TIPS") meanModel1.2b <- MCMCglmm(ZrCond ~ 1, random = ~ID+StudyCond+animal, data=condData, pl=TRUE, slice=TRUE, nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE, pr=TRUE) meanModel1.2b.1 <- meanModel1.2b meanModel1.2b.2 <- meanModel1.2b meanModel1.2b.3 <- meanModel1.2b # chain convergence hist(meanModel1.2b.1$Liab) plot(meanModel1.2b.1$VCV) # all parameters well mixed plot(meanModel1.2b.1$Sol) # intercept estimate well mixed autocorr(meanModel1.2b.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(meanModel1.2b.1$Sol) # correlation between successive samples < 0.1 for all components meanModel2bSols <- mcmc.list(list(meanModel1.2b.1$Sol, meanModel1.2b.2$Sol, meanModel1.2b.3$Sol)) plot(meanModel2bSols) gelman.diag(meanModel2bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel1.2b.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel1.2b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel1.2b.1) posterior.mode(meanModel1.2b.1$Sol) # intercept = -0.28 HPDinterval(meanModel1.2b.1$Sol) # intercept = -0.62 to 0.08 (exp(2*-0.29)-1) / (exp(2*-0.29)+1) # r = -0.28 (small effect size) # I^2 values between <- posterior.mode(meanModel1.2b.1$VCV[,2]) repeated <- posterior.mode(meanModel1.2b.1$VCV[,1]) phylo <- posterior.mode(meanModel1.2b.1$VCV[,3]) residual <- posterior.mode(meanModel1.2b.1$VCV[,5]) I2.between <- between / (between + repeated + phylo + residual + withinCond) * 100 # 11 % I2.repeated <- repeated / (between + repeated + phylo + residual + withinCond) * 100 # 8 % I2.phylo <- phylo/ (between + repeated + phylo + residual + withinCond) * 100 # 10 % ## B. Moderators ## # Differences between lab and field studies # ggplot(condData, aes(x=labfieldCon, y=ZrCond))+ geom_point(aes(fill=labfieldCon), position = position_jitterdodge())+ theme_bw() MEV <- condData$VarCond INtree <- inverseA(tree, nodes="TIPS") meanModel1.2c <- MCMCglmm(ZrCond ~ labfieldCon-1, random = ~ID+StudyCond+animal, data=condData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel1.2c) posterior.mode(meanModel1.2c$Sol) HPDinterval(meanModel1.2c$Sol) # Differences between the categories of condition # meanModel1.2d <- MCMCglmm(ZrCond ~ condcatCond-1, random = ~ID+StudyCond+animal, data=condData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel1.2d) posterior.mode(meanModel1.2d$Sol) # -0.24, -0.14, -0.40, -0.25, -0.28 (all pretty close, so unlikely to bias) HPDinterval(meanModel1.2d$Sol) table(meanModel1.2d$Sol[, 1] > meanModel1.2d$Sol[, 2]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.41, TRUE = 0.59 table(meanModel1.2d$Sol[, 1] > meanModel1.2d$Sol[, 3]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.231, TRUE = 0.769 table(meanModel1.2d$Sol[, 1] > meanModel1.2d$Sol[, 4]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.39, TRUE = 0.61 table(meanModel1.2d$Sol[, 1] > meanModel1.2d$Sol[, 5]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.47, TRUE = 0.53 table(meanModel1.2d$Sol[, 1] > meanModel1.2d$Sol[, 6]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.124, TRUE = 0.876 table(meanModel1.2d$Sol[, 2] > meanModel1.2d$Sol[, 3]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.433, TRUE = 0.567 table(meanModel1.2d$Sol[, 2] > meanModel1.2d$Sol[, 4]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.483, TRUE = 0.517 table(meanModel1.2d$Sol[, 2] > meanModel1.2d$Sol[, 5]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.546, TRUE = 0.545 table(meanModel1.2d$Sol[, 2] > meanModel1.2d$Sol[, 6]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.343, TRUE = 0.657 table(meanModel1.2d$Sol[, 3] > meanModel1.2d$Sol[, 4]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.56, TRUE = 0.44 table(meanModel1.2d$Sol[, 3] > meanModel1.2d$Sol[, 5]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.593, TRUE = 0.407 table(meanModel1.2d$Sol[, 3] > meanModel1.2d$Sol[, 6]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.359, TRUE = 0.641 table(meanModel1.2d$Sol[, 4] > meanModel1.2d$Sol[, 5]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.557, TRUE = 0.443 table(meanModel1.2d$Sol[, 4] > meanModel1.2d$Sol[, 6]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.325, TRUE = 0.675 table(meanModel1.2d$Sol[, 5] > meanModel1.2d$Sol[, 6]) / length(meanModel1.2d$Sol[, 1]) # FALSE = 0.342, TRUE = 0.658 # Differences between the amount and the probability of care # MEV <- condData$VarCond INtree <- inverseA(condTree, nodes="TIPS") meanModel1.2e <- MCMCglmm(ZrCond ~ amountprobCond-1, random = ~ID+StudyCond+animal, data=condData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel1.2e) posterior.mode(meanModel1.2e$Sol) # -0.29, -0.22 HPDinterval(meanModel1.2e$Sol) # Differences between the categories of care # MEV <- condData$VarCond INtree <- inverseA(condTree, nodes="TIPS") meanModel1.2f <- MCMCglmm(ZrCond ~ carecatCond-1, random = ~ID+StudyCond+animal, data=condData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel1.2f) posterior.mode(meanModel1.2f$Sol) # -0.12, -0.37, -0.01, -0.32, -0.07 (all pretty close, so unlikely to bias) table(meanModel1.2f$Sol[, 1] > meanModel1.2f$Sol[, 2]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.248, TRUE = 0.752 table(meanModel1.2f$Sol[, 1] > meanModel1.2f$Sol[, 3]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.741, TRUE = 0.259 table(meanModel1.2f$Sol[, 1] > meanModel1.2f$Sol[, 4]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.459, TRUE = 0.541 table(meanModel1.2f$Sol[, 1] > meanModel1.2f$Sol[, 5]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.573, TRUE = 0.427 table(meanModel1.2f$Sol[, 2] > meanModel1.2f$Sol[, 3]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.863, TRUE = 0.137 table(meanModel1.2f$Sol[, 2] > meanModel1.2f$Sol[, 4]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.674, TRUE = 0.326 table(meanModel1.2f$Sol[, 2] > meanModel1.2f$Sol[, 5]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.661, TRUE = 0.339 table(meanModel1.2f$Sol[, 3] > meanModel1.2f$Sol[, 4]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.272, TRUE = 0.728 table(meanModel1.2f$Sol[, 3] > meanModel1.2f$Sol[, 5]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.406, TRUE = 0.594 table(meanModel1.2f$Sol[, 4] > meanModel1.2f$Sol[, 5]) / length(meanModel1.2f$Sol[, 1]) # FALSE = 0.571, TRUE = 0.429 # Differences between observational and experimental studies # MEV <- condData$VarCond INtree <- inverseA(tree, nodes="TIPS") meanModel1.2g <- MCMCglmm(ZrCond ~ obsexptCond, random = ~StudyCond+common+animal, data=condData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel1.2g) posterior.mode(meanModel1.2g$Sol) # 0.32, 0.18 HPDinterval(meanModel1.2g$Sol) # Differences between brooding environments # MEV <- condData$VarCond INtree <- inverseA(condTree, nodes="TIPS") meanModel1.2h <- MCMCglmm(ZrCond ~ brooding_environment-1, random = ~ID+StudyCond+animal, data=condData, pl=TRUE, slice=TRUE, nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel1.2h) posterior.mode(meanModel1.2h$Sol) # -0.27, -0.08, -0.30, 0.03, -0.42, -0.45 HPDinterval(meanModel1.2h$Sol) table(meanModel1.2h$Sol[, 1] > meanModel1.2h$Sol[, 2]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.697, TRUE = 0.303 table(meanModel1.2h$Sol[, 1] > meanModel1.2h$Sol[, 3]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.562, TRUE = 0.438 table(meanModel1.2h$Sol[, 1] > meanModel1.2h$Sol[, 4]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.709, TRUE = 0.291 table(meanModel1.2h$Sol[, 1] > meanModel1.2h$Sol[, 5]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.353, TRUE = 0.647 table(meanModel1.2h$Sol[, 1] > meanModel1.2h$Sol[, 6]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.534, TRUE = 0.466 table(meanModel1.2h$Sol[, 2] > meanModel1.2h$Sol[, 3]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.258, TRUE = 0.742 table(meanModel1.2h$Sol[, 2] > meanModel1.2h$Sol[, 4]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.551, TRUE = 0.449 table(meanModel1.2h$Sol[, 2] > meanModel1.2h$Sol[, 5]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.103, TRUE = 0.897 table(meanModel1.2h$Sol[, 2] > meanModel1.2h$Sol[, 6]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.297, TRUE = 0.703 table(meanModel1.2h$Sol[, 3] > meanModel1.2h$Sol[, 4]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.707, TRUE = 0.293 table(meanModel1.2h$Sol[, 3] > meanModel1.2h$Sol[, 5]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.189, TRUE = 0.811 table(meanModel1.2h$Sol[, 3] > meanModel1.2h$Sol[, 6]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.442, TRUE = 0.558 table(meanModel1.2h$Sol[, 4] > meanModel1.2h$Sol[, 5]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.11, TRUE = 0.89 table(meanModel1.2h$Sol[, 4] > meanModel1.2h$Sol[, 6]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.302, TRUE = 0.698 table(meanModel1.2h$Sol[, 5] > meanModel1.2h$Sol[, 6]) / length(meanModel1.2h$Sol[, 1]) # FALSE = 0.698, TRUE = 0.302 # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 2.1: BENEFITS OF MALE CARE (Zr Offspring) WITHOUT PHYLOGENY ### ## A. Mean effect size ## # metafor # summary(rma(ZrOff ~ 1, vi=VarOff, data=offData)) # no random effects: B = 0.25, lwr = 0.07, upr = 0.44 trimfill(rma(ZrOff ~ 1, vi=VarOff, data=offData)) # 10 studies missing from the left-hand side funnel(trimfill(rma(ZrOff ~ 1, vi=VarOff, data=offData))) # funnel plot with trim and fill offData$invSE <- 1/sqrt(offData$VarOff) summary(lm(scale(ZrOff) ~ invSE, data=offData)) # Egger's test: intercept = 0.14 (p = 0.69), slope = -0.02 (p = 0.64) meanModel2.1a <- rma.mv(ZrOff ~ 1, V=VarOff, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyOff), data = offData) summary(meanModel2.1a) # B = 0.25, lwr = 0.03, upr = 0.46 (Q = 256.2, p < 0.001) forest(meanModel2.1a) # within-study variance: withinOff = 0.022 I2.between <- 0.0332 / (0.2243 + 0.0204 + 0.0332 + withinOff) * 100 # 11 % I2.repeated <- 0.0204 / (0.2243 + 0.0204 + 0.0332 + withinOff) * 100 # 7 % # MCMCglmm # MEV <- offData$VarOff meanModel2.1b <- MCMCglmm(ZrOff ~ 1, random = ~ID+StudyOff, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorA, verbose = FALSE) # run the above three times meanModel2.1b.1 <- meanModel2.1b meanModel2.1b.2 <- meanModel2.1b meanModel2.1b.3 <- meanModel2.1b # chain convergence hist(meanModel2.1b.1$Liab) plot(meanModel2.1b.1$VCV) # all parameters well mixed plot(meanModel2.1b.1$Sol) # intercept estimate well mixed autocorr(meanModel2.1b.1$VCV) # correlation between successive samples < 0.1 for all parameters autocorr(meanModel2.1b.1$Sol) # correlation between successive samples < 0.1 for intercept meanModel2bSols <- mcmc.list(list(meanModel2.1b.1$Sol, meanModel2.1b.2$Sol, meanModel2.1b.3$Sol)) plot(meanModel2bSols) gelman.diag(meanModel2bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel2.1b.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel2.1b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel2.1b.1) # effective sample = 1000 posterior.mode(meanModel2.1b.1$Sol) # Zr = 0.24 HPDinterval(meanModel2.1b.1$Sol) # lwr = 0.01, upr = 0.55 (exp(2*0.24)-1) / (exp(2*0.24)+1) # r = 0.24 (small effect size) between <- posterior.mode(meanModel2.1b.1$VCV[,2]) repeated <- posterior.mode(meanModel2.1b.1$VCV[,1]) residual <- posterior.mode(meanModel2.1b.1$VCV[,4]) I2.between <- between / (between + repeated + residual + withinOff) * 100 # 7 % I2.repeated <- repeated / (between + repeated + residual + withinOff) * 100 # 5 % ### PART 2.2: BENEFITS OF MALE CARE (Zr Offspring) WITH PHYLOGENY ### ## A. Mean effect size ## # metafor # dropTip <- tree$tip.label[which(is.na(match(tree$tip.label, offData$species)))] offTree <- drop.tip(tree, dropTip, trim.internal=T) offTree <- makeNodeLabel(offTree, method = "number") offCor <- vcv.phylo(offTree, cor=TRUE) meanModel2.2a <- rma.mv(ZrOff ~ 1, V=VarOff, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyOff, ~ 1 | species), R = list(species=offCor), data=offData) summary(meanModel2.2a) # B = 0.25, lwr = 0.03, upr = 0.46, (Q = 252.4, p < 0.001) (exp(2*0.25)-1) / (exp(2*0.25)+1) # r = 0.25 (small effect size) I2.between <- 0.0332 / (0.2243 + 0.0204 + 0.0332 + 0.0000 + withinOff) * 100 # 11 % I2.repeated <- 0.0204 / (0.2243 + 0.0204 + 0.0332 + 0.0000 + withinOff) * 100 # 7 % I2.phylo <- 0.0000 / (0.2243 + 0.0204 + 0.0332 + 0.0000 + withinOff) * 100 # 0 % # MCMCglmm # offData$animal <- offData$species MEV <- offData$VarOff INtree <- inverseA(offTree, nodes="TIPS") meanModel2.2b <- MCMCglmm(ZrOff ~ 1, random = ~ID+StudyOff+animal, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE, pr=TRUE) meanModel2.2b.1 <- meanModel2.2b meanModel2.2b.2 <- meanModel2.2b meanModel2.2b.3 <- meanModel2.2b # chain convergence hist(meanModel2.2b.1$Liab) plot(meanModel2.2b.1$VCV) # all parameters well mixed plot(meanModel2.2b.1$Sol) # intercept estimate well mixed autocorr(meanModel2.2b.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(meanModel2.2b.1$Sol) # correlation between successive samples < 0.1 for all components meanModel2bSols <- mcmc.list(list(meanModel2.2b.1$Sol, meanModel2.2b.2$Sol, meanModel2.2b.3$Sol)) plot(meanModel2bSols) gelman.diag(meanModel2bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel2.2b.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel2.2b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel2.2b.1) # effective sample = 898 posterior.mode(meanModel2.2b.1$Sol) # intercept = 0.25 HPDinterval(meanModel2.2b.1$Sol) # intercept = -0.18 to 0.76 (exp(2*0.25)-1) / (exp(2*0.25)+1) # r = 0.25 (small effect size) # I^2 values between <- posterior.mode(meanModel2.2b.1$VCV[,2]) repeated <- posterior.mode(meanModel2.2b.1$VCV[,1]) phylo <- posterior.mode(meanModel2.2b.1$VCV[,3]) residual <- posterior.mode(meanModel2.2b.1$VCV[,5]) I2.between <- between / (between + repeated + phylo + residual + withinOff) * 100 # 4 % I2.repeated <- repeated / (between + repeated + phylo + residual + withinOff) * 100 # 4 % I2.phylo <- phylo/ (between + repeated + phylo + residual + withinOff) * 100 # 5 % ## B. Moderators ## # Differences between lab and field studies # MEV <- offData$VarOff INtree <- inverseA(tree, nodes="TIPS") meanModel2.2c <- MCMCglmm(ZrOff ~ labfieldOff-1, random = ~StudyOff+common+animal, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel2.2c) posterior.mode(meanModel2.2c$Sol) # 0.17, 0.30 HPDinterval(meanModel2.2c$Sol) # Differences between the amount and the probability of care # MEV <- offData$VarOff INtree <- inverseA(tree, nodes="TIPS") meanModel2.2d <- MCMCglmm(ZrOff ~ amountprobOff-1, random = ~StudyOff+common+animal, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel2.2d) posterior.mode(meanModel2.2d$Sol) # 0.30, 0.30 HPDinterval(meanModel2.2d$Sol) # Differences between observational and experimental studies # MEV <- offData$VarOff INtree <- inverseA(tree, nodes="TIPS") meanModel2.2e <- MCMCglmm(ZrOff ~ obsexptOff, random = ~StudyOff+common+animal, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel2.2e) posterior.mode(meanModel2.2e$Sol) # 0.32, 0.18 HPDinterval(meanModel2.2e$Sol) # Differences between the categories of care # MEV <- offData$VarOff INtree <- inverseA(tree, nodes="TIPS") meanModel2.2f <- MCMCglmm(ZrOff ~ carecatOff-1, random = ~StudyOff+common+animal, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel2.2f) posterior.mode(meanModel2.2f$Sol) # -0.11, -0.00, 0.05, 0.31, 0.33 HPDinterval(meanModel2.2f$Sol) table(meanModel2.2f$Sol[, 1] > meanModel2.2f$Sol[, 2]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.547, TRUE = 0.453 table(meanModel2.2f$Sol[, 1] > meanModel2.2f$Sol[, 3]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.504, TRUE = 0.496 table(meanModel2.2f$Sol[, 1] > meanModel2.2f$Sol[, 4]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.665, TRUE = 0.335 table(meanModel2.2f$Sol[, 1] > meanModel2.2f$Sol[, 5]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.661, TRUE = 0.339 table(meanModel2.2f$Sol[, 2] > meanModel2.2f$Sol[, 3]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.464, TRUE = 0.536 table(meanModel2.2f$Sol[, 2] > meanModel2.2f$Sol[, 4]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.745, TRUE = 0.255 table(meanModel2.2f$Sol[, 2] > meanModel2.2f$Sol[, 5]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.77, TRUE = 0.23 table(meanModel2.2f$Sol[, 3] > meanModel2.2f$Sol[, 4]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.663, TRUE = 0.337 table(meanModel2.2f$Sol[, 3] > meanModel2.2f$Sol[, 5]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.665, TRUE = 0.335 table(meanModel2.2f$Sol[, 4] > meanModel2.2f$Sol[, 5]) / length(meanModel2.2f$Sol[, 1]) # FALSE = 0.524, TRUE = 0.476 # Differences between the categories of offspring success # MEV <- offData$VarOff INtree <- inverseA(tree, nodes="TIPS") meanModel2.2g <- MCMCglmm(ZrOff ~ offcatOff-1, random = ~StudyOff+common+animal, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel2.2g) posterior.mode(meanModel2.2g$Sol) # 0.32, 0.36, 0.32 HPDinterval(meanModel2.2g$Sol) table(meanModel2.2g$Sol[, 1] > meanModel2.2g$Sol[, 2]) / length(meanModel2.2g$Sol[, 1]) # FALSE = 0.428, TRUE = 0.572 table(meanModel2.2g$Sol[, 1] > meanModel2.2g$Sol[, 3]) / length(meanModel2.2g$Sol[, 1]) # FALSE = 0.274, TRUE = 0.726 table(meanModel2.2g$Sol[, 2] > meanModel2.2g$Sol[, 3]) / length(meanModel2.2g$Sol[, 1]) # FALSE = 0.472, TRUE = 0.528 # Differences between brooding environments # MEV <- offData$VarOff INtree <- inverseA(tree, nodes="TIPS") meanModel2.2h <- MCMCglmm(ZrOff ~ brooding_environment-1, random = ~ID+StudyOff+animal, data=offData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel2.2h) posterior.mode(meanModel2.2h$Sol) # -0.1, 0.63, 0.50, 0.26, -0.02 HPDinterval(meanModel2.2h$Sol) table(meanModel2.2h$Sol[, 1] > meanModel2.2h$Sol[, 2]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.672, TRUE = 0.328 table(meanModel2.2h$Sol[, 1] > meanModel2.2h$Sol[, 3]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.704, TRUE = 0.296 table(meanModel2.2h$Sol[, 1] > meanModel2.2h$Sol[, 4]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.557, TRUE = 0.443 table(meanModel2.2h$Sol[, 1] > meanModel2.2h$Sol[, 5]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.484, TRUE = 0.516 table(meanModel2.2h$Sol[, 2] > meanModel2.2h$Sol[, 3]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.504, TRUE = 0.496 table(meanModel2.2h$Sol[, 2] > meanModel2.2h$Sol[, 4]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.314, TRUE = 0.686 table(meanModel2.2h$Sol[, 2] > meanModel2.2h$Sol[, 5]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.247, TRUE = 0.753 table(meanModel2.2h$Sol[, 3] > meanModel2.2h$Sol[, 4]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.202, TRUE = 0.798 table(meanModel2.2h$Sol[, 3] > meanModel2.2h$Sol[, 5]) / length(meanModel2.2h$Sol[, 1]) # FALSE = 0.13, TRUE = 0.87 ### PART 2.3: BENEFITS OF MALE CARE (Zr Offspring) WITH SYNGNATHIDS ADDED ### ## A. 1 Syngnathid ## # keep species name, ZrOff, VarOff, ID, StudyOff offData <- offData[,c("species","ZrOff","VarOff", "StudyOff")] offData$species <- as.character(offData$species) # needs to be character to add pipefish offData$StudyOff <- as.character(offData$StudyOff) # ibid. ## Create ZrOffspring Effect Size for one Syngnathid ## pipe <- data.frame(species = "pipe1", ZrOff = 0, VarOff = 0, StudyOff = "pfX") # these are values for the columns kept in offData offData2 <- rbind(offData, pipe) # new dataframe with syngnathid (pipfish) added offData2$ZrOff[37] <- 1/2 * log((1+0.999)/(1-0.999)) # note, can't use an r value of 1 here, = inf offData2$VarOff[37] <- mean(offData2$VarOff[1:36]) # for sampling variance, use average across studies offData2$species <- factor(offData2$species) # needed to be character to add pipefish offData2$obs <- paste(offData2$species, 1:length(offData2$species), sep = "") offData2$ID <- offData2$species offData2$animal <- offData2$species # to use in MCMCglmm ## Add One Pipefish to the Phylogeny ## # add pipefish to Syngnathus_typhle tree2 <- add.species.to.genus(tree, "Syngnathus_typhle") tree2$tip.label[47] <- "pipe1" # change name to match with that added to offData2 dropTip <- tree2$tip.label[which(is.na(match(tree2$tip.label, offData$species)))] offTree2 <- drop.tip(tree2, dropTip, trim.internal=T) offTree2 <- makeNodeLabel(offTree2, method = "number") offCor <- vcv.phylo(offTree2, cor=TRUE) ## Run the models from PART 2.2 ## # ZrOff within-study variance wjOff <- 1/offData2$VarOff kOff <- length(offData2$ZrOff) # number of estimates withinOff <- sum(wjOff * (kOff -1)) / (sum(wjOff)^2 - sum(wjOff^2)) # 0.022 # metafor # meanModel2.3a <- rma.mv(ZrOff ~ 1, V=VarOff, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyOff, ~ 1 | species), R = list(species=offCor), data=offData2) summary(meanModel2.3a) # B = 0.44, lwr = 0.07, upr = 0.82, (Q = 475.56, p < 0.001) (exp(2*0.44)-1) / (exp(2*0.44)+1) # r = 0.41 (medium effect size) I2.phylo <- 0.0000 / (0.3004 + 0.0627 + 0.3821 + 0.0000 + withinOff) * 100 # 0 % # MCMCglmm # priorB <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1,nu=0.02), G2=list(V=1,nu=0.02), G3=list(V=1,nu=0.02))) offData2$animal <- offData2$species MEV <- offData2$VarOff INtree <- inverseA(offTree2, nodes="TIPS") meanModel2.3b <- MCMCglmm(ZrOff ~ 1, random = ~ID+StudyOff+animal, data=offData2, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) meanModel2.3b.1 <- meanModel2.3b meanModel2.3b.2 <- meanModel2.3b meanModel2.3b.3 <- meanModel2.3b # chain convergence hist(meanModel2.3b.1$Liab) plot(meanModel2.3b.1$VCV) # all parameters well mixed plot(meanModel2.3b.1$Sol) # intercept estimate well mixed autocorr(meanModel2.3b.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(meanModel2.3b.1$Sol) # correlation between successive samples < 0.1 for all components meanModel2bSols <- mcmc.list(list(meanModel2.3b.1$Sol, meanModel2.3b.2$Sol, meanModel2.3b.3$Sol)) plot(meanModel2bSols) gelman.diag(meanModel2bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel2.3b.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel2.3b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel2.3b.1) # effective sample = 1000 posterior.mode(meanModel2.3b.1$Sol[,1]) # intercept = 0.43 HPDinterval(meanModel2.3b.1$Sol[,1]) # intercept = -0.52 to 1.33 (exp(2*0.43)-1) / (exp(2*0.43)+1) # r = 0.41 (medium effect size) # I^2 values between <- posterior.mode(meanModel2.3b.1$VCV[,2]) repeated <- posterior.mode(meanModel2.3b.1$VCV[,1]) phylo <- posterior.mode(meanModel2.3b.1$VCV[,3]) residual <- posterior.mode(meanModel2.3b.1$VCV[,5]) I2.between <- between / (between + repeated + phylo + residual + withinOff) * 100 # 4 % I2.repeated <- repeated / (between + repeated + phylo + residual + withinOff) * 100 # 5 % I2.phylo <- phylo/ (between + repeated + phylo + residual + withinOff) * 100 # 8 % ## B. 99 Syngnathids ## ## Create ZrOffspring Effect Size for 99 Syngnathids ## pipes <- data.frame(species = paste(rep("pipe"), 2:100, sep=""), ZrOff = rep(offData$ZrOff[37], 99), VarOff = rep(offData$VarOff[37], 99), StudyOff = paste(rep("pfX"), 2:100, sep="")) # diff study for each value assumes each comes from a different study offData <- zrData[!is.na(zrData$ZrOff),] # keep species name, ZrOff, VarOff, ID, StudyOff offData <- offData[,c("species","ZrOff","VarOff", "StudyOff")] offData$species <- as.character(offData$species) # needs to be character to add pipefish offData$StudyOff <- as.character(offData$StudyOff) # ibid. offData2 <- rbind(offData, pipes) offData2$species <- factor(offData2$species) # needed to be character to add pipefish offData2$obs <- paste(offData2$species, 1:length(offData2$species), sep = "") offData2$ID <- offData2$species offData2$animal <- offData2$species # to use in MCMCglmm ## Add 99 syngnathidae to the Phylogeny ## pAdds <- as.character(pipes$species) tree3 <- tree for(i in 1:length(pAdds)) {tree3 <- add.species.to.genus(tree3, pAdds[i], genus = "Syngnathus", where="random")} plotTree(tree3) tree3$tip.label # now has 147 species dropTip <- tree3$tip.label[which(is.na(match(tree3$tip.label, offData2$species)))] offTree3 <- drop.tip(tree3, dropTip, trim.internal=T) offTree3 <- makeNodeLabel(offTree3, method = "number") is.ultrametric(offTree3) offTree3 <- chronoMPL(offTree3) offCor <- vcv.phylo(offTree3, cor=TRUE) ## Run the model ## # ZrOff within-study variance (need to recalculate as more effect sizes now) wjOff <- 1/offData2$VarOff kOff <- length(offData2$ZrOff) # number of estimates withinOff <- sum(wjOff * (kOff -1)) / (sum(wjOff)^2 - sum(wjOff^2)) # 0.0399 # metafor # meanModel2.3c <- rma.mv(ZrOff ~ 1, V=VarOff, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyOff, ~ 1 | species), R = list(species=offCor), data=offData2) summary(meanModel2.3c) # B = 0.83, lwr = -0.97, upr = 2.62, (Q = 11080.3, p < 0.001) (exp(2*0.82)-1) / (exp(2*0.82)+1) # r = 0.68 (large effect size) I2.phylo <- 2.2825 / (2.2825 + 0.1787 + 0.0000 + 0.0000 + withinOff) * 100 # 0 % # MCMCglmm # offData2$animal <- offData2$species MEV <- offData2$VarOff INtree <- inverseA(offTree3, nodes="TIPS") meanModel2.3d <- MCMCglmm(ZrOff ~ 1, random = ~ID+StudyOff+animal, data=offData2, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) meanModel2.3d.1 <- meanModel2.3d meanModel2.3d.2 <- meanModel2.3d meanModel2.3d.3 <- meanModel2.3d # chain convergence hist(meanModel2.3d.1$Liab) plot(meanModel2.3d.1$VCV) # all parameters well mixed plot(meanModel2.3d.1$Sol) # intercept estimate well mixed autocorr(meanModel2.3d.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(meanModel2.3d.1$Sol) # correlation between successive samples < 0.1 for all components meanModel2bSols <- mcmc.list(list(meanModel2.3d.1$Sol, meanModel2.3d.2$Sol, meanModel2.3d.3$Sol)) plot(meanModel2bSols) gelman.diag(meanModel2bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel2.3d.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel2.3d.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel2.3d.1) # effective sample = 1000 posterior.mode(meanModel2.3d.1$Sol[,1]) # intercept = 0.61 HPDinterval(meanModel2.3d.1$Sol[,1]) # intercept = -0.84 to 2.92 (exp(2*0.46)-1) / (exp(2*0.46)+1) # r = 0.43 (medium effect size) # I^2 values between <- posterior.mode(meanModel2.3d.1$VCV[,2]) repeated <- posterior.mode(meanModel2.3d.1$VCV[,1]) phylo <- posterior.mode(meanModel2.3d.1$VCV[,3]) residual <- posterior.mode(meanModel2.3d.1$VCV[,5]) I2.between <- between / (between + repeated + phylo + residual + withinOff) * 100 # 1 % I2.repeated <- repeated / (between + repeated + phylo + residual + withinOff) * 100 # 0 % I2.phylo <- phylo/ (between + repeated + phylo + residual + withinOff) * 100 # 92 % # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 3.1: BENEFITS OF MALE CARE (Zr Preference1) WITHOUT PHYLOGENY ### ## A. Zr Preference1 Mean effect size ## # metafor # summary(rma(ZrPref1 ~ 1, vi=VarPref1, data=prefData1)) # no random effects: B = 0.47, lwr = 0.32, upr = 0.62 trimfill(rma(ZrPref1 ~ 1, vi=VarPref1, data=prefData1)) # 12 studies missing from the left-hand side funnel(trimfill(rma(ZrPref1 ~ 1, vi=VarPref1, data=prefData1))) # funnel plot with trim and fill prefData1$invSE <- 1/sqrt(prefData1$VarPref1) summary(lm(scale(ZrPref1) ~ invSE, data=prefData1)) # Egger's test: intercept = 0.17 (p = 0.44), slope = -0.04 (p = 0.33) meanModel3.1a <- rma.mv(ZrPref1 ~ 1, V=VarPref1, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyPref1), data = prefData1) summary(meanModel3.1a) # B = 0.50, lwr = 0.26, upr = 0.74 (Q = 485.45, p < 0.001) forest(meanModel3.1a) # within-study variance: withinCond = 0.029 I2.between <- 0.0863 / (0.0682 + 0.1883 + 0.0863 + withinPref1) * 100 # 23 % I2.repeated <- 0.1883 / (0.0682 + 0.1883 + 0.0863 + withinPref1) * 100 # 51 % # MCMCglmm # MEV <- prefData1$VarPref1 meanModel3.1b <- MCMCglmm(ZrPref1 ~ 1, random = ~ID+StudyPref1, data=prefData1, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorA, verbose = FALSE) # run the above three times meanModel3.1b.1 <- meanModel3.1b meanModel3.1b.2 <- meanModel3.1b meanModel3.1b.3 <- meanModel3.1b # chain convergence hist(meanModel3.1b.1$Liab) plot(meanModel3.1b.1$VCV) # all parameters well mixed plot(meanModel3.1b.1$Sol) # intercept estimate well mixed autocorr(meanModel3.1b.1$VCV) # correlation between successive samples < 0.1 for all parameters autocorr(meanModel3.1b.1$Sol) # correlation between successive samples < 0.1 for intercept meanModel3bSols <- mcmc.list(list(meanModel3.1b.1$Sol, meanModel3.1b.2$Sol, meanModel3.1b.3$Sol)) plot(meanModel3bSols) gelman.diag(meanModel3bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel3.1b.1$VCV) # units passed halfwidth heidel.diag(meanModel3.1b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel3.1b.1) # effective sample = 1000 posterior.mode(meanModel3.1b.1$Sol) # Zr = 0.49 HPDinterval(meanModel3.1b.1$Sol) # lwr = 0.26, upr = 0.74 (exp(2*0.49)-1) / (exp(2*0.49)+1) # r = 0.45 (medium effect size) between <- posterior.mode(meanModel3.1b.1$VCV[,2]) repeated <- posterior.mode(meanModel3.1b.1$VCV[,1]) residual <- posterior.mode(meanModel3.1b.1$VCV[,4]) I2.between <- between / (between + repeated + residual + withinPref1) * 100 # 38 % I2.repeated <- repeated / (between + repeated + residual + withinPref1) * 100 # 11 % ### PART 3.2: BENEFITS OF MALE CARE (Zr Preference1) WITH PHYLOGENY ### ## A. Zr Offspring Mean effect size ## # metafor # dropTip <- tree$tip.label[which(is.na(match(tree$tip.label, prefData1$species)))] prefTree1 <- drop.tip(tree, dropTip, trim.internal=T) prefTree1 <- makeNodeLabel(prefTree1, method = "number") prefCor1 <- vcv.phylo(prefTree1, cor=TRUE) meanModel3.2a <- rma.mv(ZrPref1 ~ 1, V=VarPref1, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyPref1, ~ 1 | species), R = list(species=prefCor1), data=prefData1) summary(meanModel3.2a) # B = 0.50, lwr = 0.26, upr = 0.74, (Q = 485.45, p < 0.001) (exp(2*0.5)-1) / (exp(2*0.5)+1) # r = 0.46 (large effect size) I2.between <- 0.0863 / (0.0682 + 0.1883 + 0.0863 + 0.0000 + withinPref1) * 100 # 23 % I2.repeated <- 0.1883 / (0.0682 + 0.1883 + 0.0863 + 0.0000 + withinPref1) * 100 # 51 % I2.phylo <- 0.0000 / (0.0682 + 0.1883 + 0.0863 + 0.0000 + withinPref1) * 100 # 0 % # MCMCglmm # prefData1$animal <- prefData1$species priorB <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1,nu=0.02),G2=list(V=1,nu=0.02),G3=list(V=1,nu=0.02))) MEV <- prefData1$VarPref1 INtree <- inverseA(prefTree1, nodes="TIPS") meanModel3.2b <- MCMCglmm(ZrPref1 ~ 1, random = ~ID+StudyPref1+animal, data=prefData1, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE, pr=TRUE) meanModel3.2b.1 <- meanModel3.2b meanModel3.2b.2 <- meanModel3.2b meanModel3.2b.3 <- meanModel3.2b # chain convergence hist(meanModel3.2b.1$Liab) plot(meanModel3.2b.1$VCV) # all parameters well mixed plot(meanModel3.2b.1$Sol) # intercept estimate well mixed autocorr(meanModel3.2b.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(meanModel3.2b.1$Sol) # correlation between successive samples < 0.1 for all components meanModel3bSols <- mcmc.list(list(meanModel3.2b.1$Sol, meanModel3.2b.2$Sol, meanModel3.2b.3$Sol)) plot(meanModel3bSols) gelman.diag(meanModel3bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel3.2b.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel3.2b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel3.2b.1) # effective sample = 1000 posterior.mode(meanModel3.2b.1$Sol) # intercept = 0.51 HPDinterval(meanModel3.2b.1$Sol) # intercept = 0.04 to 1.38 (exp(2*0.51)-1) / (exp(2*0.51)+1) # r = 0.47 (medium effect size) # I^2 values between <- posterior.mode(meanModel3.2b.1$VCV[,2]) repeated <- posterior.mode(meanModel3.2b.1$VCV[,1]) phylo <- posterior.mode(meanModel3.2b.1$VCV[,3]) residual <- posterior.mode(meanModel3.2b.1$VCV[,5]) I2.between <- between / (between + repeated + phylo + residual + withinPref1) * 100 # 17 % I2.repeated <- repeated / (between + repeated + phylo + residual + withinPref1) * 100 # 23 % I2.phylo <- phylo/ (between + repeated + phylo + residual + withinPref1) * 100 # 11 % ## B. Moderators ## # Differences between lab and field studies # MEV <- prefData1$VarPref1 INtree <- inverseA(tree, nodes="TIPS") meanModel3.2c <- MCMCglmm(ZrPref1 ~ labfieldPref1, random = ~StudyPref1+common+animal, data=prefData1, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel3.2c) posterior.mode(meanModel3.2c$Sol) # 0.56, 0.52 HPDinterval(meanModel3.2c$Sol) # Differences between studies where brood size was natural or manipulated # MEV <- prefData1$VarPref1 INtree <- inverseA(tree, nodes="TIPS") meanModel3.2d <- MCMCglmm(ZrPref1 ~ manipnatPref1, random = ~StudyPref1+common+animal, data=prefData1, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel3.2d) posterior.mode(meanModel3.2d$Sol) # 0.51, 0.49 HPDinterval(meanModel3.2d$Sol) # Differences between observational and experimental studies # MEV <- prefData1$VarPref1 INtree <- inverseA(tree, nodes="TIPS") meanModel3.2e <- MCMCglmm(ZrPref1 ~ obsexptPref1, random = ~StudyPref1+common+animal, data=prefData1, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel3.2e) posterior.mode(meanModel3.2e$Sol) # 0.54, 0.68 HPDinterval(meanModel3.2e$Sol) # Differences between the categories of female Preference # MEV <- prefData1$VarPref1 INtree <- inverseA(tree, nodes="TIPS") meanModel3.2f <- MCMCglmm(ZrPref1 ~ prefcatPref1-1, random = ~StudyPref1+common+animal, data=prefData1, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel3.2f) posterior.mode(meanModel3.2f$Sol) # 0.63, 0.60, 0.63 HPDinterval(meanModel3.2f$Sol) table(meanModel3.2f$Sol[, 1] > meanModel3.2f$Sol[, 2]) / length(meanModel3.2f$Sol[, 1]) # FALSE = 0.45, TRUE = 0.55 table(meanModel3.2f$Sol[, 1] > meanModel3.2f$Sol[, 3]) / length(meanModel3.2f$Sol[, 1]) # FALSE = 0.46, TRUE = 0.54 table(meanModel3.2f$Sol[, 2] > meanModel3.2f$Sol[, 3]) / length(meanModel3.2f$Sol[, 1]) # FALSE = 0.52, TRUE = 0.48 # Differences between egg removal and natural egg absence # MEV <- prefData1$VarPref1 INtree <- inverseA(tree, nodes="TIPS") meanModel3.2g <- MCMCglmm(ZrPref1 ~ eggremovePref1, random = ~StudyPref1+common+animal, data=prefData1, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel3.2g) posterior.mode(meanModel3.2g$Sol) # 0.64, 0.55 HPDinterval(meanModel3.2g$Sol) # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 4.1: BENEFITS OF MALE CARE (Zr Preference2) WITHOUT PHYLOGENY ### ## A. Zr Preference Mean effect size ## # metafor # summary(rma(ZrPref2 ~ 1, vi=VarPref2, data=prefData2)) # no random effects: B = 0.32, lwr = 0.10, upr = 0.53 trimfill(rma(ZrPref2 ~ 1, vi=VarPref2, data=prefData2)) # 0 studies missing funnel(trimfill(rma(ZrPref2 ~ 1, vi=VarPref2, data=prefData2))) # funnel plot with trim and fill prefData2$invSE <- 1/sqrt(prefData2$VarPref2) summary(lm(scale(ZrPref2) ~ invSE, data=prefData2)) # Egger's test: intercept = -0.14 (p = 0.76), slope = 0.03 (p = 0.74) meanModel4.1a <- rma.mv(ZrPref2 ~ 1, V=VarPref2, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyPref2), data = prefData2) summary(meanModel4.1a) # B = 0.17, lwr = -0.16, upr = 0.50 (Q = 194.514, p < 0.001) forest(meanModel4.1a) # within-study variance: withinCond = 0.043 I2.between <- 0.0000 / (0.0947 + 0.0000 + 0.2309 + withinPref2) * 100 # 0 % I2.repeated <- 0.2309 / (0.0947 + 0.0000 + 0.2309 + withinPref2) * 100 # 62 % # MCMCglmm # MEV <- prefData2$VarPref2 meanModel4.1b <- MCMCglmm(ZrPref2 ~ 1, random = ~ID+StudyPref2, data=prefData2, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorA, verbose = FALSE) # run the above three times meanModel4.1b.1 <- meanModel4.1b meanModel4.1b.2 <- meanModel4.1b meanModel4.1b.3 <- meanModel4.1b # chain convergence hist(meanModel4.1b.1$Liab) plot(meanModel4.1b.1$VCV) # all parameters well mixed plot(meanModel4.1b.1$Sol) # intercept estimate well mixed autocorr(meanModel4.1b.1$VCV) # correlation between successive samples < 0.1 for all parameters autocorr(meanModel4.1b.1$Sol) # correlation between successive samples < 0.1 for intercept meanModel3bSols <- mcmc.list(list(meanModel4.1b.1$Sol, meanModel4.1b.2$Sol, meanModel4.1b.3$Sol)) plot(meanModel3bSols) gelman.diag(meanModel3bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel4.1b.1$VCV) # units passed halfwidth heidel.diag(meanModel4.1b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel4.1b.1) # effective sample = 1000 posterior.mode(meanModel4.1b.1$Sol) # Zr = 0.27 HPDinterval(meanModel4.1b.1$Sol) # lwr = -0.15, upr = 0.62 (exp(2*0.27)-1) / (exp(2*0.27)+1) # r = 0.26 (medium effect size) between <- posterior.mode(meanModel4.1b.1$VCV[,2]) repeated <- posterior.mode(meanModel4.1b.1$VCV[,1]) residual <- posterior.mode(meanModel4.1b.1$VCV[,4]) I2.between <- between / (between + repeated + residual + withinPref2) * 100 # 5 % I2.repeated <- repeated / (between + repeated + residual + withinPref2) * 100 # 43 % ### PART 4.2: BENEFITS OF MALE CARE (Zr Preference2) WITH PHYLOGENY ### ## A. Zr Preference2 Mean effect size ## # metafor # dropTip <- tree$tip.label[which(is.na(match(tree$tip.label, prefData2$species)))] prefTree2 <- drop.tip(tree, dropTip, trim.internal=T) prefTree2 <- makeNodeLabel(prefTree2, method = "number") prefTree2<-chronoMPL(prefTree2) prefCor2 <- vcv.phylo(prefTree2, cor=TRUE) meanModel4.2a <- rma.mv(ZrPref2 ~ 1, V=VarPref2, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyPref2, ~ 1 | species), R = list(species=prefCor2), data=prefData2) summary(meanModel4.2a) # B = 0.14, lwr = -0.35, upr = 0.63, (Q = 194.51, p < 0.001) (exp(2*0.14)-1) / (exp(2*0.14)+1) # r = 0.14 (large effect size) I2.between <- 0.0145 / (0.0984 + 0.0145 + 0.0530 + 0.2617 + withinPref2) * 100 # 3 % I2.repeated <- 0.0530 / (0.0984 + 0.0145 + 0.0530 + 0.2617 + withinPref2) * 100 # 11 % I2.phylo <- 0.2627 / (0.0984 + 0.0145 + 0.0530 + 0.2617 + withinPref2) * 100 # 56 % # MCMCglmm # prefData2$animal <- prefData2$species priorB <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1,nu=0.02),G2=list(V=1,nu=0.02),G3=list(V=1,nu=0.02))) MEV <- prefData2$VarPref2 INtree <- inverseA(prefTree2, nodes="TIPS") meanModel4.2b <- MCMCglmm(ZrPref2 ~ 1, random = ~ID+StudyPref2+animal, data=prefData2, pl=TRUE, slice=TRUE,nitt = 1200000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE, pr=TRUE) meanModel4.2b.1 <- meanModel4.2b meanModel4.2b.2 <- meanModel4.2b meanModel4.2b.3 <- meanModel4.2b # chain convergence hist(meanModel4.2b.1$Liab) plot(meanModel4.2b.1$VCV) # all parameters well mixed plot(meanModel4.2b.1$Sol) # intercept estimate well mixed autocorr(meanModel4.2b.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(meanModel4.2b.1$Sol) # correlation between successive samples < 0.1 for all components meanModel3bSols <- mcmc.list(list(meanModel4.2b.1$Sol, meanModel4.2b.2$Sol, meanModel4.2b.3$Sol)) plot(meanModel3bSols) gelman.diag(meanModel3bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel4.2b.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel4.2b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel4.2b.1) # effective sample = 1100 posterior.mode(meanModel4.2b.1$Sol) # intercept = 0.19 HPDinterval(meanModel4.2b.1$Sol) # intercept = -0.32 to 0.70 (exp(2*0.19)-1) / (exp(2*0.19)+1) # r = 0.19 (small effect size) # I^2 values between <- posterior.mode(meanModel4.2b.1$VCV[,2]) repeated <- posterior.mode(meanModel4.2b.1$VCV[,1]) phylo <- posterior.mode(meanModel4.2b.1$VCV[,3]) residual <- posterior.mode(meanModel4.2b.1$VCV[,5]) I2.between <- between / (between + repeated + phylo + residual + withinPref2) * 100 # 7 % I2.repeated <- repeated / (between + repeated + phylo + residual + withinPref2) * 100 # 7 % I2.phylo <- phylo/ (between + repeated + phylo + residual + withinPref2) * 100 # 10 % ## B. Moderators ## # Differences between lab and field studies # MEV <- prefData2$VarPref2 INtree <- inverseA(tree, nodes="TIPS") meanModel4.2c <- MCMCglmm(ZrPref2 ~ labfieldPref2-1, random = ~StudyPref2+common+animal, data=prefData2, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel4.2c) posterior.mode(meanModel4.2c$Sol) HPDinterval(meanModel4.2c$Sol) # Differences between studies where brood size was natural or manipulated # MEV <- prefData2$VarPref2 INtree <- inverseA(tree, nodes="TIPS") meanModel4.2d <- MCMCglmm(ZrPref2 ~ manipnatPref2-1, random = ~StudyPref2+common+animal, data=prefData2, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel4.2d) posterior.mode(meanModel4.2d$Sol) # 0.12, 0.11 HPDinterval(meanModel4.2d$Sol) # Differences between observational and experimental studies # MEV <- prefData2$VarPref2 INtree <- inverseA(tree, nodes="TIPS") meanModel4.2e <- MCMCglmm(ZrPref2 ~ obsexptPref2-1, random = ~StudyPref2+common+animal, data=prefData2, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel4.2e) posterior.mode(meanModel4.2e$Sol) # 0.32, 0.13 HPDinterval(meanModel4.2e$Sol) # Differences between the categories of female preference # MEV <- prefData2$VarPref2 INtree <- inverseA(tree, nodes="TIPS") meanModel4.2f <- MCMCglmm(ZrPref2 ~ prefcatPref2-1, random = ~StudyPref2+common+animal, data=prefData2, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel4.2f) posterior.mode(meanModel4.2f$Sol) # 0.35, 0.25, 0.28 HPDinterval(meanModel4.2f$Sol) table(meanModel4.2f$Sol[, 1] > meanModel4.2f$Sol[, 2]) / length(meanModel4.2f$Sol[, 1]) # FALSE = 0.241, TRUE = 0.759 table(meanModel4.2f$Sol[, 1] > meanModel4.2f$Sol[, 3]) / length(meanModel4.2f$Sol[, 1]) # FALSE = 0.26, TRUE = 0.74 table(meanModel4.2f$Sol[, 2] > meanModel4.2f$Sol[, 3]) / length(meanModel4.2f$Sol[, 1]) # FALSE = 0.582, TRUE = 0.418 # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 5.1: THE EFFECT OF BROOD SIZE ON MALE CARE (Zr Brood size) WITHOUT PHYLOGENY ### ## A. Mean effect size ## # metafor # summary(rma(ZrBrood ~ 1, vi=VarBrood, data=broodData)) # no random effects: B = 0.36, lwr = 0.30, upr = 0.43 trimfill(rma(ZrBrood ~ 1, vi=VarBrood, data=broodData)) # 0 studies missing from the left-hand side funnel(trimfill(rma(ZrBrood ~ 1, vi=VarBrood, data=broodData))) # funnel plot with trim and fill broodData$invSE <- 1/sqrt(broodData$VarBrood) summary(lm(scale(ZrBrood) ~ invSE, data=broodData)) # Egger's test: intercept = 0.13 (p = 0.56), slope = -0.02 (p = 0.50) meanModel5.1a <- rma.mv(ZrBrood ~ 1, V=VarBrood, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyBrood), data = broodData) summary(meanModel5.1a) # B = 0.41, lwr = 0.31, upr = 0.51 (Q = 474.43, p < 0.001) forest(meanModel5.1a) # within-study variance: withinCond = 0.019 I2.between <- 0.0228 / (0.0228 + 0.0297 + 0.0134 + withinBrood) * 100 # 27 % I2.repeated <- 0.0297 / (0.0228 + 0.0297 + 0.0134 + withinBrood) * 100 # 35 % # MCMCglmm # MEV <- broodData$VarBrood meanModel5.1b <- MCMCglmm(ZrBrood ~ 1, random = ~ID+StudyBrood, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorA, verbose = FALSE) # run the above three times meanModel5.1b.1 <- meanModel5.1b meanModel5.1b.2 <- meanModel5.1b meanModel5.1b.3 <- meanModel5.1b # chain convergence hist(meanModel5.1b.1$Liab) plot(meanModel5.1b.1$VCV) # all parameters well mixed plot(meanModel5.1b.1$Sol) # intercept estimate well mixed autocorr(meanModel5.1b.1$VCV) # correlation between successive samples < 0.1 for all parameters autocorr(meanModel5.1b.1$Sol) # correlation between successive samples < 0.1 for intercept meanModel4bSols <- mcmc.list(list(meanModel5.1b.1$Sol, meanModel5.1b.2$Sol, meanModel5.1b.3$Sol)) plot(meanModel4bSols) gelman.diag(meanModel4bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel5.1b.1$VCV) # units passed halfwidth heidel.diag(meanModel5.1b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel5.1b.1) # effective sample = 893 posterior.mode(meanModel5.1b.1$Sol) # Zr = 0.36 HPDinterval(meanModel5.1b.1$Sol) # lwr = 0.31, upr = 0.52 (exp(2*0.36)-1) / (exp(2*0.36)+1) # r = 0.35 (medium effect size) between <- posterior.mode(meanModel5.1b.1$VCV[,2]) repeated <- posterior.mode(meanModel5.1b.1$VCV[,1]) residual <- posterior.mode(meanModel5.1b.1$VCV[,4]) I2.between <- between / (between + repeated + residual + withinBrood) * 100 # 24 % I2.repeated <- repeated / (between + repeated + residual + withinBrood) * 100 # 38 % ### PART 5.2: THE EFFECT OF BROOD SIZE ON MALE CARE (Zr Brood Size) WITH PHYLOGENY ### ## A. Mean effect size ## # metafor # dropTip <- tree$tip.label[which(is.na(match(tree$tip.label, broodData$species)))] broodTree <- drop.tip(tree, dropTip, trim.internal=T) broodTree <- makeNodeLabel(broodTree, method = "number") broodCor <- vcv.phylo(broodTree, cor=TRUE) meanModel5.2a <- rma.mv(ZrBrood ~ 1, V=VarBrood, random = list(~ 1 | obs, ~ 1 | ID, ~ 1 | StudyBrood, ~ 1 | species), R = list(species=broodCor), data=broodData) summary(meanModel5.2a) # B = 0.48, lwr = 0.23, upr = 0.72, (Q = 476.5, p = 0.0001) (exp(2*0.48)-1) / (exp(2*0.48)+1) # r = 0.45 (medium effect size) I2.between <- 0.0232 / (0.0232 + 0.0207 + 0.0255 + 0.0139 + withinBrood) * 100 # 23 % I2.repeated <- 0.0207 / (0.0232 + 0.0207 + 0.0255 + 0.0139 + withinBrood) * 100 # 20 % I2.phylo <- 0.0255 / (0.0232 + 0.0207 + 0.0255 + 0.0139 + withinBrood) * 100 # 25 % # MCMCglmm # broodData$animal <- broodData$species MEV <- broodData$VarBrood INtree <- inverseA(broodTree, nodes="TIPS") meanModel5.2b <- MCMCglmm(ZrBrood ~ 1, random = ~ID+StudyBrood+animal, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE, pr=TRUE) meanModel5.2b.1 <- meanModel5.2b meanModel5.2b.2 <- meanModel5.2b meanModel5.2b.3 <- meanModel5.2b # chain convergence hist(meanModel5.2b.1$Liab) plot(meanModel5.2b.1$VCV) # all parameters well mixed plot(meanModel5.2b.1$Sol) # intercept estimate well mixed autocorr(meanModel5.2b.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(meanModel5.2b.1$Sol) # correlation between successive samples < 0.1 for all components meanModel4bSols <- mcmc.list(list(meanModel5.2b.1$Sol, meanModel5.2b.2$Sol, meanModel5.2b.3$Sol)) plot(meanModel4bSols) gelman.diag(meanModel4bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(meanModel5.2b.1$VCV) # all parameters passed halfwidth heidel.diag(meanModel5.2b.1$Sol) # intercept passed halfwidth # model parameters summary(meanModel5.2b.1) # effective sample = 1000 posterior.mode(meanModel5.2b.1$Sol) # intercept = 0.47 HPDinterval(meanModel5.2b.1$Sol) # intercept = 0.19 to 0.86 (exp(2*0.47)-1) / (exp(2*0.47)+1) # r = 0.44 (medium effect size) # I^2 values between <- posterior.mode(meanModel5.2b.1$VCV[,2]) repeated <- posterior.mode(meanModel5.2b.1$VCV[,1]) phylo <- posterior.mode(meanModel5.2b.1$VCV[,3]) residual <- posterior.mode(meanModel5.2b.1$VCV[,5]) I2.between <- between / (between + repeated + phylo + residual + withinBrood) * 100 # 20 % I2.repeated <- repeated / (between + repeated + phylo + residual + withinBrood) * 100 # 15 % I2.phylo <- phylo/ (between + repeated + phylo + residual + withinBrood) * 100 # 18 % ## B. Moderators ## # Differences between lab and field studies # MEV <- broodData$VarBrood INtree <- inverseA(tree, nodes="TIPS") meanModel5.2c <- MCMCglmm(ZrBrood ~ labfieldBrood-1, random = ~StudyBrood+common+animal, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel5.2c) posterior.mode(meanModel5.2c$Sol) HPDinterval(meanModel5.2c$Sol) # Differences between observational and experimental studies # MEV <- broodData$VarBrood INtree <- inverseA(tree, nodes="TIPS") meanModel5.2d <- MCMCglmm(ZrBrood ~ obsexptBrood-1, random = ~StudyBrood+common+animal, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel5.2d) posterior.mode(meanModel5.2d$Sol) HPDinterval(meanModel5.2d$Sol) # Differences between the amount and the probability of care # MEV <- broodData$VarBrood INtree <- inverseA(tree, nodes="TIPS") meanModel5.2e <- MCMCglmm(ZrBrood ~ amountprobBrood-1, random = ~StudyBrood+common+animal, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel5.2e) posterior.mode(meanModel5.2e$Sol) HPDinterval(meanModel5.2e$Sol) # Differences between studies where brood size was natural or manipulated # MEV <- broodData$VarBrood INtree <- inverseA(tree, nodes="TIPS") meanModel5.2f <- MCMCglmm(ZrBrood ~ manipnatBrood-1, random = ~StudyBrood+common+animal, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel5.2f) posterior.mode(meanModel5.2f$Sol) HPDinterval(meanModel5.2f$Sol) # Differences between the categories of brood size # MEV <- broodData$VarBrood INtree <- inverseA(tree, nodes="TIPS") meanModel5.2g <- MCMCglmm(ZrBrood ~ sizecatBrood-1, random = ~StudyBrood+common+animal, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel5.2g) posterior.mode(meanModel5.2g$Sol) HPDinterval(meanModel5.2g$Sol) table(meanModel5.2g$Sol[, 1] > meanModel5.2g$Sol[, 2]) / length(meanModel5.2g$Sol[, 1]) # FALSE = 0.307, TRUE = 0.693 table(meanModel5.2g$Sol[, 1] > meanModel5.2g$Sol[, 3]) / length(meanModel5.2g$Sol[, 1]) # FALSE = 0.274, TRUE = 0.726 table(meanModel5.2g$Sol[, 2] > meanModel5.2g$Sol[, 3]) / length(meanModel5.2g$Sol[, 1]) # FALSE = 0.379, TRUE = 0.621 # Differences between the categories of care # MEV <- broodData$VarBrood INtree <- inverseA(tree, nodes="TIPS") meanModel5.2h <- MCMCglmm(ZrBrood ~ carecatBrood-1, random = ~StudyBrood+common+animal, data=broodData, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorB, ginverse = list(animal = INtree$Ainv), verbose = FALSE) summary(meanModel5.2h) posterior.mode(meanModel5.2h$Sol) HPDinterval(meanModel5.2h$Sol) table(meanModel5.2h$Sol[, 1] > meanModel5.2h$Sol[, 2]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.38, TRUE = 0.62 table(meanModel5.2h$Sol[, 1] > meanModel5.2h$Sol[, 3]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.14, TRUE = 0.86 table(meanModel5.2h$Sol[, 1] > meanModel5.2h$Sol[, 4]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.365, TRUE = 0.635 table(meanModel5.2h$Sol[, 1] > meanModel5.2h$Sol[, 5]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.257, TRUE = 0.743 table(meanModel5.2h$Sol[, 2] > meanModel5.2h$Sol[, 3]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.209, TRUE = 0.791 table(meanModel5.2h$Sol[, 2] > meanModel5.2h$Sol[, 4]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.441, TRUE = 0.559 table(meanModel5.2h$Sol[, 2] > meanModel5.2h$Sol[, 5]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.357, TRUE = 0.643 table(meanModel5.2h$Sol[, 3] > meanModel5.2h$Sol[, 4]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.686, TRUE = 0.314 table(meanModel5.2h$Sol[, 3] > meanModel5.2h$Sol[, 5]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.697, TRUE = 0.303 table(meanModel5.2h$Sol[, 4] > meanModel5.2h$Sol[, 5]) / length(meanModel5.2h$Sol[, 1]) # FALSE = 0.431, TRUE = 0.569 # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 6: FEMALE PREFERENCE FOR CARING MALES ### ## Dataframe with only species present in both condData and prefData1 # remove species not in other dataset prefSpecies3<-prefSpecies1[which((prefSpecies1$species %in% broodSpecies$species)==T),] broodSpecies2<-broodSpecies[which((broodSpecies$species %in% prefSpecies1$species)==T),] # combine data set into one data frame PrefBrood<-data.frame(prefSpecies3,broodSpecies2) ## 6.1 WITHOUT PHYLOGENY ## # metafor # Model6.1a <- rma.mv(ZrPref1 ~ ZrBrood, V=VarPref1, random = list(~ 1 | obs), data = PrefBrood) summary(Model6.1a) # B = 1.36, lwr = 0.46, upr = 2.27 (Q = 81.42, p < 0.001) Model6.1a <- rma.mv(ZrPref1 ~ ZrBrood + log(Nbrood), V=VarPref1, random = list(~ 1 | obs), data = PrefBrood) summary(Model6.1a) # B = 1.27, lwr = 0.33, upr = 2.21 (Q = 81.42, p < 0.001) # MCMCglmm # priorA <- list(R=list(V=1,nu=0.02), G=list(G1=list(V=1,nu=0.02), G2=list(V=1,nu=0.02))) MEV <- PrefBrood$VarPref Model6.1b <- MCMCglmm(ZrPref1 ~ ZrBrood + log(Nbrood), data=PrefBrood, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorA, verbose = FALSE) # run the above three times Model6.1b.1 <- Model6.1b Model6.1b.2 <- Model6.1b Model6.1b.3 <- Model6.1b # chain convergence hist(Model6.1b.1$Liab) plot(Model6.1b.1$VCV) # residual var well mixed plot(Model6.1b.1$Sol) # intercept and slope estimates well mixed autocorr(Model6.1b.1$VCV) # correlation between successive samples < 0.1 for all parameters autocorr(Model6.1b.1$Sol) # correlation between successive samples < 0.1 for intercept Model6.1bSols <- mcmc.list(list(Model6.1b.1$Sol, Model6.1b.2$Sol, Model6.1b.3$Sol)) plot(Model6.1bSols) gelman.diag(Model6.1bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(Model6.1b.1$VCV) # units passed halfwidth heidel.diag(Model6.1b.1$Sol) # intercept and slopes passed halfwidth # model parameters summary(Model6.1b.1) # effective sample = 1000 for all three parameters posterior.mode(Model6.1b.1$Sol) # ZrBrood = 1.21 HPDinterval(Model6.1b.1$Sol) # ZrBrood lwr = 0.22, ZrBrood upr = 2.22 ## 6.2 WITH PHYLOGENY ## # Metafor # dropTip <- tree$tip.label[which(is.na(match(tree$tip.label, PrefBrood$species)))] pdTree <- drop.tip(tree, dropTip, trim.internal=T) pdTree <- makeNodeLabel(pdTree, method = "number") pdCor <- vcv.phylo(pdTree, cor=TRUE) Model6.2a <- rma.mv(ZrPref1 ~ ZrBrood + log(Nbrood), V=VarPref1, random = list(~ 1 | obs, ~ 1 | species), R = list(species=pdCor), data=PrefBrood) summary(Model6.2a) # B = 1.27, lwr = 0.33, upr = 2.21 (Q = 81.42, p < 0.001) # MCMCglmm # PrefBrood$animal <- PrefBrood$species priorC <- list(R=list(R1=list(V=1,nu=0.02), R2=list(V=1,nu=0.02)), G=list(G1=list(V=1,nu=0.02))) MEV <- PrefBrood$VarPref1 INtree <- inverseA(pdTree, nodes="TIPS") Model6.2b <- MCMCglmm(ZrPref1 ~ ZrBrood + log(Nbrood), random = ~animal, data=PrefBrood, pl=TRUE, slice=TRUE,nitt = 1100000, burnin = 100000, thin = 1000, mev = MEV, prior = priorC, ginverse = list(animal = INtree$Ainv), verbose = FALSE) Model6.2b.1 <- Model6.2b Model6.2b.2 <- Model6.2b Model6.2b.3 <- Model6.2b # chain convergence hist(Model6.2b.1$Liab) plot(Model6.2b.1$VCV) # all parameters well mixed plot(Model6.2b.1$Sol) # intercept and slope estimates well mixed autocorr(Model6.2b.1$VCV) # correlation between successive samples < 0.1 for all components autocorr(Model6.2b.1$Sol) # correlation between successive samples < 0.1 for all components Model6.2bSols <- mcmc.list(list(Model6.2b.1$Sol, Model6.2b.2$Sol, Model6.2b.3$Sol)) plot(Model6.2bSols) gelman.diag(Model6.2bSols) # upper CI = 1.00 for intercept suggesting convergence heidel.diag(Model6.2b.1$VCV) # all parameters passed halfwidth heidel.diag(Model6.2b.1$Sol) # intercept passed halfwidth # model parameters summary(Model6.2b.1) posterior.mode(Model6.2b.1$Sol) # ZrBrood = 1.27 HPDinterval(Model6.2b.1$Sol) # ZrBrood lwr = 0.44, ZrBrood upr = 2.47 # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 7: FIGURES ### ## Figure 2: ZrCondition and ZrOffspring forest plots ## condSpecies$common<-gsub("_", " ", condSpecies$common) condSpecies$Numstudy<-paste("(", condSpecies$Nstudy, ")", sep="") condSpecies<-transform(condSpecies, common_study=paste(common, Numstudy, sep=" ")) condEffects<-group_by(condData, species) condEffects<-summarise(condEffects, Neffects=length(ZrCond)) condSpecies$Neffects<-condEffects$Neffects offSpecies$common<-gsub("_", " ", offSpecies$common) offSpecies$Numstudy<-paste("(", offSpecies$Nstudy, ")", sep="") offSpecies<-transform(offSpecies, common_study=paste(common, Numstudy, sep=" ")) offEffects<-group_by(offData, species) offEffects<-summarise(offEffects, Neffects=length(ZrOff)) offSpecies$Neffects<-offEffects$Neffects fp1<-ggplot(condSpecies, aes(reorder(common_study,-ZrCond),ZrCond))+ geom_point(fill="grey",aes(size=Neffects),shape=16,col="black")+ ylim(-2,1)+ geom_hline(yintercept=-0.28, linetype=2)+ geom_errorbar(ymin=condSpecies$lciCond, ymax=condSpecies$uciCond,size=0.2,width=0.3,col="black")+ geom_rect(aes(ymin=-0.62, ymax=0.08, xmin=0, xmax=Inf), alpha=0.01)+ geom_hline(yintercept = 0, linetype=1)+ ggtitle("(a)")+ coord_flip()+ xlab("")+ ylab(expression(italic("Zr")[condition]))+ theme_bw(base_size=11)+ theme(legend.position="none", plot.title=element_text(hjust=-0.78)) fp2<-ggplot(offSpecies, aes(reorder(common_study,-ZrOff),ZrOff))+ geom_point(aes(size=Neffects),col="Black")+ geom_errorbar(aes(ymin=lciOff, ymax=uciOff),size=0.2,width=0.3,col="black")+ ylim(-1.5,2)+ geom_hline(yintercept = 0, linetype=1)+ geom_hline(yintercept = 0.25, linetype=2)+ geom_rect(aes(ymin=-0.18, ymax=0.76, xmin=0, xmax=Inf), alpha=0.01)+ ggtitle("(b)")+ coord_flip()+ xlab("")+ ylab(expression(italic("Zr")[offspring]))+ theme_bw(base_size=11)+ theme(legend.position="none",plot.title=element_text(hjust=-0.8)) grid.arrange(fp1,fp2, ncol=1, nrow=2) ## Figure 3: ZrPreference1 and ZrPreference 2 forest plot ## prefSpecies1$common<-gsub("_", " ", prefSpecies1$common) prefSpecies1$Numstudy<-paste("(", prefSpecies1$Nstudy, ")", sep="") prefSpecies1<-transform(prefSpecies1, common_study=paste(common, Numstudy, sep=" ")) pref1Effects<-group_by(prefData1, species) pref1Effects<-summarise(pref1Effects, Neffects=length(ZrPref1)) prefSpecies1$Neffects<-pref1Effects$Neffects prefSpecies2$common<-gsub("_", " ", prefSpecies2$common) prefSpecies2$Numstudy<-paste("(", prefSpecies2$Nstudy, ")", sep="") prefSpecies2<-transform(prefSpecies2, common_study=paste(common, Numstudy, sep=" ")) pref2Effects<-group_by(prefData2, species) pref2Effects<-summarise(pref2Effects, Neffects=length(ZrPref2)) prefSpecies2$Neffects<-pref2Effects$Neffects fp3<-ggplot(prefSpecies1, aes(reorder(common_study,-ZrPref1),ZrPref1))+ geom_point(aes(size=Neffects))+ ylim(-1.2,3)+ geom_errorbar(ymin=prefSpecies1$lciPref1, ymax=prefSpecies1$uciPref1,size=0.2,width=0.3,col="black")+ geom_hline(yintercept = 0, linetype=1)+ geom_hline(yintercept = 0.51, linetype=2)+ geom_rect(aes(ymin=0.04, ymax=1.38, xmin=0, xmax=Inf), alpha=0.01)+ ggtitle("(a)")+ coord_flip()+ xlab("")+ ylab(expression(italic("Zr")[preference1]))+ theme_bw(base_size=11)+ theme(legend.position="none", plot.title=element_text(hjust=-0.8)) fp4<-ggplot(prefSpecies2, aes(reorder(common_study,-ZrPref2),ZrPref2))+ geom_point(aes(size=Neffects))+ ylim(-1.2,3)+ geom_errorbar(ymin=prefSpecies2$lciPref2, ymax=prefSpecies2$uciPref2,size=0.2,width=0.3,col="black")+ geom_hline(yintercept = 0, linetype=1)+ geom_hline(yintercept = 0.19, linetype=2)+ geom_rect(aes(ymin=-0.32, ymax=0.7, xmin=0, xmax=Inf), alpha=0.02)+ ggtitle("(b)")+ coord_flip()+ xlab("")+ ylab(expression(italic("Zr")[preference2]))+ theme_bw(base_size=11)+ theme(legend.position="none", plot.title=element_text(hjust=-0.8)) grid.arrange(fp3,fp4, ncol=1, nrow=2) ## Figure 4: ZrBrood size forest plot ## broodSpecies$common<-gsub("_", " ", broodSpecies$common) broodSpecies$Numstudy<-paste("(", broodSpecies$Nstudy, ")", sep="") broodSpecies<-transform(broodSpecies, common_study=paste(common, Numstudy, sep=" ")) broodEffects<-group_by(broodData, species) broodEffects<-summarise(broodEffects, Neffects=length(ZrBrood)) broodSpecies$Neffects<-broodEffects$Neffects ggplot(broodSpecies, aes(reorder(common_study,-ZrBrood),ZrBrood))+ geom_point(aes(size=Neffects),col="Black")+ ylim(-1,2)+ geom_hline(yintercept=0.47, linetype=2)+ geom_errorbar(ymin=broodSpecies$lciBrood, ymax=broodSpecies$uciBrood,size=0.2,width=0.3,col="black")+ geom_hline(yintercept = 0, linetype=1)+ geom_rect(aes(ymin=-0.19, ymax=0.86, xmin=0, xmax=Inf), alpha=0.01)+ coord_flip()+ xlab("")+ ylab(expression(italic("Zr")[brood~size]))+ theme_bw(base_size=11)+ theme(legend.position="none") ## Figure 5: ZrPreference ~ ZrBrood size scatterplot ## require("ggrepel") set.seed(20) PrefBrood$common<-gsub("_", " ", PrefBrood$common) ggplot(data=PrefBrood,aes(x=ZrBrood, y=ZrPref1))+ geom_smooth(method="lm",se=T,col="black",size=0.8)+ geom_point(stat="identity", size=3.5, col="black")+ ggtitle("")+ geom_text_repel(aes(label=common), force=4)+ xlim(0.25,1.25)+ xlab(expression(italic("Zr")[brood~size]))+ ylab(expression(italic("Zr")[preference1]))+ theme(axis.title=element_text(size=16))+ theme_bw(base_size=16) ## Figure S3: Funnel plots ## par(mfrow=c(3,2), oma=c(3,3,2,2), mar=c(4,4,0,0)) plot(ZrCond~Ncond, data=condStudy, xlab="", ylab="") mtext("A", adj=0) abline(h=0, lty=2) abline(h=-0.28, lty=1) plot(ZrOff~Noff, data=offStudy, xlab="", ylab="") mtext("B", adj=0) abline(h=0, lty=2) abline(h=0.25, lty=1) plot(ZrPref1~Npref1, data=prefStudy1, xlab="", ylab="") mtext("C", adj=0) abline(h=0, lty=2) abline(h=0.56, lty=1) plot(ZrPref2~Npref2, data=prefStudy2, xlab="", ylab="") mtext("D", adj=0) abline(h=0, lty=2) abline(h=0.19, lty=1) plot(ZrBrood~Nbrood, data=broodStudy, xlab="", ylab="") mtext("E", adj=0) abline(h=0, lty=2) abline(h=0.47, lty=1) mtext("Sample size", outer=T, side=1) mtext("Standardized effect size (Zr)", outer=T, side=2) # Figure S4: Brooding environments # condSpecies$brooding_environment<-gsub("_", "\n", condSpecies$brooding_environment) offSpecies$brooding_environment<-gsub("_", "\n", offSpecies$brooding_environment) require(gridExtra) plot1<-ggplot(condSpecies, aes(x=brooding_environment, y=ZrCond))+ geom_boxplot()+ geom_point(size=1, aes(col=brooding_environment), show.legend = F)+ labs(x="Brooding Environment", y="Zr Condition")+ ylab(expression(italic("Zr")[Condition]))+ ggtitle("A")+ theme_bw(base_size=10) plot3<-ggplot(offSpecies, aes(x=brooding_environment, y=ZrOff))+ geom_boxplot()+ geom_point(size=1, aes(colour=brooding_environment), show.legend = F)+ labs(y="Zr Offspring", x="Brooding Environment")+ ylab(expression(italic("Zr")[Offspring]))+ ggtitle("B")+ theme_bw(base_size=10) grid.arrange(plot1, plot3, ncol=2, nrow=1) # Figure S5: Laboratory and field studies # plot1<-ggplot(condData, aes(x=labfieldCon, y=ZrCond))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=labfieldCon, col=labfieldCon), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Condition", x="")+ ggtitle("A")+ theme_bw(base_size=10) plot2<-ggplot(offData, aes(x=labfieldOff, y=ZrOff))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=labfieldOff, col=labfieldOff), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Offspring", x="")+ ggtitle("B")+ theme_bw(base_size=10) plot3<-ggplot(prefData1, aes(x=labfieldPref1, y=ZrPref1))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=labfieldPref1, col=labfieldPref1), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference1", x="")+ ggtitle("C")+ theme_bw(base_size=10) plot4<-ggplot(prefData2, aes(x=labfieldPref2, y=ZrPref2))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=labfieldPref2, col=labfieldPref2), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference2", x="")+ ggtitle("D")+ theme_bw(base_size=10) plot5<-ggplot(broodData, aes(x=labfieldBrood, y=ZrBrood))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=labfieldBrood, col=labfieldBrood), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Brood Size", x="")+ ggtitle("E")+ theme_bw(base_size=10) grid.arrange(plot1, plot2, plot3, plot4, plot5, nrow=2) # Figure S6: Observational and experimental studies # plot1<-ggplot(condData, aes(x=obsexptCond, y=ZrCond))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=obsexptCond, col=obsexptCond), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Condition", x="")+ ggtitle("A")+ theme_bw(base_size=10) plot2<-ggplot(offData, aes(x=obsexptOff, y=ZrOff))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=obsexptOff, col=obsexptOff), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Offspring", x="")+ ggtitle("B")+ theme_bw(base_size=10) plot3<-ggplot(prefData1, aes(x=obsexptPref1, y=ZrPref1))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=obsexptPref1, col=obsexptPref1), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference1", x="")+ ggtitle("C")+ theme_bw(base_size=10) plot4<-ggplot(prefData2, aes(x=obsexptPref2, y=ZrPref2))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=obsexptPref2, col=obsexptPref2), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference2", x="")+ ggtitle("D")+ theme_bw(base_size=10) plot5<-ggplot(broodData, aes(x=obsexptBrood, y=ZrBrood))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=obsexptBrood, col=obsexptBrood), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Brood Size", x="")+ ggtitle("E")+ theme_bw(base_size=10) grid.arrange(plot1, plot2, plot3, plot4, plot5, nrow=2) # Figure S7: Continuous and binary care studies # plot1<-ggplot(condData, aes(x=contbinCond, y=ZrCond))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=contbinCond, col=contbinCond), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Condition", x="")+ ggtitle("A")+ theme_bw(base_size=10) plot2<-ggplot(offData, aes(x=contbinOff, y=ZrOff))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=contbinOff, col=contbinOff), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Offspring", x="")+ ggtitle("B")+ theme_bw(base_size=10) plot3<-ggplot(broodData, aes(x=contbinBrood, y=ZrBrood))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=contbinBrood, col=contbinBrood), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Brood Size", x="")+ ggtitle("C")+ theme_bw(base_size=10) grid.arrange(plot1, plot2, plot3, nrow=1) # Figure S8: Natural brood size and manipulated brood size studies # plot2<-ggplot(prefData1, aes(x=manipnatPref1, y=ZrPref1))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=manipnatPref1, col=manipnatPref1), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference1", x="")+ ggtitle("A")+ theme_bw(base_size=10) plot3<-ggplot(prefData2, aes(x=manipnatPref2, y=ZrPref2))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=manipnatPref2, col=manipnatPref2), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference2", x="")+ ggtitle("B")+ theme_bw(base_size=10) plot4<-ggplot(broodData, aes(x=manipnatBrood, y=ZrBrood))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=manipnatBrood, col=manipnatBrood), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Brood Size", x="")+ ggtitle("C")+ theme_bw(base_size=10) grid.arrange(plot2, plot3, plot4, nrow=1) # Figure S9: Natural brood absence and brood removal studies # plot1<-ggplot(prefData1, aes(x=eggremovePref1, y=ZrPref1))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=eggremovePref1, col=eggremovePref1), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference1", x="Brood removal")+ theme_bw(base_size=10) # Figure S10: Zr condition: condition measures # condData$condcatCond<-gsub("_", " ", condData$condcatCond) ggplot(condData, aes(x=condcatCond, y=ZrCond))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=condcatCond, col=condcatCond), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Condition", x="Condition measure")+ theme_bw(base_size=10) # Figure S11: Zr condition: care measures # condData$carecatCond<-gsub("_", " ", condData$carecatCond) ggplot(condData, aes(x=carecatCond, y=ZrCond))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=carecatCond, col=carecatCond), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Condition", x="Care measure")+ scale_x_discrete(labels=c("brood vs. no brood", "care period", "defence", "fanning", "number of broods"))+ theme_bw(base_size=10) # Figure S12: Zr offspring: offspring survival measures # offData$offcatoff<-gsub("_", " ", offData$offcatOff) ggplot(offData, aes(x=offcatOff, y=ZrOff))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=offcatOff, col=offcatOff), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr offspring", x="Offspring survival measure")+ scale_x_discrete(labels=c("number of eggs", "probability of clutch surviving", "proportion of eggs", "fanning", "number of broods"))+ theme_bw(base_size=10) # Figure S13: Zr offspring: care measures # offData$carecatOff<-gsub("_", " ", offData$carecatOff) ggplot(offData, aes(x=carecatOff, y=ZrOff))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=carecatOff, col=carecatOff), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr offspring", x="Care measure")+ scale_x_discrete(labels=c("nest attendance", "defence", "care duration", "fanning", "with vs without male"))+ theme_bw(base_size=10) # Figure S14: Zr preference1: female preference measures # prefData1$prefcatPref1<-gsub("_", " ", prefData1$prefcatPref1) ggplot(prefData1, aes(x=prefcatPref1, y=ZrPref1))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=prefcatPref1, col=prefcatPref1), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference1", x="Preference measure")+ scale_x_discrete(labels=c("association preference", "brood size", "spawning choice"))+ theme_bw(base_size=10) # Figure S15: Zr preference2: female preference measures # prefData2$prefcatPref2<-gsub("_", " ", prefData1$prefcatPref2) ggplot(prefData2, aes(x=prefcatPref2, y=ZrPref2))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=prefcatPref2, col=prefcatPref2), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Preference2", x="Preference measure")+ scale_x_discrete(labels=c("association preference", "brood size", "spawning choice"))+ theme_bw(base_size=10) # Figure S16: Zr brood size: brood size measures # broodData$sizecatBrood<-gsub("_", " ", broodData$sizecatBrood) ggplot(broodData, aes(x=sizecatBrood, y=ZrBrood))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=sizecatBrood, col=sizecatBrood), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Brood size", x="Brood size measure")+ scale_x_discrete(labels=c("brood area", "absolute egg number", "egg number category"))+ theme_bw(base_size=10) # Figure S17: Zr brood size: care measures # broodData$carecatBrood<-gsub("_", " ", broodData$carecatBrood) ggplot(broodData, aes(x=carecatBrood, y=ZrBrood))+ geom_boxplot(width=0.5)+ geom_point(size=1, aes(fill=carecatBrood, col=carecatBrood), position = position_jitterdodge(), show.legend=F)+ labs(y="Zr Brood size", x="Care measure")+ scale_x_discrete(labels=c("brood cared for\nvs. abandoned", "brood cared for\nvs. cannibalised", "defence", "fanning", "other care category"))+ theme_bw(base_size=10) # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ # ### PART 8: PHYLOGENETIC TREE ### ## import tree from Betancur-R et al. 2017. (BMC Evolutionary Biology) as a .nexus file ## is.ultrametric(tree2) # Remove text after third '_' so just left with taxonomic names tree2$tip.label <- gsub("^([^_]*_[^_]*_[^_]*)_.*$", "\\1", tree2$tip.label) # Remove text before second '_' so just left with genus and species tree2$tip.label <- unlist(regmatches(tree2$tip.label, gregexpr("(?<=_).*", tree2$tip.label, perl=TRUE))) ## replace missing species with species in same genus ## tree2$tip.label[which(tree2$tip.label == "Acanthemblemaria_aspera")] <- "Acanthemblemaria_crockeri" # replace Acanthemblemaria_aspera tree2$tip.label[which(tree2$tip.label == "Ostorhinchus_cookii")] <- "Apogon_doederleini" # replace Ostorhinchus_cookii tree2$tip.label[which(tree2$tip.label == "Ostorhinchus_lateralis")] <- "Apogon_notatus" # replace Ostorhinchus_lateralis tree2$tip.label[which(tree2$tip.label == "Chromis_atripectoralis")] <- "Chromis_notata" # replace Chromis_atripectoralis tree2$tip.label[which(tree2$tip.label == "Chrysiptera_taupou")] <- "Chrysiptera_cyanea" # replace Chrysiptera_taupou tree2$tip.label[which(tree2$tip.label == "Etheostoma_juliae")] <- "Etheostoma_crossopterum" # replace Etheostoma_juliae tree2$tip.label[which(tree2$tip.label == "Etheostoma_simoterum")] <- "Etheostoma_flabellare" # replace Etheostoma_simoterum tree2$tip.label[which(tree2$tip.label == "Etheostoma_vitreum")] <- "Etheostoma_olmstedi" # replace Etheostoma_vitreum tree2$tip.label[which(tree2$tip.label == "Eviota_albolineata")] <- "Eviota_prasina" # replace Eviota_albolineata tree2$tip.label[which(tree2$tip.label == "Lepomis_cyanellus")] <- "Lepomis_gibbosus" # replace Lepomis_cyanellus tree2$tip.label[which(tree2$tip.label == "Pimephales_vigilax")] <- "Pimephales_promelas" # replace Pimephales_vigilax tree2$tip.label[which(tree2$tip.label == "Stegastes_diencaeus")] <- "Stegastes_adustus" # replace Stegastes_diencaeus tree2$tip.label[which(tree2$tip.label == "Stegastes_fuscus")] <- "Stegastes_rectifraenum" # replace Stegastes_fuscus tree2$tip.label[which(tree2$tip.label == "Cottus_carolinae")] <- "Cottus_gobio" # replace Cottus_carolinae ## replace missing species with species in same family ## tree2$tip.label[which(tree2$tip.label == "Hypleurochilus_sp")] <- "Aidablennius_sphynx" # replace Hypleurochilus_sp tree2$tip.label[which(tree2$tip.label == "Lepidogobius_lepidus")] <- "Gobiusculus_flavescens" # replace Coryphopterus_glaucofraenum tree2$tip.label[which(tree2$tip.label == "Batrachoides_pacifici")] <- "Halobatrachus_didactylus" # replace Batrachoides_pacifici tree2$tip.label[which(tree2$tip.label == "Pleurogrammus_monopterygius")] <- "Oxylebius_pictus" # replace Pleurogrammus_monopterygius tree2$tip.label[which(tree2$tip.label == "Blenniella_chrysospilos")] <- "Rhabdoblennius_nitidus" # replace Blenniella_chrysospilos tree2$tip.label[which(tree2$tip.label == "Lutjanus_griseus")] <- "Symphodus_tinca" # replace Lutjanus_griseus tree2$tip.label[which(tree2$tip.label == "Hypsoblennius_hentz")] <- "Salaria_fluviatilis" # replace Hypsoblennius_hentz tree2$tip.label[which(tree2$tip.label == "Leptocottus_armatus")] <- "Artedius_fenestralis" # replace Leptocottus_armatus ## trim the tree ## dropTip <- tree2$tip.label[which(is.na(match(tree2$tip.label, zrData$species)))] # 39 species at this stage zrTree <- drop.tip(tree2, dropTip, trim.internal=T) zrTree <- makeNodeLabel(zrTree, method = "number") ## add missing species to existing genera ## zrTree2 <- add.species.to.genus(zrTree, "Micropterus_dolomieu") zrTree2 <- add.species.to.genus(zrTree2, "Salaria_pavo") zrTree2 <- add.species.to.genus(zrTree2, "Cottus_nozawae") zrTree2 <- add.species.to.genus(zrTree2, "Artedius_harringtoni") # 43 species at this stage ## re-label the nodes ## zrTree2 <- makeNodeLabel(zrTree2, method = "number") plot(zrTree2, show.node.label=TRUE) ## add Rhinogobius, half way up Gobiusculus branch which(zrTree$tip.label=="Gobiusculus_flavescens") # 34 (edge number of this tip) which(zrTree2$edge[,2] == 34) # 70 (row with this edge length) zrTree2$edge.length[70] # 22.29 zrTree2 <- bind.tip(zrTree2, "Rhinogobius_brunneus", edge.length=NULL, where=38, position=zrTree2$edge.length[70]/2) ## re-label the nodes ## zrTree2 <- makeNodeLabel(zrTree2, method = "number") plot(zrTree2, show.node.label=TRUE) ## add Gymnogobius, half way up Rhinogobius branch which(zrTree2$tip.label=="Rhinogobius_brunneus") # 39 (edge number of this tip) which(zrTree2$edge[,2] == 39) # 79 (row with this edge length) zrTree2$edge.length[79] # 11.14 zrTree2 <- bind.tip(zrTree2, "Gymnogobius_isaza", edge.length=NULL, where=39, position=zrTree2$edge.length[79]/2) ## re-label the nodes ## zrTree2 <- makeNodeLabel(zrTree2, method = "number") plot(zrTree2, show.node.label=TRUE) ## add Pomatoschitus, half way up Gobiusculus branch which(zrTree2$tip.label=="Gobiusculus_flavescens") # 38 (edge number of this tip) which(zrTree2$edge[,2] == 38) # 78 (row with this edge length) zrTree2$edge.length[78] # 11.14 zrTree2 <- bind.tip(zrTree2, "Pomatoschistus_minutus", edge.length=NULL, where=38, position=zrTree2$edge.length[78]/2) ## re-label the nodes ## zrTree2 <- makeNodeLabel(zrTree2, method = "number") plot(zrTree2, show.node.label=TRUE) ## add missing species to newly created genera ## zrTree2 <- add.species.to.genus(zrTree2, "Rhinogobius_flumineus") zrTree2 <- add.species.to.genus(zrTree2, "Pomatoschistus_microps") ## re-label the nodes ## zrTree2 <- makeNodeLabel(zrTree2, method = "number") plot(zrTree2, show.node.label=TRUE) # ------------------------------------------------------------------------------------------------------ # # ------------------------------------------------------------------------------------------------------ #