###Model of floral visitor abundance in three blueberry types #Libraries library(bayesplot) library(brms) library(rstan) library(plyr) library(dplyr) library(emmeans) library(gridExtra) library(modelr) library(scales) library(tidyr) library(tidybayes) library(ggplot2) blue_abun=read.csv("floral visitor abundance JAPPL-2020-00197.csv") #Priors blue.abun.n.prior=prior(normal(0,1.5),class="b")+ prior(normal(0,1),class="Intercept")+ prior(normal(0,1),class="sd",group="Year")+ prior(normal(0,1),class="sd",group="Year:Site") #fit model #model terms #Abundance = count of visitors #VAR = Blueberry type #Species = floral visitor species (other = t. carbonaria in rabbiteye, southern highbush, b. terrestris in northern highbush) #Time = survey time of day #duration = length of survey #Year = survey year #Site = survey block abun.brm=brm(Abundance~VAR*Species+Time+offset(log(duration))+(1|Year/Site), family=zero_inflated_negbinomial(), prior=blue.abun.n.prior, cores=4, iter=3000, control=list(adapt_delta=0.999, max_treedepth=20), seed=123, data=blue_abun) pp_check(abun.brm,nsamples = 100) bayes_R2(abun.brm) #compute HDIs and group differences abun_plot <- plot(conditional_effects(abun.brm, conditions=data.frame(duration=1), re_formula=NA)) abun_plot <- abun_plot[[4]]$data abun_plot$cld_group <- paste(abun_plot$Species,abun_plot$VAR, sep = "_") #Compact letter differences between floral visitor species and blueberry types ####code for calculating CLDs get_p_matrix <- function(df, only_low = TRUE) { # pre-define matrix out <- matrix(-1, nrow = ncol(df), ncol = ncol(df), dimnames = list(colnames(df), colnames(df))) for (i in seq_len(ncol(df))) { for (j in seq_len(ncol(df))) { out[i, j] <- mean(df[,i] < df[,j]) } } if (only_low) out[out > .5] <- 1- out[out > .5] out } cld_pmatrix <- function(model, pars, level = 0.05) { p_matrix <- get_p_matrix(model) lp_matrix <- (p_matrix < (level/2) | p_matrix > (1-(level/2))) cld <- multcompView::multcompLetters(lp_matrix)$Letters cld } #### abun_hdi <- abun.brm %>% emmeans( ~ Species|VAR,transform="response") %>% gather_emmeans_draws() %>% mutate(cld_group = paste(Species,VAR, sep = "_")) abun_cld=abun_hdi %>% ungroup() %>% ## to get rid of unneeded columns select(.value, cld_group, .draw) %>% spread(cld_group, .value) %>% select(-.draw) %>% ## we need to get rid of all columns not containing draws cld_pmatrix() abun_cld<- as.data.frame(abun_cld) abun_cld$cld_group <- rownames(abun_cld) rownames(abun_cld) <- NULL colnames(abun_cld)[1] <- "cld" ###Dataframe for plotting results and making plot abun_plot_data<-merge(abun_plot,abun_cld,by=c("cld_group")) abun_plot_data$VAR=as.character(abun_plot_data$VAR) colnames(abun_plot_data)=c("cld_group","Blueberry species","Pollinator species","duration", "Abundance","Time","Year","Site","cond__", "effect1__","effect2__","estimate","se","lower","upper","cld") abun_plot_data$`Blueberry species`=revalue(abun_plot_data$`Blueberry species`,c("SH" = "Southern highbush", "RE" = "Rabbiteye", "NH" = "Northern highbush")) abundy_melt5=blue_abun colnames(abundy_melt5)[3]=c("Blueberry species") colnames(abundy_melt5)[8]=c("Pollinator species") abundy_melt5$`Blueberry species`=revalue(abundy_melt5$`Blueberry species`,c("SH" = "Southern highbush", "NH" = "Northern highbush", "RE" = "Rabbiteye")) abundy_plot_df=rbind.fill(abun_plot_data,abundy_melt5) abundy_plot_df$`Pollinator species` <- as.character(abundy_plot_df$`Pollinator species`) abundy_plot_df[(abundy_plot_df$`Pollinator species`%in%"Other" & abundy_plot_df$`Blueberry species`%in%"Rabbiteye"), c("Pollinator species")]="Tetragonula carbonaria" abundy_plot_df[(abundy_plot_df$`Pollinator species`%in%"Other" & abundy_plot_df$`Blueberry species`%in%"Southern highbush"), c("Pollinator species")]="Tetragonula carbonaria" abundy_plot_df[(abundy_plot_df$`Pollinator species`%in%"Other" & abundy_plot_df$`Blueberry species`%in%"Northern highbush"), c("Pollinator species")]="Bombus terrestris" abundy_plot_df$`Blueberry species` <- factor(abundy_plot_df$`Blueberry species`, levels = c("Northern highbush", "Rabbiteye","Southern highbush")) abundy_plot_df$`Pollinator species` <- revalue(abundy_plot_df$`Pollinator species`,c("Bombus terrestris"="Bumblebee", "Apis mellifera" = "Honeybee", "Tetragonula carbonaria"= "Stingless bee")) abundy_plot_df$pollsp2 <- abundy_plot_df$`Pollinator species` abundy_plot_df$pollsp2 <- revalue(abundy_plot_df$pollsp2,c("Honeybee" = "A", "Stingless bee" = "B", "Bumblebee" = "C")) #Making plot abun_cols=c("Bumblebee" = "#6a3d9a","C" = "#cab2d6", "Stingless bee" = "#33a02c","B" = "#b2df8a", "Honeybee" = "#1f78b4","A" = "#a6cee3") abundance_ggplot <- ggplot(abundy_plot_df,aes(x=`Pollinator species`,y=estimate,col=pollsp2))+ geom_text(aes(y=c(15),label = cld,col=`Pollinator species`),size=4, position = position_dodge(0.9),show.legend = FALSE, colour="black")+ theme_bw()+ ylim(0,15)+ geom_point(aes(y=Abundance/duration),pch = 1, position = position_jitterdodge(dodge.width = 1),alpha=1)+ geom_point(aes(col=`Pollinator species`),position = position_dodge(width = 1),size=2)+ geom_errorbar(aes(ymin = lower, ymax = upper,col=`Pollinator species`),pch=1, position = position_dodge(width = 0.5),width=0.1)+ theme(plot.title = element_text(hjust = 0.5), legend.box.background = element_rect(colour = "black"), legend.title = element_text(),aspect.ratio = 1, strip.background = element_blank(), text=element_text(), axis.ticks = element_blank(), axis.title.y =element_text(size=12), axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,angle = 45, hjust = 1), strip.text = element_text(size=14,family="Helvetica",face="bold"), panel.spacing = unit(0.25,"lines"))+ scale_color_manual(values=abun_cols)+ facet_wrap(~`Blueberry species`,scales="free_x")+ theme(legend.position = "none")+ ylab("Floral visitor abundance")+ xlab(NULL) abundance_ggplot