#Fruit weight analyses of pollination experiments ####### library(bayesplot) library(brms) library(plyr) library(dplyr) library(emmeans) library(gridExtra) library(modelr) library(scales) library(tidyr) library(tidybayes) library(ggplot2) library(stringr) #dataframe #save as csv BLUEWGT <- read.csv("fruit weight pollination experiments JAPPL-2020-00197.csv") #set reference levels BLUEWGT$Species <- factor(BLUEWGT$Species,levels=c("Snowchaser","Blue Rose","Powder Blue")) BLUEWGT$TREATMENT <- factor(BLUEWGT$TREATMENT,levels=c("No pollination","Cross pollination","Open pollination","Self pollination")) ###MODEL #model terms #Fresh.wgt = Fresh fruit weight #Species = Blueberry type (Snowchaser = Southern highbush, Powder Blue = Rabbiteye, Blue Rose = Northern Highbush) #TREATMENT = Pollination treatment: Open, Cross, Self, No pollination #Year = sampling year #Block = Sampling block #RP = Plant ID within block (Plant Row _ Plant Number) #priors gam.priors <- prior(normal(0, 2),class="Intercept")+ prior(normal(0, 2),class="b")+ prior(normal(0,0.25),class="sd")+ prior(normal(0,0.5),class="sd",group="Year")+ prior(normal(0,1),dpar="shape")+ #northern highbush and rabbiteye prior(normal(0,2),class="Intercept",dpar="shape") #southern highbush BLUE.fruit.wgt=brm(bf(Fresh.wgt~Species*TREATMENT+ (1|Year/Block/RP),shape ~ Species, family=Gamma(link="log")), cores=4, prior=gam.priors, iter=2000, control=list(adapt_delta=0.999, max_treedepth=15), seed=123, data=BLUEWGT) pp_check(BLUE.fruit.wgt,nsamples=100) bayes_R2(BLUE.fruit.wgt) ############################### #####PAIRWISE DIFFERENCES##### ############################# ##FRUIT WEIGHT HDI PLOT fruit_wgt_hdi <- BLUE.fruit.wgt %>% emmeans( ~ Species*TREATMENT) %>% gather_emmeans_draws() %>% mean_hdci(.value=exp(.value)) %>% mutate(cld_group = paste(Species, TREATMENT, sep = "_")) #Compact letter differences blue_fw_pw <- BLUE.fruit.wgt %>% emmeans( ~ Species|TREATMENT) %>% gather_emmeans_draws() %>% mutate(intera = paste(Species,TREATMENT, sep = "_")) blue_fw_cld=blue_fw_pw %>% ungroup() %>% ## to get rid of unneeded columns select(.value, intera, .draw) %>% spread(intera, .value) %>% select(-.draw) %>% ## we need to get rid of all columns not containing draws cld_pmatrix() blue_fw_cld <- as.data.frame(blue_fw_cld) blue_fw_cld$cld_group <- rownames(blue_fw_cld) rownames(blue_fw_cld) <- NULL #percentage differences between treatments #northern highbush # (fruit_wgt_hdi[11,3]- fruit_wgt_hdi[10,3])/ fruit_wgt_hdi[11,3] ##28% #Rabbiteye #cross vs. open #rabbiteye (fruit_wgt_hdi[11,3]- fruit_wgt_hdi[10,3])/ fruit_wgt_hdi[11,3] ##28% #cross vs. self #rabbiteye (fruit_wgt_hdi[11,3]- fruit_wgt_hdi[12,3])/ fruit_wgt_hdi[11,3] ##32 #snowchaser #cross vs. self (fruit_wgt_hdi[3,3]- fruit_wgt_hdi[4,3])/ fruit_wgt_hdi[3,3] ##25% #cross vs. open (fruit_wgt_hdi[2,3]- fruit_wgt_hdi[4,3])/ fruit_wgt_hdi[2,3] ##26% #BR #cross vs. open (fruit_wgt_hdi[7,3]- fruit_wgt_hdi[8,3])/ fruit_wgt_hdi[7,3] ##20% #Plot dataframe fruit_wgt_plot_merge=merge(fruit_wgt_hdi,blue_fw_cld,by="cld_group") BLUEWGT3=BLUEWGT colnames(BLUEWGT3)[12]="Treatment" colnames(fruit_wgt_plot_merge)[3]="Treatment" fruit_wgt_plot=rbind.fill(fruit_wgt_plot_merge,BLUEWGT3) colnames(fruit_wgt_plot) fruit_wgt_plot$Species <- revalue(fruit_wgt_plot$Species,c("Blue Rose"="Northern highbush", "Powder Blue" = "Rabbiteye", "Snowchaser" = "Southern highbush")) fruit_wgt_plot$Species=factor(fruit_wgt_plot$Species,levels=c("Northern highbush","Rabbiteye", "Southern highbush")) table(fruit_wgt_plot$Treatment) fruit_wgt_plot$Treatment=factor(fruit_wgt_plot$Treatment,levels=c("No pollination","Cross pollination", "Open pollination","Self pollination")) fruit_wgt_plot$Treatment2 <- fruit_wgt_plot$Treatment fruit_wgt_plot$Treatment2 <- revalue(fruit_wgt_plot$Treatment2,c("No pollination" = "A", "Cross pollination" = "B", "Open pollination" = "C", "Self pollination" = "D")) fw_plot_cols=c("A" = "#a6cee3", "No pollination" = "#1f78b4", "B" = "#b2df8a", "Cross pollination" = "#33a02c", "C" = "#fb9a99", "Open pollination" = "#e31a1c", "D" = "#cab2d6", "Self pollination" ="#6a3d9a") fruit_wgt_overall_plot <- ggplot(fruit_wgt_plot,aes(col=Treatment2,y=.value,x=Treatment))+ geom_text(aes(y=4.2,label = blue_fw_cld), color="black", position = position_dodge(0.9),show.legend = FALSE)+ geom_point(aes(x=Treatment,y=Fresh.wgt,col=Treatment2), position = position_jitterdodge(dodge.width = 1.2,jitter.width = 0.5,jitter.height=0.1), alpha=0.5, pch=1)+ theme_bw()+ geom_point(aes(col=Treatment),position = position_dodge(width=1),size=2.5)+ geom_errorbar(aes(ymin=.lower,ymax=.upper,col=Treatment),position = position_dodge(width=1),width=0.3)+ scale_color_manual(values=fw_plot_cols,guide=FALSE)+ 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=14), 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"))+ facet_wrap(~Species,scales="free_x")+ labs(x = NULL, y = "Fruit weight (g)", color = "Treatment") fruit_wgt_overall_plot