# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # POLLINATOR VISITATION AND HABITAT TYPE # # Section "Pollinators and habitat type" # Figs. 6 and 7 # # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ library(dplyr) library(lme4) library(blmeco) library(ggplot2) library(ggeffects) library(emmeans) mythem <- theme_bw(22) + theme(legend.position="none", legend.key = element_blank(), strip.background = element_blank(), panel.grid = element_blank(), panel.border= element_rect(size=rel(1.75),color='black'), strip.text.x=element_text(vjust = 1)) master.data <- read.csv("master.data.csv", stringsAsFactors = FALSE) %>% filter(!insect.order %in% c("Neuroptera", "Hemiptera", "Orthoptera")) %>% mutate(ano = factor(ano), fecha = as.Date(fecha, format="%Y-%m-%d"), flores.scaled = scale(flores), day.scaled = scale(day.of.year), day.scaled2 = scale(day.of.year)^2, habitat.type=factor(habitat.type)) # 1. ALL POLLINATORS COMBINED ------------ orden <- "All" tmp1 <- master.data %>% filter(insect.order==orden) # PATCH VISITATION fit1 <- glmer(event ~ habitat.type + day.scaled + day.scaled2 + flores.scaled + (1|especie) + (1|sitio) + (1|ano), data = tmp1, family ='binomial') print(dispersion_glmer(fit1)) print(car::Anova(fit1, test="Chisq")) predicted.event <- ggeffects::ggpredict(fit1, terms = "habitat.type", type="fe", ci.lvl = 0.95) predicted.event$response <- factor("Patch visitation") predicted.event$group <- factor(orden) assign(paste0("predicted.event.",orden), predicted.event) assign(paste0("fit1.", orden), fit1) # FLOWER VISITATION tmp1$olre <- 1:nrow(tmp1) fit2 <- glmer(visrate ~ habitat.type + day.scaled + day.scaled2 + flores.scaled + (1|especie) + (1|sitio) + (1|ano) + (1|olre), data = tmp1, family ='binomial', weights = flores) print(dispersion_glmer(fit2)) print(car::Anova(fit2, test="Chisq")) predicted.visrate <- ggeffects::ggpredict(fit2, terms = "habitat.type", type="fe", ci.lvl = 0.95) predicted.visrate$response <- factor("Flower visitation") predicted.visrate$group <- factor(orden) assign(paste0("predicted.visrate.",orden), predicted.visrate) assign(paste0("fit2.", orden), fit2) # Marginal effects ------ FIG. 6 # combine into single dataframe predicted.all.combined <- rbind(predicted.event.All, predicted.visrate.All) habitats <- sort(unique(master.data$habitat.type)) # FIG. 6 ------------------- fig <- ggplot(predicted.all.combined, aes(factor(x), predicted)) + geom_linerange(aes(ymin = conf.low, ymax = conf.high), color="darkgrey", size=1) + geom_point(size=3.5, color="blue") + facet_wrap (~response, scales = "free_y") + scale_y_log10() fig <- fig + scale_x_discrete(labels=habitats) + mythem + theme(axis.text.x = element_text(size=rel(1),angle=45,hjust=1)) + labs(x = "Habitat type", y = "Visitation probability") print(fig) # Uncomment to save Fig 6 as pdf file # ggsave(file = "Fig.6.pdf", height=7, width = 10) # 2. BY ORDER ----------------- # Interaction Habitat x Insect Order # PATCH VISITATION ordenes <- c("Coleoptera", "Diptera", "Hymenoptera", "Lepidoptera") tmp2 <- master.data %>% filter(insect.order %in% ordenes) fit3 <- glmer(event ~ habitat.type + insect.order + habitat.type: insect.order + day.scaled + day.scaled2 + flores.scaled + (1|especie) + (1|sitio) + (1|ano), data = tmp2, family ='binomial') fit4 <- glmer(event ~ habitat.type + insect.order + day.scaled + day.scaled2 + flores.scaled + (1|especie) + (1|sitio) + (1|ano), data = tmp2, family ='binomial') anova(fit3, fit4) print(dispersion_glmer(fit3)) print(car::Anova(fit3, test="Chisq")) predict.patch <- ggpredict(fit3, ~ habitat.type | insect.order, type="fe") predict.patch$response <- "Patch visitation" # Heterogeneity among habitats separately by order emm3 <- emmeans(fit3, ~habitat.type|insect.order, type="response") test(pairs(emm3), joint = TRUE, adjust="Tukey") # FLOWER VISITATION fit5 <- glmer(visrate ~ habitat.type + insect.order + habitat.type: insect.order + day.scaled + day.scaled2 + flores.scaled + (1|especie) + (1|sitio) + (1|ano), data = tmp2, family ='binomial', weights=flores) fit6 <- glmer(visrate ~ habitat.type + insect.order + day.scaled + day.scaled2 + flores.scaled + (1|especie) + (1|sitio) + (1|ano), data = tmp2, family ='binomial', weights=flores) anova(fit5, fit6) print(dispersion_glmer(fit5)) print(car::Anova(fit5, test="Chisq")) predict.visrate <- ggpredict(fit5, ~ habitat.type | insect.order, type="fe") predict.visrate$response <- "Flower visitation" # Heterogeneity between habitats, separately by insect order emm5 <- emmeans(fit5, ~habitat.type|insect.order, type="response") test(pairs(emm5), joint = TRUE, adjust="Tukey") # FIG. 7 ------------------- predict.all <- rbind(predict.patch, predict.visrate) predict.all$response <- factor(predict.all$response, levels=c("Patch visitation", "Flower visitation"), ordered =TRUE ) fig1 <- ggplot(predict.all, aes(x=factor(x), y=predicted, color= group)) + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(width=0.35)) + facet_wrap (~response, scales = "free_y") + scale_y_log10() + ggthemes::scale_color_gdocs() fig1 <- fig1 + scale_x_discrete(labels=habitats) + mythem + theme(axis.text.x = element_text(size=rel(1),angle=45,hjust=1)) + labs(x = "Habitat type", y = "Visitation probability") fig1 <-fig1 + theme(legend.position="top", legend.margin=margin(0,0,0,0),legend.title = element_blank(),legend.key.height=unit(2.5,"line"), legend.box.margin=margin(c(-10,-10,-10,-10))) print(fig1) # uncomment to save the figure as a pdf file # ggsave(file = "Fig.7.pdf", height=7, width = 10) # clean up # rm(list=ls(pattern="^tmp")) # rm(list=ls(pattern="^fit")) # rm(list=ls(pattern="^predic"))