# Meta-analysis of fledging success for nest predation dataset # set working directory # read in data x <- read.csv("nest_predation_dataset.csv", stringsAsFactors = FALSE) nrow(x) # [1] 515 # rename column for analysis x$ApparentSuccess <- x$PercentageSuccessfulNests # Extra variables x$nApparentSuccess <- round(x$SampleSize * x$ApparentSuccess) # approx number of successes # Did a fudge here to avoid zero variance when success = 0 # Set a minimum success of 0.01 range(x$ApparentSuccess, na.rm = TRUE) # [1] 0.0000 0.9286 x$ApparentSuccess[x$ApparentSuccess < 0.01] <- 0.01 # x$ApparentSuccess[x$ApparentSuccess < 0.05] <- 0.05 # ** results virtually identical as for 0.01 # Standard error of proportion x$SampleSize # no NA's, good x$ASse <- sqrt(x$ApparentSuccess * (1 - x$ApparentSuccess) / x$SampleSize) head(x[ , c("ApparentSuccess", "SampleSize")]) # ApparentSuccess SampleSize # 1 0.6790197 125 # 2 0.2675000 243 # 3 0.2900000 797 # 4 0.4956000 900 # 5 0.5220000 306 # 6 0.5750000 675 # drop rows without apparent success data x <- x[!is.na(x$ApparentSuccess), ] dim(x) # 369 # -------------------------------------- # Analyses of means using metafor package # Default method is random effects using REML library(metafor) library(ggplot2) theme_set(theme_classic()) # this works better here than setting "+theme_classic()" on each command, which causes +theme() not to work # Plot the data # loess in ggplot can be weighted by exected variance, calculated for random effects meta-analysis as 1/(SE^2 + Tau^2) # Here I used Tau^2 = 0.03, which is roughly the value in model fits below. #pdf("apparent success - loess fit.pdf") ggplot(x, aes(x = Lat, y = ApparentSuccess, weight = 1/(ASse^2 + 0.03))) + geom_point(size = 2, col = "firebrick", alpha = 0.5) + geom_smooth(method = loess, size = 1) + labs(x = "Latitude", y = "Apparent success") + theme(aspect.ratio = 0.80) #dev.off() # --- # Fit a single mean apparent success rate z <- rma.uni(ApparentSuccess ~ 1, vi = ASse^2, data = x) summary(z) AIC(z) # -168.1847 x$predict <- predict(z)$pred #pdf("apparentsuccess_interceptonly.pdf") ggplot(x, aes(x = Lat, y = ApparentSuccess, weight = 1/(ASse^2 + z$tau2))) + geom_point(size = 2, col = "firebrick", alpha = 0.5) + geom_smooth(method = loess, size = 1) + geom_line( aes( x = Lat, y = predict), linetype = "dashed", size = 1) + labs(x = "Latitude", y = "Apparent success") + theme(aspect.ratio = 0.80) #dev.off() # terrible fit # --- # Linear regression model with latitude - equal slopes # Fit is a pair of straight lines outward from the equator having the same slope z <- rma.uni(ApparentSuccess ~ abs(Lat), vi = ASse^2, data = x) x$predict <- predict(z)$pred AIC (z) # -213.2483 summary(z) # Using the Tau^2 from the parametric fit (it is always about 0.03) #pdf("apparentsuccess_linearregression.pdf") ggplot(x, aes(x = Lat, y = ApparentSuccess, weight = 1/(ASse^2 + z$tau2))) + geom_point(size = 2, col = "firebrick", alpha = 0.5) + geom_smooth(method = loess, size = 1) + geom_line( aes( x = Lat, y = predict), linetype = "dashed", size = 1) + labs(x = "Latitude", y = "Apparent success") + theme(aspect.ratio = 0.80) #dev.off() # bad fit, especially in S hemisphere # --- # Breakpoint regression - equal slopes around 23.4 deg latitude z <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat) > 23.4) * (abs(Lat) - 23.4)), vi = ASse^2, data = x) x$predict <- predict(z)$pred AIC(z) # -229.8574 summary(z) ggplot(x, aes(x = Lat, y = ApparentSuccess, weight = 1/(ASse^2 + z$tau2))) + geom_point(size = 2, col = "firebrick", alpha = 0.5) + geom_smooth(method = loess, size = 1) + geom_line( aes( x = Lat, y = predict), linetype = "dashed", size = 1) + labs(x = "Latitude", y = "Apparent success") + theme(aspect.ratio = 0.80) # Cleaner lines using a grid of xpoints xpts <- seq(min(x$Lat), max(x$Lat), by = 0.1) yhat <- coef(z)[1] + I( coef(z)[2]*as.integer(abs(xpts) > 23.4)*(abs(xpts) - 23.4) ) pdf("apparentsuccess_breakpoint.pdf", height=5, width=5) ggplot() + geom_point(data = x, aes(x = Lat, y = ApparentSuccess), size = 2, col = "firebrick", alpha = 0.5) + geom_smooth(data = x, aes(x = Lat, y = ApparentSuccess), method = loess, size = 1) + geom_line(data = data.frame(xpts, yhat), mapping = aes( x = xpts, y = yhat), linetype = "dashed", color = "black", size = 1) + labs(x = "Latitude", y = "Fledging success") + theme(aspect.ratio = 0.80) dev.off() # good fit except for S hemisphere, where data looks flat # so fit an additional model that lets slope vary between S & N Hemisphere because data look different for S Hemisphere z <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat)>23.4)*(abs(Lat)-23.4)) + I(as.integer(Lat>23.4)*(Lat-23.4)), vi = ASse^2, data = x) x$predict <- predict(z)$pred AIC(z) # -234.536 summary(z) ggplot(x, aes(x = Lat, y = ApparentSuccess, weight = 1/(ASse^2 + z$tau2))) + geom_point(size = 2, col = "firebrick", alpha = 0.5) + geom_smooth(method = loess, size = 1) + geom_line( aes( x = Lat, y = predict), linetype = "dashed", size = 1) + labs(x = "Latitude", y = "Apparent success") + theme(aspect.ratio = 0.80) # Cleaner lines using a grid of x-points pdf("apparentsuccess_breakpoint_unequal.pdf", width = 5, height = 5) xpts <- seq(min(x$Lat), max(x$Lat), by = 0.1) yhat <- coef(z)[1] + I(coef(z)[2] * as.integer(abs(xpts) > 23.4)*(abs(xpts) - 23.4)) + I(coef(z)[3]*as.integer(xpts > 23.4)*(xpts - 23.4)) ggplot() + geom_point(data = x, aes(x = Lat, y = ApparentSuccess), size = 2, col = "firebrick", alpha = 0.5) + geom_smooth(data = x, aes(x = Lat, y = ApparentSuccess), method = loess, size = 1) + geom_line(data = data.frame(xpts, yhat), mapping = aes( x = xpts, y = yhat), linetype = "dashed", color = "black", size = 1) + labs(x = "latitude", y = "fledging success") + theme(aspect.ratio = 0.80) dev.off() # reviewers suggested testing impact of covariates z1 <- rma.uni(ApparentSuccess ~ 1 + demelevation, vi = ASse^2, data = x) z2 <- rma.uni(ApparentSuccess ~ 1 + NestShape, vi = ASse^2, data = x) z3 <- rma.uni(ApparentSuccess ~ 1 + Habitat, vi = ASse^2, data = x) z4 <- rma.uni(ApparentSuccess ~ abs(Lat) + demelevation, vi = ASse^2, data = x) z5 <- rma.uni(ApparentSuccess ~ abs(Lat) + NestShape, vi = ASse^2, data = x) z6 <- rma.uni(ApparentSuccess ~ abs(Lat) + Habitat, vi = ASse^2, data = x) z7 <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat) > 23.4) * (abs(Lat) - 23.4)) + demelevation, vi = ASse^2, data = x) z8 <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat) > 23.4) * (abs(Lat) - 23.4)) + NestShape, vi = ASse^2, data = x) z9 <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat) > 23.4) * (abs(Lat) - 23.4)) + Habitat, vi = ASse^2, data = x) z10 <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat)>23.4)*(abs(Lat)-23.4)) + I(as.integer(Lat>23.4)*(Lat-23.4)) + demelevation, vi = ASse^2, data = x) z11 <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat)>23.4)*(abs(Lat)-23.4)) + I(as.integer(Lat>23.4)*(Lat-23.4)) + NestShape, vi = ASse^2, data = x) z12 <- rma.uni(ApparentSuccess ~ 1 + I(as.integer(abs(Lat)>23.4)*(abs(Lat)-23.4)) + I(as.integer(Lat>23.4)*(Lat-23.4)) + Habitat, vi = ASse^2, data = x) AIC(z1, z2, z3, z4, z5, z6, z7, z8, z9, z10, z11, z12) #df AIC #z1 3 -165.0989 #z2 3 -168.0928 #z3 4 -168.5993 #z4 4 -209.6706 #z5 4 -210.3158 #z6 5 -207.7217 #z7 4 -226.2698 #z8 4 -226.4228 #z9 5 -223.3509 #z10 5 -230.9176 #z11 5 -230.9139 #z12 6 -228.9453 ############ # same analysis but incorporating phylogeny and species ID ############### # Phylogenetic VCV library(ape) library(phytools) library(metafor) tree <- read.nexus("mcc_pruned_tree_app_success.tre") branches <- compute.brlen(tree) cor <- vcv(branches, cor = T) # run for equal slopes breakpoint model phy.mod.re <- rma.mv(yi = ApparentSuccess, mods=c(I(as.integer(abs(Lat) > 23.4) * (abs(Lat) - 23.4))), V = ASse^2, random = list(~1|JetzTipLabel, ~1|Citation, ~1|species_id), R = list(JetzTipLabel = cor), method = "REML", data = x) summary(phy.mod.re) # and also run this for the unequal slopes breakpoint model phy.mod.re <- rma.mv(ApparentSuccess ~ I(as.integer(abs(Lat)>23.4)*(abs(Lat)-23.4)) + I(as.integer(Lat>23.4)*(Lat-23.4)), V = ASse^2, random = list(~1|JetzTipLabel, ~1|species_id), R = list(JetzTipLabel = cor), method = "REML", data = x) summary(phy.mod.re) # pretty much identical for both # Looks like none of the random effect are really all that important nrow(x) # 369 jetz1<-read.nexus("predation_trees.tre") # prune tree to 369 species in final data set. Latin names must match those in the tre file sp<-x$JetzTipLabel pruned.tree<-drop.tip(jetz1[[1]],jetz1[[1]]$tip.label[match(sp, jetz1[[1]]$tip.label)]) other_names<-pruned.tree$tip.label # repeat to prune all the phylogenies in the list final.tree<-lapply(jetz1,drop.tip,tip=other_names) class(final.tree) tree<-final.tree[[1]] write.nexus(tree, file = "pruned_trees_app_success.tre") # same number of tips as in nomenclature? length(unique(x$JetzTipLabel)) length(unique(tree$tip.label)) #Calculate branch lengths and set as VCV matrix for rma.mv branches <- compute.brlen(tree) cor <- vcv(branches, cor = T) # Run rma models # Random effects model only # Include phylogeny as the correlation matrix while accounting for multiple species in some studies # run for equal slopes breakpoint model phy.mod.re <- rma.mv(yi = ApparentSuccess, mods=c(I(as.integer(abs(Lat) > 23.4) * (abs(Lat) - 23.4))), V = ASse^2, random = list(~1|JetzTipLabel, ~1|Citation, ~1|species_id), R = list(JetzTipLabel = cor), method = "REML", data = x) summary(phy.mod.re) # and also run this for the unequal slopes breakpoint model phy.mod.re <- rma.mv(ApparentSuccess ~ I(as.integer(abs(Lat)>23.4)*(abs(Lat)-23.4)) + I(as.integer(Lat>23.4)*(Lat-23.4)), V = ASse^2, random = list(~1|JetzTipLabel, ~1|species_id), R = list(JetzTipLabel = cor), method = "REML", data = x) summary(phy.mod.re) # pretty much identical for both # Looks like none of the random effect are really all that important! Proceed with rma.uni?