--- title: "Predator.Parasite.MetaAnalysis" author: "D Daversa" date: "16/11/2017" output: word_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = F) ``` ##This code was started on 11 July 2017 by D Daversa at the University of Liverpool ##Project: "Broadening the Ecology of Fear: Non-lethal effects arise from diverse responses to predators and parasites" ## Description: a review and analysis of studies that have quantified trait changes elicited by predators and parasites in the same resource species and set of experiments **DATA PREPARATION** Here we will upload the full dataset from studies meeting the criteria of our review. We will do a bit of organization of the columns, and remove studies that we had to exclude due to data quality or accessibility issues. We will then calculate the standardized mean difference using the metafor package as a check of our personal calculations in the main dataset. Finally, we will make a data subset for studies using tadpoles as the focal prey/hosts for further analysis. ```{r, include=FALSE} ##read in database of predator + parasite studies read.csv("/Users/ddaversa/Desktop/R/predator_parasite_review/NLE_review_rawdata.csv")->raw.data ##load the data raw.data<-droplevels(raw.data[which(raw.data$study_id != ""),]) ##remove excluded studies and entries raw.data<-droplevels(raw.data[which(raw.data$outcome == "included"),]) ##turn entry ID into a factor raw.data$entry_id<-factor(as.character(raw.data$entry_id)) ##calculate effect size and variance using metafor functions library(metafor) raw.data<-escalc(n1i = n_treat, n2i = n_control, m1i = treat_mean, m2i = control_mean, sd1i = treat_sd, sd2i = control_sd, data = raw.data, measure = "SMD", append = TRUE) raw.data$yi<-as.numeric(as.character(raw.data$yi)) raw.data$yi<--(raw.data$yi) ##reorder interaction_type levels raw.data$interaction_type<-factor(raw.data$interaction_type, levels = c("predator-prey", "host-parasite", "combined")) raw.data<-raw.data[which(raw.data$trait_gen != "stress"),] ##compute adjusted sampling variance raw.data$vi_2<-with(raw.data, 1/n_treat + 1/n_control + yi^2/(2*n_tot)) ##reverse signs of studies whose predicted effect are positive values (e.g. refuge use, distance to cues) for (x in 1:nrow(raw.data)) { if (raw.data$predicted_effect[x] == "negative") { raw.data$hedges_d[x]<--raw.data$hedges_d[x] raw.data$yi[x]<--raw.data$yi[x] } } ##Calculate confidence intervals of each study raw.data$ci.lower<-raw.data$yi-sqrt(raw.data$vi_2)*1.96 raw.data$ci.upper<-raw.data$yi+sqrt(raw.data$vi_2)*1.96 ##create classification of positive, negative, or no observed effects raw.data$observed_effect<-"none" for (i in 1:nrow(raw.data)) { if (raw.data$ci.lower[i] > 0) { if (raw.data$ci.upper[i] > 0) { raw.data$observed_effect[i]<-"positive" } } if (raw.data$ci.lower[i] < 0) { if (raw.data$ci.upper[i] < 0) { raw.data$observed_effect[i]<-"negative" } } } ##create subset of amphibian resources raw.data.sub<-droplevels(subset(raw.data, resource_group == "amphibian")) ##reorder consumer state levels raw.data.sub$consumer_state<-factor(raw.data.sub$consumer_state, levels = c("questing", "attacking", "consuming")) ##create subsets parasite.data<-droplevels(subset(raw.data.sub, interaction_type == "host-parasite")) predator.data<-droplevels(subset(raw.data.sub, interaction_type == "predator-prey")) combined.data<-droplevels(subset(raw.data.sub, interaction_type == "combined")) ``` With our dataframe prepared, we will perform a series of steps that use meta-analytic approaches to test the influence of different factors on the strength of trait responses. First, we will do a bit more organization of the dataset, and calculate co-variance matrices for the madel that account for multi-arm experimental designs. ```{r } library(lme4) library(MASS) library(clubSandwich) ##drop stress measures - not enough samples raw.data.sub2<-droplevels(subset(raw.data.sub, trait_gen_b != "stress")) ##relevel trait levels raw.data.sub2$trait_gen_b<-factor(raw.data.sub2$trait_gen_b, levels = c("activity", "space use","morphology/physiology")) ##create co-variance matrices ##calculate covariation matrix to account for same control groups for multiple treatments calc.v <- function(x) { v <- matrix(1/x$n_treat[1] + outer(x$yi, x$yi, "*")/(2*x$n_tot[1]), nrow=nrow(x), ncol=nrow(x)) diag(v) <- x$vi_2 v } ##calculate covariation matrix to account for same control groups for multiple treatments V<-bldiag(lapply(split(raw.data, list(raw.data$study_id, raw.data$grouping_factor), drop = T), calc.v)) V.sub <- bldiag(lapply(split(raw.data.sub, list(raw.data.sub$study_id, raw.data.sub$grouping_factor), drop = T), calc.v)) V.para<- bldiag(lapply(split(parasite.data, list(parasite.data$study_id, parasite.data$grouping_factor), drop = T), calc.v)) V.pred<- bldiag(lapply(split(predator.data, list(predator.data$study_id, predator.data$grouping_factor), drop = T), calc.v)) V.combined<- bldiag(lapply(split(combined.data, list(combined.data$study_id, combined.data$grouping_factor), drop = T), calc.v)) V.sub2<- bldiag(lapply(split(raw.data.sub2, list(raw.data.sub2$study_id, raw.data.sub2$grouping_factor), drop = T), calc.v)) ``` We will now do preliminary tests to determine whether or not to include interactions between different factors. ```{r} ##Full model lmm1<-rma.mv(yi, V.sub2, mods = ~interaction_type + analysis_scale + trait_gen_b + interaction_type:analysis_scale + interaction_type:trait_gen_b, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm1, digits = 3) ##________________________________ ## TEST FOR INTERACTIONS ## ##________________________________ ##drop analysis scale:interaction type lmm1b<-rma.mv(yi, V.sub2, mods = ~interaction_type + analysis_scale + trait_gen_b + interaction_type:trait_gen_b, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm1b, digits = 3) anova(lmm1, lmm1b, test = "Chisq") ##drop trait:interaction type lmm1c<-rma.mv(yi, V.sub2, mods = ~interaction_type + analysis_scale + trait_gen_b + interaction_type:analysis_scale, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm1c, digits = 3) anova(lmm1, lmm1c, test = "Chisq") ``` Here, we update the base model to use in our analyses, dropping the non-significant interaction ```{r} ##Update model with all interactions dropped lmm.sub<-rma.mv(yi, V.sub2, mods = ~interaction_type*analysis_scale + trait_gen_b, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm.sub, digits = 3) ``` Now, we will test the the main effect of interaction type on effect sizes on tadpole traits. This will tell us whether, on average, the strength of trait respones to parasites differ from the strength of trait responses to predators, and their combined presence. ```{r} ##likelihood ratio test of main effect of interaction type lmm.sub1a<-rma.mv(yi, V.sub2, mods = ~analysis_scale + trait_gen_b, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm.sub1a, digits = 3) ##Likelihood ratio test against based model, using X2 distribution anova(lmm.sub, lmm.sub1a, test = "Chisq") ##compare factor levels lmm.sub1b<-rma.mv(yi, V.sub2, mods = ~interaction_type, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm.sub1b, digits = 3) ##Estimate the contribution of factor levels to variance lmm.sub1c<-rma.mv(yi, V.sub2, mods = ~interaction_type-1, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "REML") print(lmm.sub1c, digits = 3) ## predators and combined presence elicit signficant responses ##Robust sandwich estimation lmm.sub.1d<-rma.mv(yi, vi, mods = ~interaction_type-1, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "REML") robust(lmm.sub.1d, cluster = raw.data$study_id, adjust = T) ``` Now, we will test the main effect of analysis scale and the type of trait measured on the strength of responeses. Analysis scale distinguishes between studies that measured responses at the group-level (e.g. proportion of individuals moving) from from that measured responses in individuals. ```{r} ##______________________________ ## ANALYSIS SCALE ##______________________________ ##likelihood ratio test of main effect lmm.sub.2a<-rma.mv(yi, V.sub2, mods = ~interaction_type + trait_gen_b, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") anova(lmm.sub, lmm.sub.2a, test = "Chisq") ##compare factor levels lmm.sub.2b<-rma.mv(yi, V.sub2, mods = ~analysis_scale, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm.sub.2b, digits = 3) ##Estimate the contribution of factor levels to variance lmm.sub.2c<-rma.mv(yi, V.sub2, mods = ~interaction_type:analysis_scale-1, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "REML") print(lmm.sub.2c, digits = 3) ##robust sandwich estimation lmm.sub.2d<-rma.mv(yi, vi, mods = ~interaction_type:analysis_scale-1, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "REML") robust(lmm.sub.2d, cluster = raw.data.sub2$study_id, adjust = T, digits = 3) ##______________________________ ## trait type ##______________________________ ##drop main effect of trait type lmm.sub.3a<-rma.mv(yi, V.sub2, mods = ~interaction_type + analysis_scale, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm.sub.3a, digits = 3) anova(lmm.sub, lmm.sub.3a, test = "Chisq") #compare factor levels lmm.sub.3b<-rma.mv(yi, V.sub2, mods = ~trait_gen_b, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "ML") print(lmm.sub.3b, digits = 3) ##Estimate the contribution of factor levels to variance lmm.sub.3c<-rma.mv(yi, V.sub2, mods = ~interaction_type:trait_gen_b-1, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "REML") print(lmm.sub.3c, digits = 3) #robust sandwich estimation lmm.sub.3d<-rma.mv(yi, vi, mods = ~interaction_type:trait_gen_b-1, random = ~1|study_id/entry_id, data = raw.data.sub2, method = "REML") robust(lmm.sub.3d, cluster = raw.data.sub2$study_id, adjust = T, digits = 3) ``` Here, we will repeat the above analysis using only data from questing states of predators or parasites. This controls for the predator-prey data being limited to questing predator states. ```{r} ##create subset of questing state questing.data<-droplevels(subset(raw.data.sub2, consumer_state == "questing")) ##calculate covariation matrix to account for same control groups for multiple treatments V.quest<- bldiag(lapply(split(questing.data, list(questing.data$study_id, questing.data$grouping_factor), drop = T), calc.v)) ##model without interaction type questing.null<-rma.mv(yi, V.quest, random = ~1|study_id/entry_id, data = questing.data, method = "ML") ##fit model with interaction type as fixed effect lmm.quest.1<-rma.mv(yi, V.quest, mods = ~interaction_type, random = ~1|study_id/entry_id, data = questing.data, method = "ML") print(lmm.quest.1, digits = 3) par(mfrow=c(2,1)) ##profile likelihood profile(lmm.quest.1, sigma2=1) profile(lmm.quest.1, sigma2=2) ## LRT of significance of interaction type anova(questing.null, lmm.quest.1) ##compare factor levels lmm.quest.1b<-rma.mv(yi, V.quest, mods = ~interaction_type, random = ~1|study_id/entry_id, data = questing.data, method = "ML") print(lmm.quest.1b, digits = 3) ##assess the influence of each factor level on effect size lmm.quest.1c<-rma.mv(yi, V.quest, mods = ~interaction_type-1, random = ~1|study_id/entry_id, data = questing.data, method = "ML") print(lmm.quest.1c, digits = 3) #robust sandwich estimation lmm.quest.1d<-rma.mv(yi, V.quest, mods = ~interaction_type-1, random = ~1|study_id/entry_id, data = questing.data, method = "ML") robust(lmm.quest.1d, cluster=questing.data$study_id) ``` Here, we will analyse the parasite data to determine whether the class of trait reponses - avoid, counter, combat - or parasite strategy as per Lafferty and Kuris (2002) influence the effect sizes on host traits ```{r } ##subset data parasite.sub<-droplevels(raw.data.sub2[which(raw.data.sub2$interaction_type == "host-parasite"),]) #full model lmm.psite.1<-rma.mv(yi, vi,mods = ~consumer_state + consumer_type, random = ~1|study_id/entry_id, data = parasite.sub, method = "ML") print(lmm.psite.1, digits = 3) ##__________________________ ## EFFECTS OF PARASITE STATE ##___________________________ ##model with parasite state dropped lmm.psite.1b<-rma.mv(yi, vi, mods = ~consumer_type, random = ~1|study_id/entry_id, data = parasite.sub, method = "ML") print(lmm.psite.1b, digits = 3) ##likihood ratio test anova(lmm.psite.1, lmm.psite.1b, test = "Chisq") ##assess the influence of each factor level on effect size lmm.psite.1c<-rma.mv(yi, vi, mods = ~consumer_state-1, random = ~1|study_id/entry_id, data = parasite.sub, method = "REML") print(lmm.psite.1c, digits = 3) ##___________________________ ##EFFECTS OF PARASITE STRATEGY ##___________________________ ##fit model with parasite type dropped lmm.psite.2<-rma.mv(yi, vi, mods = ~consumer_state, random = ~1|study_id/entry_id, data = parasite.sub, method = "ML") print(lmm.psite.2, digits = 3) ##likihood ratio test anova(lmm.psite.1, lmm.psite.2, test = "Chisq") ##assess the influence of each factor level on effect size lmm.psite.2b<-rma.mv(yi, vi, mods = ~consumer_type-1, random = ~1|study_id/entry_id, data = parasite.sub, method = "REML") print(lmm.psite.2b, digits = 3) ``` Compute the I2 statistic to assess between- vs. within-study heterogeneity ```{r} ##____________________________ ##(taken from metafor website analysis examples (http://www.metafor-project.org/doku.php/tips:i2_multilevel_multivariate) ##___________________________ W <- diag(1/raw.data.sub2$vi) X <- model.matrix(lmm.sub) P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W 100 * sum(lmm.sub$sigma2) / (sum(lmm.sub$sigma2) + (lmm.sub$k-lmm.sub$p)/sum(diag(P))) 100 * lmm.sub$sigma2 / (sum(lmm.sub$sigma2) + (lmm.sub$k-lmm.sub$p)/sum(diag(P))) ``` **FIGURES** *Figure 2* ```{r FOREST PLOTS } library(ggplot2) ##relevel trait type factor levels raw.data.sub2$trait_gen_b<-factor(raw.data.sub2$trait_gen_b, levels = c("activity", "space use", "morphology/physiology")) ##set color pallete forest.colors<-c("#6699CC", "#990000", "#333333") ##create subsets parasite.data<-droplevels(subset(raw.data.sub2, interaction_type == "host-parasite")) predator.data<-droplevels(subset(raw.data.sub2, interaction_type == "predator-prey")) combined.data<-droplevels(subset(raw.data.sub2, interaction_type == "combined")) ##PARASITES effects.means.parasite<-ggplot(parasite.data, aes(x = reorder(entry_id, yi), y = yi), group = trait_gen_b_b) effects.means.parasite<-effects.means.parasite + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.parasite<-effects.means.parasite + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.parasite<-effects.means.parasite + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.parasite<-effects.means.parasite + scale_colour_manual("trait class", values = forest.colors) effects.means.parasite<-effects.means.parasite + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 20) effects.means.parasite<-effects.means.parasite + theme(legend.position = "none",text = element_text(family = "Times", size = 20), axis.text=element_text(size=18), axis.text.y=element_blank(), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.parasite<-effects.means.parasite + scale_y_continuous(limits = c(-7,7), breaks=seq(-7, 7, 1)) effects.means.parasite<-effects.means.parasite + ggtitle("Parasites") effects.means.parasite ggsave(file="pred.v.parasite.Fig2 Forest Plot_parasite.pdf", plot=effects.means.parasite, width=7, height=10) ##PREDATORS effects.means.predator<-ggplot(predator.data, aes(reorder(entry_id, yi), yi), group = trait_gen_b_b) effects.means.predator<-effects.means.predator + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.predator<-effects.means.predator + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.predator<-effects.means.predator + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.predator<-effects.means.predator + scale_colour_manual("trait", values = forest.colors) effects.means.predator<-effects.means.predator + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 20) effects.means.predator<-effects.means.predator + theme( legend.position = "none",text = element_text(family = "Times", size = 20), axis.text=element_text(size=18), axis.text.y=element_blank(), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.predator<-effects.means.predator + scale_y_continuous(limits = c(-7,7), breaks=seq(-7, 7, 1)) effects.means.predator<-effects.means.predator + ggtitle("Predators") effects.means.predator ggsave(file="pred.v.predator.Fig2 Forest Plot_predator.pdf", plot=effects.means.predator, width=7, height=10) ##COMBINED effects.means.combined<-ggplot(combined.data, aes(reorder(entry_id, yi), yi), group = trait_gen_b_b) effects.means.combined<-effects.means.combined + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.combined<-effects.means.combined + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.combined<-effects.means.combined + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.combined<-effects.means.combined + scale_colour_manual("trait", values = forest.colors) effects.means.combined<-effects.means.combined + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 20) effects.means.combined<-effects.means.combined + theme( legend.position = "none", text = element_text(family = "Times", size = 20), axis.text=element_text(size=18), axis.text.y=element_blank(), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.combined<-effects.means.combined + scale_y_continuous(limits = c(-7,7), breaks=seq(-7, 7, 1)) effects.means.combined<-effects.means.combined + ggtitle("Combined") effects.means.combined ggsave(file="pred.v.predator.Fig2 Forest Plot_combined.pdf", plot=effects.means.combined, width=7, height=10) ##CREATE DUMMY FIGURE FOR THE LEGEND effects.means.dummy<-ggplot(combined.data, aes(reorder(entry_id, yi), yi), group = trait_gen_b_b) effects.means.dummy<-effects.means.dummy + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.dummy<-effects.means.dummy + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.dummy<-effects.means.dummy + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.dummy<-effects.means.dummy + scale_colour_manual("trait", values = forest.colors) effects.means.dummy<-effects.means.dummy + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 16) effects.means.dummy<-effects.means.dummy + theme( text = element_text(family = "Times", size = 14), axis.text=element_text(size=12), axis.text.y=element_blank(), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.dummy<-effects.means.dummy + scale_y_continuous(limits = c(-7,7), breaks=seq(-7, 7, 1)) effects.means.dummy<-effects.means.dummy + ggtitle("Combined") effects.means.dummy ##create combined plot with all three forest plots library(grid) pdf(file = paste("Pred.V.parasite.Fig2.pdf"), width = 12, height = 14) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) grid.newpage() pushViewport(viewport(layout = grid.layout(2, 2))) print(effects.means.predator, vp = vplayout(1, 1)) print(effects.means.parasite, vp = vplayout(1, 2)) print(effects.means.combined, vp = vplayout(2, 1)) dev.off() ##SAVE DUMMY FIGURE TO ADD LATER IN POWERPOINT ggsave(file="pred.v.parasite.Fig2.legend.pdf", plot=effects.means.dummy, width=10, height=8) ``` *FIGURE 3A-B* Stregnth of trait responses to predators, parasites, and their combined presence, and the influence of consumer state on the strength of trait responses ```{r PLOT THE MEAN EFFECT SIZE OF EACH INTERACTION TYPE, include = F} library(ggplot2) effects.by.interaction<-ggplot(raw.data.sub2, aes(interaction_type, yi), group = interaction_type) effects.by.interaction<-effects.by.interaction + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = .2, size = .8, aes(group = interaction_type), color = "black") + stat_summary(fun.data = "mean_se", geom = "crossbar", aes(fill = interaction_type), width = 0.2, size = .2, position=position_dodge(width=0.1)) effects.by.interaction<-effects.by.interaction + scale_fill_manual(values=c("predator-prey"="#6699CC", "host-parasite"="#CC0000", combined = "grey")) effects.by.interaction<-effects.by.interaction + labs(x = "Interaction type", y = "Mean response magnitude\n Hedge's d ± 95% CI)") + theme_bw(base_size = 16) effects.by.interaction<-effects.by.interaction + theme(text = element_text(family = "Times", size = 22), axis.text=element_text(size=20), legend.position = c(.5, .8), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.x=element_blank()) effects.by.interaction<-effects.by.interaction + coord_cartesian(ylim = c(-0.5,2)) + geom_segment(aes(x = 0.4, xend=3.6,y=-0,yend=0),color = "grey65", linetype = "dashed") effects.by.interaction<-effects.by.interaction + scale_y_continuous(breaks = seq(-0.5, 2, .5)) + annotate("text", x = 0.6, y = 2.0, label = "(a)", size = 6) effects.by.interaction<-effects.by.interaction + scale_x_discrete(labels = c("predator-\nprey", "host-\nparasite", "combined")) effects.by.interaction ggsave(file="pred.v.parasite.Fig3.pdf", plot=effects.by.interaction, width=8, height=6) ``` ```{r PLOT EFFECTS BY CONSUMER STATE} raw.data.sub2$consumer_state<-factor(raw.data.sub2$consumer_state, levels = c("questing", "attacking", "consuming")) ##add dummy rows with predator-prey entries for attacking and consuming states. This will standardize the column widths raw.data.sub2.alt<-raw.data.sub2[,c(5,19, 54)] pred.attack.dummy<-c("predator-prey", "attacking", 0.00) pred.consume.dummy<-c("predator-prey", "consuming", 0.00) raw.data.sub2.alt<-rbind(raw.data.sub2.alt, pred.attack.dummy, pred.consume.dummy) raw.data.sub2.alt$yi<-as.numeric(raw.data.sub2.alt$yi) library(ggplot2) effects.by.state<-ggplot(raw.data.sub2.alt[which(raw.data.sub2.alt$interaction_type != "combined"),], aes(consumer_state, yi), group = interaction_type) effects.by.state<-effects.by.state + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = .4, size = .8, aes(group = interaction_type), color = "black", position = position_dodge(width = .5)) + stat_summary(fun.data = "mean_se", geom = "crossbar", aes(fill = interaction_type), width = 0.4, size = .2, position = position_dodge(width = .5)) effects.by.state<-effects.by.state + scale_fill_manual(values=c("predator-prey"="#6699CC", "host-parasite"="#CC0000", "combined" = "grey")) effects.by.state<-effects.by.state + labs(x = "Response type", y = "Mean response magnitude \n (Hedge's d ±- 95% CI)") + theme_bw(base_size = 16) effects.by.state<-effects.by.state + theme(text = element_text(family = "Times", size = 22), axis.text=element_text(size=20), panel.grid.minor = element_blank(), legend.title=element_blank(), legend.position= c(.5, .8), panel.grid.major = element_blank(), axis.ticks.x=element_blank()) effects.by.state<-effects.by.state + coord_cartesian(ylim = c(-1.2,2)) + geom_segment(aes(x = 0.4, xend=3.6,y=-0,yend=0),color = "grey65", linetype = "dashed") effects.by.state<-effects.by.state + scale_y_continuous(breaks = seq(-1, 2, .5)) + scale_x_discrete(labels = c("avoid", "counter", "combat")) + annotate("text", x = 0.6, y = 2.0, label = "(b)", size = 6) effects.by.state ggsave(file="pred.v.parasite.Fig3B.pdf", plot=effects.by.state, width=10, height=8) ``` Integrate the above two plots into a single, two-panel figure using the 'grid' package ```{r } library(grid) pdf(file = paste("Pred.parasite.Fig3 mean responses by factors.pdf"), width = 12, height = 6) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(effects.by.interaction, vp = vplayout(1, 1)) print(effects.by.state, vp = vplayout(1, 2)) dev.off() ``` **SUPPORTING INFORMATION** *Summary statistics* To get a general summary of effect directions, we first calculate counts of positive and negative effects. Results are reported in the supplementary material. ```{r} #PREDATORS #positive = 20 nrow(predator.data[which(predator.data$observed_effect == "positive"),]) #negative = 1 nrow(predator.data[which(predator.data$observed_effect == "negative"),]) #none = 21 nrow(predator.data[which(predator.data$observed_effect == "none"),]) #PARASITES #positive = 9 nrow(parasite.data[which(parasite.data$observed_effect == "positive"),]) #negative = 4 nrow(parasite.data[which(parasite.data$observed_effect == "negative"),]) #none = 24 nrow(parasite.data[which(parasite.data$observed_effect == "none"),]) #COMBINED #positive = 15 nrow(combined.data[which(combined.data$observed_effect == "positive"),]) #negative = 3 nrow(combined.data[which(combined.data$observed_effect == "negative"),]) #none = 21 nrow(combined.data[which(combined.data$observed_effect == "none"),]) ``` Here, will will calculate and plot the number of data entries across the different factors tested in our analyses above ```{r pressure, echo=FALSE, fig.width=10, fig.height=12} ##total number of entries per interaction type table(raw.data$interaction_type) ##number of studies used length(unique(raw.data$study_id)) ##number of host-parasite entries per consumer state table(parasite.data$consumer_state) ##CONSUMER STATE ##________________ ##create table that lists the count of entries by consumer state and interaction type entries.table<-aggregate(entry_id ~ interaction_type + consumer_state, data = raw.data.sub2, FUN = length) entries.table[7,]<-c("predator-prey", "attacking", 0) entries.table[8,]<-c("predator-prey", "consuming", 0) entries.table$entry_id<-as.numeric(entries.table$entry_id) entries.table$consumer_state<-factor(entries.table$consumer_state, levels = c("questing", "attacking", "consuming")) # ``` ```{r} ##create barplot for count of entries for each consumer state library(ggplot2) ``` ```{r} #create barplot of count of entries by consumer state counts.by.state<-ggplot(entries.table[which(entries.table$interaction_type != "combined"),], aes(x = consumer_state, y = entry_id), group = interaction_type) counts.by.state<-counts.by.state + geom_bar(stat = "identity", aes(fill = interaction_type), color = "black", width = .5, position = "dodge") counts.by.state<-counts.by.state + scale_fill_brewer(palette = "Greys") counts.by.state<-counts.by.state + labs(x = "consumer state", y = "number of entries") + theme_bw(base_size = 18) counts.by.state<-counts.by.state + theme(legend.position = c(.5, .8), text = element_text(family = "Times", size = 18), axis.text=element_text(angle = 10, vjust = 0.6, size = 18), axis.text.x=element_text(size = 14), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.x=element_blank()) counts.by.state<-counts.by.state + annotate("text", x = .55, y = 40, label = "a.)") + coord_cartesian(ylim = c(0,45)) ##TRAIT CLASS## ##_____________ ##create a table that counts the number of studies for each trait class entries.table.traits<-aggregate(entry_id ~ trait_gen_b + interaction_type, data = raw.data.sub2, FUN = length) entries.table.traits$entry_id<-as.numeric(entries.table.traits$entry_id) ##create barplot of count of entries by trait class counts.by.trait.gen<-ggplot(entries.table.traits, aes(x = trait_gen_b, y = entry_id), group = interaction_type) counts.by.trait.gen<-counts.by.trait.gen + geom_bar(stat = "identity", aes(fill = interaction_type), color = "black", width = .5, position = "dodge") counts.by.trait.gen<-counts.by.trait.gen + scale_fill_brewer(palette = "Greys") counts.by.trait.gen<-counts.by.trait.gen + labs(x = "trait type", y = "number of entries") + theme_bw(base_size = 18) counts.by.trait.gen<-counts.by.trait.gen + theme(legend.position = c(.5, .7), text = element_text(family = "Times", size = 18), axis.text=element_text(angle = 10, vjust = 0.6, size = 18), axis.text.x=element_text(size = 14), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.x=element_blank()) counts.by.trait.gen<-counts.by.trait.gen + annotate("text", x = .5, y = 42, label = "b.)") + coord_cartesian(ylim = c(0,45)) ##scale# ##_____________ ##create a table that counts the number of studies for each trait class entries.table.scale<-aggregate(entry_id ~ analysis_scale + interaction_type, data = raw.data.sub2, FUN = length) entries.table.scale$entry_id<-as.numeric(entries.table.scale$entry_id) ##create barplot of count of entries by resource group counts.by.scale<-ggplot(entries.table.scale, aes(x = analysis_scale, y = entry_id), group = interaction_type) counts.by.scale<-counts.by.scale + geom_bar(stat = "identity", aes(fill = interaction_type), color = "black", width = .5, position = "dodge") counts.by.scale<-counts.by.scale + scale_fill_brewer(palette = "Greys") counts.by.scale<-counts.by.scale + labs(x = "analysis scale", y = "number of entries") + theme_bw(base_size = 18) counts.by.scale<-counts.by.scale + theme(legend.position = c(.5, .8), text = element_text(family = "Times", size = 18), axis.text=element_text(size=20), axis.text.x=element_text(angle = 20, vjust = 0.6, size = 18), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.x=element_blank()) counts.by.scale<-counts.by.scale + annotate("text", x = .5, y = 42, label = "c.)") + coord_cartesian(ylim = c(0,45)) ##combine plates into one figure ##Combine Figures library(grid) pdf(file = paste("Pred.vs.Parasite.FigS2 summary statistics.pdf"), width = 14, height = 10) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) grid.newpage() pushViewport(viewport(layout = grid.layout(2, 2))) print(counts.by.state, vp = vplayout(1, 1)) print(counts.by.trait.gen, vp = vplayout(1, 2)) print(counts.by.scale, vp = vplayout(2, 1)) dev.off() ``` Here, we will perform multi-model inference using the glmulti package for validation of our main analyses ```{r } library(lme4) library(MASS) ##create co-variance matrices ##calculate covariation matrix to account for same control groups for multiple treatments calc.v <- function(x) { v <- matrix(1/x$n_treat[1] + outer(x$yi, x$yi, "*")/(2*x$n_tot[1]), nrow=nrow(x), ncol=nrow(x)) diag(v) <- x$vi_2 v } ##calculate covariation matrix to account for same control groups for multiple treatments V<-bldiag(lapply(split(raw.data, list(raw.data$study_id, raw.data$grouping_factor), drop = T), calc.v)) V.sub2<- bldiag(lapply(split(raw.data.sub2, list(raw.data.sub2$study_id, raw.data.sub2$grouping_factor), drop = T), calc.v)) ##__________________________ ## MODEL SELECTION ##_________________________ ##Load glmulti and define function for fitting random-effects model library(glmulti) ##NOTE: it took me forever to get this package to load with the latest R version. The problem was that a dependency, rJava, was not properly directed to java files in my harddrive. I had to update where R looks for the Java files by updated Jave Runtime Environment to the latest 64 bit version, and then, CRUCIALLY, entering the following in command line (terminal): sudo R CMD javareconf rma.glmulti <- function(formula, data, ...) { rma.mv(formula, V.sub2, data=data, random = ~1|study_id/entry_id, method="ML", ...) } full.model<- glmulti(yi ~ interaction_type*analysis_scale + trait_gen_b, data= raw.data.sub2, level=1, fitfunction=rma.glmulti, crit="aicc", confsetsize=128) ##level denotes whether interaction terms are included (2 for yes, 1 for no) print(full.model) ###summary results plot(full.model) ###plot models by aicc value ##list the models in order of AICc ranking tmp <- weightable(full.model) tmp ##plot variable importances plot(full.model, type="s") ##details of full model summary(full.model@objects[[1]]) ##______________________________ ##MODEL AVERAGING ##_____________________________ setOldClass("rma.mv") setMethod('getfit', 'rma.mv', function(object, ...) { if (object$test=="z") { cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=100000) } else { cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=object$k-object$p) } }) round(coef(full.model), 4) preds <- list() for (j in 1:full.model@nbmods) { model <- full.model@objects[[j]] vars <- names(coef(model))[-1] if (length(vars) == 0) { preds[[j]] <- predict(model) } else { preds[[j]] <- predict(model, newmods=x[vars]) } } weights <- weightable(full.model)$weights yhat <- sum(weights * sapply(preds, function(x) x$pred)) round(yhat, 3) ``` Now, we will test for publication bias using 3 standard methods: Rank correlation testing, Egger's regression testing, and Rosenthal fail-safe numbering. We will perform the tests with the full dataset, as well as the two model subsets - parasite data subset and the questing state data subset - used in our main analyses. ```{r } ## FULL DATASET ##_________________ ##basic model base.mm.1<-rma(yi, vi_2, mods = ~interaction_type*analysis_scale + trait_gen_b, data = raw.data) print(base.mm.1) ##rank correlation test ranktest(base.mm.1) ##egger's regression test regtest(base.mm.1) ##compute Rosenthal's fail-safe number fsn(yi, vi_2, type = "Rosenthal", alpha = .05, data = raw.data, digits = 3) ##QUESTING MODEL ##_______________ ##create subset of questing state questing.data<-droplevels(subset(raw.data.sub, consumer_state == "questing")) ##base model base.mm.2<-rma(yi, vi_2, mods = ~interaction_type, data = questing.data) print(base.mm.2) ##rank correlation test ranktest(base.mm.2) ##egger's regression test regtest(base.mm.2) ##compute Rosenthal's fail-safe number fsn(yi, vi_2, type = "Rosenthal", alpha = .05, data = questing.data, digits = 3) ##PARASITE MODEL ##_______________ ##base model base.mm.3<-rma(yi, vi_2, data = parasite.data) print(base.mm.3) ##rank correlation test ranktest(base.mm.3) ##egger's regression test regtest(base.mm.3) ##compute Rosenthal's fail-safe number fsn(yi, vi_2, type = "Rosenthal", alpha = .05, data = parasite.data, digits = 3) ##_________________________________ ##Generate Figure of funnel plots ##________________________________ par(mar=c(5,4,1,2)) par(mfrow=c(2,2)) funnel(base.mm.1, main = "Overall model", xlab = "Effect size") funnel.2<-funnel(base.mm.2, main = "Questing model", xlab = "Effect size") funnel(base.mm.3, main = "Parasite model", xlab = "Effect size") ``` **Figure S6** Alternative Forest plot of trait response magnitudes forest plot sorted by study (shown in supplementary material) ```{r} library(ggplot2) ##relevel trait type factor levels raw.data.sub2$trait_gen_b<-factor(raw.data.sub2$trait_gen_b, levels = c("activity", "space use", "morphology/physiology")) ##set color pallete forest.colors<-c("#6699CC", "#990000", "#333333") ##create subsets parasite.data<-droplevels(subset(raw.data.sub2, interaction_type == "host-parasite")) predator.data<-droplevels(subset(raw.data.sub2, interaction_type == "predator-prey")) combined.data<-droplevels(subset(raw.data.sub2, interaction_type == "combined")) ##PARASITES effects.means.parasite<-ggplot(parasite.data, aes(x = reorder(entry_id, study_no), y = yi), group = trait_gen_b) effects.means.parasite<-effects.means.parasite + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.parasite<-effects.means.parasite + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.parasite<-effects.means.parasite + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.parasite<-effects.means.parasite + scale_colour_manual("trait class", values = forest.colors) effects.means.parasite<-effects.means.parasite + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 18) effects.means.parasite<-effects.means.parasite + theme(legend.position = "none",text = element_text(family = "Times", size = 16), axis.text=element_text(size=14), axis.text.y=element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.grid.major = element_blank(), legend.title=element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.parasite<-effects.means.parasite + scale_y_continuous(limits = c(-10,7), breaks=seq(-7, 7, 1), labels = seq(-7, 7, 1)) effects.means.parasite<-effects.means.parasite + geom_text(aes(y=-7,label=study_id),size= 4) effects.means.parasite<-effects.means.parasite + ggtitle("Parasites") effects.means.parasite ggsave(file="pred.v.parasite.FigS6a_SI Forest Plot_parasites.pdf", plot=effects.means.parasite, width=7, height=10) ##PREDATORS effects.means.predator<-ggplot(predator.data, aes(reorder(entry_id, study_no), yi), group = trait_gen_b_b) effects.means.predator<-effects.means.predator + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.predator<-effects.means.predator + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.predator<-effects.means.predator + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.predator<-effects.means.predator + scale_colour_manual("trait", values = forest.colors) effects.means.predator<-effects.means.predator + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 18) effects.means.predator<-effects.means.predator + theme( legend.position = "none",text = element_text(family = "Times", size = 16), axis.text=element_text(size=14), axis.text.y=element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.grid.major = element_blank(), legend.title=element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.predator<-effects.means.predator + scale_y_continuous(limits = c(-7,7), breaks=seq(-7, 7, 1)) effects.means.predator<-effects.means.predator + geom_text(aes(y=-5,label=study_id),size= 4) effects.means.predator<-effects.means.predator + ggtitle("Predators") effects.means.predator ggsave(file="pred.v.parasite.FigS6a_SI Forest Plot_predators.pdf", plot=effects.means.predator, width=7, height=10) ##COMBINED effects.means.combined<-ggplot(combined.data, aes(reorder(entry_id, study_no), yi), group = trait_gen_b_b) effects.means.combined<-effects.means.combined + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.combined<-effects.means.combined + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.combined<-effects.means.combined + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.combined<-effects.means.combined + scale_colour_manual("trait", values = forest.colors) effects.means.combined<-effects.means.combined + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 18) effects.means.combined<-effects.means.combined + theme( legend.position = "none", text = element_text(family = "Times", size = 16), axis.text=element_text(size=14), axis.text.y=element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.grid.major = element_blank(), legend.title=element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.combined<-effects.means.combined + scale_y_continuous(limits = c(-7,7), breaks=seq(-7, 7, 1)) effects.means.combined<-effects.means.combined + geom_text(aes(y=-5,label=study_id),size= 4) effects.means.combined<-effects.means.combined + ggtitle("Combined") effects.means.combined ggsave(file="pred.v.parasite.FigS6a_SI Forest Plot_combined.pdf", plot=effects.means.combined, width=7, height=10) ##CREATE DUMMY FIGURE FOR THE LEGEND effects.means.dummy<-ggplot(combined.data, aes(reorder(entry_id, study_no), yi), group = trait_gen_b_b) effects.means.dummy<-effects.means.dummy + geom_point(size = 2, aes(color = trait_gen_b)) effects.means.dummy<-effects.means.dummy + geom_hline(yintercept = 0, color = "black", size = 0.5) effects.means.dummy<-effects.means.dummy + geom_errorbar(aes(ymin= yi - 1.96*sqrt(vi_2), ymax= yi + 1.96*sqrt(vi), color = trait_gen_b)) + coord_flip() effects.means.dummy<-effects.means.dummy + scale_colour_manual("trait", values = forest.colors) effects.means.dummy<-effects.means.dummy + labs(x = "", y = "Response magnitude (Hedge's d ± 95% CI)") + theme_bw(base_size = 16) effects.means.dummy<-effects.means.dummy + theme( text = element_text(family = "Times", size = 14), axis.text=element_text(size=12), axis.text.y=element_blank(), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), plot.title = element_text(hjust = 0.5)) effects.means.dummy<-effects.means.dummy + scale_y_continuous(limits = c(-7,7), breaks=seq(-7, 7, 1)) effects.means.dummy<-effects.means.dummy + ggtitle("Combined") effects.means.dummy ``` Integate the 3 forest plots using the grid package ```{r} library(grid) pdf(file = paste("Pred.V.parasite.FigS6a.pdf"), width = 10, height = 12) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) grid.newpage() pushViewport(viewport(layout = grid.layout(2, 2))) print(effects.means.predator, vp = vplayout(1, 1)) print(effects.means.parasite, vp = vplayout(1, 2)) print(effects.means.combined, vp = vplayout(2, 1)) dev.off() ``` **Figure S7.** Plot effect sizes of trait responses for (a) different trait types and (b) analysis scales ```{r} library(ggplot2 ) ##(a) Trait type ##relevel trait classes raw.data.sub2$trait_gen_b<-factor(raw.data.sub2$trait_gen_b, levels = c("activity", "space use", "morphology/physiology")) library(ggplot2) effects.by.trait<-ggplot(raw.data.sub2, aes(trait_gen_b, yi), group = interaction_type) effects.by.trait<-effects.by.trait + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "black", width = .4, size = .9, position = position_dodge(width = .5), aes(group = interaction_type)) + stat_summary(fun.data = "mean_se", geom = "crossbar", aes(fill = interaction_type), width = 0.4, size = .4, position = position_dodge(width = .5)) effects.by.trait<-effects.by.trait + scale_fill_manual(values=c("predator-prey"="#6699CC", "host-parasite"="#CC0000", "combined" = "grey")) effects.by.trait<-effects.by.trait + labs(x = "trait type", y = "Mean response magnitude \n (Hedge's d ± 95% CI)") + theme_bw(base_size = 16) effects.by.trait<-effects.by.trait + theme(text = element_text(family = "Times", size = 22), axis.text=element_text(size=20), legend.position = c(.5, .8), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.x=element_blank()) effects.by.trait<-effects.by.trait + coord_cartesian(ylim = c(-1,2)) + geom_segment(aes(x = 0.4, xend=3.6,y=-0,yend=0),color = "grey65", linetype = "dashed") effects.by.trait<-effects.by.trait + scale_y_continuous(breaks = seq(-1, 2, .5)) effects.by.trait<-effects.by.trait + scale_x_discrete(labels = c("activity", "space\nuse", "morphology/\nphysiology")) + annotate("text", x = 0.55, y = 2.0, label = "(a)", size = 6) effects.by.trait ggsave(file="pred.v.parasite.FigS7b.pdf", plot=effects.by.trait, width=8, height=6) ##(b) analysis scale effects.by.scale<-ggplot(raw.data.sub2, aes(analysis_scale, yi), group = interaction_type) effects.by.scale<-effects.by.scale + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "black", width = .3, size = .8, position = position_dodge(width = .5), aes(group = interaction_type)) + stat_summary(fun.data = "mean_se", geom = "crossbar", aes(fill = interaction_type), width = 0.3, size = .3, position = position_dodge(width = .5)) effects.by.scale<-effects.by.scale + scale_fill_manual(values=c("predator-prey"="#6699CC", "host-parasite"="#CC0000", "combined" = "grey")) effects.by.scale<-effects.by.scale + labs(x = "analysis scale", y = "Mean response magnitude\n (Hedge's d ± 95% CI)") + theme_bw(base_size = 16) effects.by.scale<-effects.by.scale + theme(text = element_text(family = "Times", size = 22), axis.text=element_text(size=20), legend.position = c(.7, .8), panel.grid.minor = element_blank(), legend.title=element_blank(), panel.grid.major = element_blank(), axis.ticks.x=element_blank()) effects.by.scale<-effects.by.scale + coord_cartesian(ylim = c(-1,2)) + geom_segment(aes(x = 0.4, xend=2.6,y=-0,yend=0),color = "grey65", linetype = "dashed") effects.by.scale<-effects.by.scale + scale_y_continuous(breaks = seq(-1, 2, .5)) + annotate("text", x = 0.50, y = 2.0, label = "(b)", size = 6) effects.by.scale ggsave(file="pred.v.parasite.FigS7b.pdf", plot=effects.by.scale, width=8, height=6) ``` Integrate the above two plots into a single, two-panel figure using the 'grid' package ```{r } library(grid) pdf(file = paste("Pred.parasite.FigS7 mean responses by more factors.pdf"), width = 12, height = 6) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(effects.by.trait, vp = vplayout(1, 1)) print(effects.by.scale, vp = vplayout(1, 2)) dev.off() ```