# ------------------------------------------------------- # Script 1 # TESTS #1 AND #2 IN TABLE 2 # Section "Individual plant species: patch visitation" # ------------------------------------------------------- # This script runs with these versions of # required packages. May not work with later versions. require(car); # ver. 3.0-0 require(broom); # ver. 0.4.5 require(purrr); # ver. 0.2.5 require(ggplot2); # ver. 3.0.0 require(ggthemes); # ver. 4.0.0 require(dplyr) ; # ver. 0.7.6 my_se <- function(x) sqrt(var(x,na.rm=TRUE)/length(which(!is.na(x)))) working.dataset <- read.csv("working.dataset.csv", stringsAsFactors = FALSE) #------------------------------- # PATCH VISITATION PROBABILITY, # BY INDIVIDUAL PLANT SPECIES, # ALL INSECT ORDERS COMBINED # # TEST #1 IN TABLE 2 #------------------------------- event.all <- working.dataset %>% group_by(especie, sitio, ano.as.factor, ano, mes, fecha, censo.num) %>% summarize(flores=mean(flores), visitas =sum(visitas)) %>% mutate(event = ifelse(visitas> 0, 1, 0)) resumen.modelos.patch <- event.all %>% split(.$especie) %>% map_df(~ tidy(car::Anova(glm(event ~ ano.as.factor + scale(flores), family="binomial", data=., maxit = 500), type=2, test.statistic="LR")[1,3])) slopes <- event.all %>% split(.$especie) %>% map_df(~ tidy(coefficients (glm(event ~ ano + flores, family="binomial", data=., maxit = 500))[2])) %>% dplyr::select(-names) names(slopes) <- "slope" resumen.modelos.patch$especie <- sort(unique(event.all$especie)) names(resumen.modelos.patch) <- c("pvalue", "especie") resumen.modelos.patch$slope <- slopes$slope ; rm(slopes) # Benjamini Hochberg correction p values resumen.modelos.patch$pvalueBH <- p.adjust(resumen.modelos.patch$pvalue, method='BH') label.p <- data.frame(especie=sort(unique(event.all$especie)), pvalue=resumen.modelos.patch$pvalueBH, ano = 1999.5, patchvisit=0.82) label.p$signif.code <- ifelse(label.p$pvalue < 0.05, "Sign", "Nonsign") # freqs slope signs in species with significant heterogeneity resumen.modelos.patch$signif <- label.p$signif.code resumen.modelos.patch$sign <- ifelse(resumen.modelos.patch$slope>=0, 'Positive', 'Negative') with(resumen.modelos.patch, table(sign, signif, dnn=c("Sign of slope", "Stat. significance"))) # Significance of heterogeneity vs number sampling dates numdates <- working.dataset %>% group_by(especie) %>% summarize (dates = length(unique(fecha))) numdates <- full_join(label.p, numdates, by="especie") %>% group_by(signif.code) numdates %>% summarize(meandates=mean(dates), sedates=my_se(dates)) with(numdates, kruskal.test(dates~factor(signif.code))) ; rm(numdates) # FIG. 3 --------- patchvisit <- event.all %>% group_by(especie,ano.as.factor, ano) %>% summarise(patchvisit=mean(event), sepatchvisit=my_se(event)) patchvisit$order <- 'All' fi <- ggplot(patchvisit, aes(x=ano, y=patchvisit, group=especie)) + geom_point() + geom_line() + facet_wrap(~especie, nrow = 13) fi <- fi + geom_rect(data = label.p, aes(fill =signif.code),xmin = -Inf,xmax = Inf, ymin = -Inf,ymax = Inf,alpha = 0.3, inherit.aes = FALSE) + scale_fill_manual(values=c(NA, "grey40")) fi <- fi + geom_errorbar(aes(ymin = patchvisit - 2*sepatchvisit, ymax = patchvisit + 2*sepatchvisit), width = 0.5) + geom_point(size=1.5)+ theme_bw() fi <- fi + theme(strip.text = element_text(face = "italic"), strip.background = element_blank()) + labs(x='Year', y ='Patch visitation probability') + theme(panel.spacing.y = unit(0.10, "lines")) + scale_y_continuous(breaks=c(0, 0.5, 1.0)) + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.spacing.y=unit(-0.2, "lines"), panel.spacing.x=unit(0.2, "lines")) fi <- fi + geom_text(data=label.p, aes(x=ano, y=patchvisit, label = paste0(format.pval(pvalue, digits=1))), inherit.aes = FALSE, size=3) print(fi) # uncomment to save figure in pdf file # ggsave(file='Fig.3.Patch.Visitation.by.species.pdf', width = 8.5, height=10.5, unit ='in') rm(fi, label.p) #------------------------------- # PATCH VISITATION PROBABILITY, # BY INDIVIDUAL PLANT SPECIES, # BY INSECT ORDERS # # TEST #2 IN TABLE 2 #------------------------------- tmp2 <- working.dataset %>% group_by (especie, sitio, ano.as.factor, ano, mes, fecha, censo.num, flores, order) %>% summarize(visitas = sum(visitas)) %>% mutate(event = ifelse(visitas>0, 1, 0)) # COLEOPTERA event.col <- tmp2 %>% mutate(event=ifelse(!order=='Coleoptera' | is.na(order), 0, 1)) %>% summarize(event=sum(event)) event.col$order='Coleoptera' resumen.col <- event.col %>% split(.$especie) %>% map_df(~ tidy(car::Anova(glm(event ~ ano.as.factor + flores, family="binomial", data=., maxit = 500))[1,3])) %>% mutate(x = ifelse(x==1, NA, x)) names(resumen.col) <- "p_col" resumen.col$especie <- sort(unique(event.col$especie)) resumen.col$order <- 'Coleoptera' resumen.col$p_colBH <- p.adjust(resumen.col$p_col, method="BH") patchvisit.col <- event.col %>% group_by(especie,ano.as.factor,ano) %>% summarise(patchvisit=mean(event), sepatchvisit=my_se(event)) patchvisit.col$order <- 'Coleoptera' slopes.col <- event.col %>% split(.$especie) %>% map_df(~ tidy(coefficients (glm(event ~ ano + flores, family="binomial", data=., maxit = 500))[2])) %>% dplyr::select(-names) table(slopes.col$x[resumen.col$p_colBH<0.05]>0, dnn="Coleoptera - Slope > 0") # DIPTERA event.dip <- tmp2 %>% mutate(event = ifelse(!order=='Diptera' | is.na(order), 0, 1)) %>% summarize (event=sum(event)) event.dip$order='Diptera' resumen.dip <- event.dip %>% split(.$especie) %>% map_df(~ tidy(car::Anova(glm(event ~ ano.as.factor + flores, family="binomial", data=., maxit = 500))[1,3])) %>% mutate(x = ifelse(x==1, NA, x)) names(resumen.dip) <- "p_dip" resumen.dip$especie <- sort(unique(event.dip$especie)) resumen.dip$order <- 'Diptera' resumen.dip$p_dipBH <- p.adjust(resumen.dip$p_dip, method="BH") patchvisit.dip <- event.dip %>% group_by(especie,ano.as.factor,ano) %>% summarise(patchvisit=mean(event), sepatchvisit=my_se(event)) patchvisit.dip$order <- 'Diptera' slopes.dip <- event.dip %>% split(.$especie) %>% map_df(~ tidy(coefficients (glm(event ~ ano + flores, family="binomial", data=., maxit = 500))[2])) %>% dplyr::select(-names) table(slopes.dip$x[resumen.dip$p_dipBH<0.05]>0, dnn="Diptera - Slope > 0") # HYMENOPTERA event.hym <- tmp2 %>% mutate(event = ifelse(!order=='Hymenoptera' | is.na(order), 0, 1)) %>% summarize (event=sum(event)) event.hym$order='Hymenoptera' resumen.hym <- event.hym %>% split(.$especie) %>% map_df(~ tidy(car::Anova(glm(event ~ ano.as.factor + flores, family="binomial", data=., maxit = 500))[1,3])) %>% mutate(x = ifelse(x==1, NA, x)) names(resumen.hym) <- "p_hym" resumen.hym$especie <- sort(unique(event.hym$especie)) resumen.hym$order <- 'Hymenoptera' resumen.hym$p_hymBH <- p.adjust(resumen.hym$p_hym, method="BH") patchvisit.hym <- event.hym %>% group_by(especie,ano.as.factor,ano) %>% summarise(patchvisit=mean(event), sepatchvisit=my_se(event)) patchvisit.hym$order <- 'Hymenoptera' slopes.hym <- event.hym %>% split(.$especie) %>% map_df(~ tidy(coefficients (glm(event ~ ano + flores, family="binomial", data=., maxit = 500))[2])) %>% dplyr::select(-names) table(slopes.hym$x[resumen.hym$p_hymBH<0.05]>0, dnn="Hymenoptera - Slope > 0") # LEPIDOPTERA event.lep <- tmp2 %>% mutate(event = ifelse(!order=='Lepidoptera' | is.na(order), 0, 1)) %>% summarize (event=sum(event)) event.lep$order='Lepidoptera' resumen.lep <- event.lep %>% split(.$especie) %>% map_df(~ tidy(car::Anova(glm(event ~ ano.as.factor + flores, family="binomial", data=., maxit = 500))[1,3]))%>% mutate(x = ifelse(x==1, NA, x)) names(resumen.lep) <- "p_lep" resumen.lep$especie <- sort(unique(event.lep$especie)) resumen.lep$order <- 'Lepidoptera' resumen.lep$p_lepBH <- p.adjust(resumen.lep$p_lep, method="BH") patchvisit.lep <- event.lep %>% group_by(especie,ano.as.factor,ano) %>% summarise(patchvisit=mean(event), sepatchvisit=my_se(event)) patchvisit.lep$order <- 'Lepidoptera' slopes.lep <- event.lep %>% split(.$especie) %>% map_df(~ tidy(coefficients (glm(event ~ ano + flores, family="binomial", data=., maxit = 500))[2])) %>% dplyr::select(-names) table(slopes.lep$x[resumen.lep$p_lepBH<0.05]>0, dnn="Lepidoptera - Slope > 0") # FIG. S1, APPENDIX S5 --------- patchvisit.all <- rbind(patchvisit.col, patchvisit.dip, patchvisit.hym, patchvisit.lep) names(patchvisit.all)[6] <- "Order" fi2 <- ggplot(patchvisit.all, aes(x=ano, y=patchvisit, group=Order, color=Order)) + geom_line() + theme_bw() + facet_wrap(~especie, nrow = 13) fi2 <- fi2 + theme(strip.text = element_text(face = "italic"), strip.background = element_blank()) + labs(x='Year', y ='Patch visitation probability') + theme(panel.spacing.y = unit(0.10, "lines")) + scale_y_continuous(breaks=c(0, 0.5, 1.0)) + theme(legend.position="top") fi2 <- fi2 + scale_color_gdocs() print(fi2) # uncomment to save figure in pdf file # ggsave(file='App.S5.Fig.S1.Patch.Visitation.by.species.by.order.pdf', width = 9.5, height=12.5, unit ='in') rm(fi2) # TABLE S1, APPENDIX S5 --------- resumen.pvalues.all.orders <- cbind(dplyr::select(resumen.col, especie,p_colBH), resumen.dip[,"p_dipBH"], resumen.hym[,"p_hymBH"], resumen.lep[, "p_lepBH"]) resumen.pvalues.all.orders[nrow(resumen.pvalues.all.orders)+1,2:5] <- apply(as.matrix(resumen.pvalues.all.orders[,2:5] < 0.05), 2, function(x) sum(x, na.rm=TRUE)) resumen.pvalues.all.orders[1:65, ncol(resumen.pvalues.all.orders)+1] <- apply(as.matrix(resumen.pvalues.all.orders[1:65,2:5] < 0.05), 1, function(x) sum(x,na.rm=TRUE)) names(resumen.pvalues.all.orders)[6] <- "signif.orders" write.table(format(resumen.pvalues.all.orders, digits=3), file="Table.S1.Appendix.S5.txt",quote = FALSE, sep='\t',row.names = FALSE) ; # External file required in Script 4 save(patchvisit.col,patchvisit.dip,patchvisit.hym,patchvisit.lep,file="patch.visitation.by.order.RObj") # File required in Script 4 # uncomment to clean up working environment # rm(list=(ls(pattern="^resumen"))) # rm(list=(ls(pattern="^patchvisit"))) # rm(list=(ls(pattern="^event"))) # rm(list=(ls(pattern="^tmp"))) # rm(list=(ls(pattern="^slopes")))