### R-Script 'Honeybees possess a structurally diverse and functionally redundant set of queen pheromones' ### Proceedings of the Royal Society (2019) ### by Sarah A. Princen, Ricardo Caliari Oliveira, Ulrich R. Ernst, Jocelyn G. Millar, Jelle S. van Zweden & Tom Wenseleers ### Analysis of queen pheromone treatment on the proportion of workers with developed ovaries (Fig. 2) # load required packages library(devtools) library(afex) library(lme4) library(car) library(effects) library(doBy) library(lattice) library(emmeans) library(ggplot2) library(ggthemes) library(export) library(gtools) # read in data data=read.csv("Princen_etal_datafig2.csv") data$Treatment=factor(data$Treatment,levels=c("Acetone","9-ODA","9-HDA","HOB","HVA","Cut_esters","Cut_alkenes","TG_esters","TG_acids","QMP")) head(data) # Variables: # Nextbox = nestbox ID # Replicate = replicate nr. # Worker = individual worker ID # Treatment = pheromone treatments (Acetone = solvent-only negative control) # Ovaries_developed = whether or not each worker had developed ovaries # Analysis of effects of pheromone treatments on proportion of workers with developed ovaries #### set_sum_contrasts() # we use sum contrast coding (effect coding) fit=glmer(Ovaries_developed~(1|Replicate)+(1|Nestbox)+Treatment, family=binomial(link=logit), data=data, control=glmerControl(optimizer="bobyqa")) AIC(fit) # 1600.495, better than random slope model with (Treatment|Replicate) (AIC=1665.40) summary(fit) # Fixed effects: # Estimate Std. Error z value Pr(>|z|) # (Intercept) -0.7554 0.3940 -1.918 0.055173 . # Treatment9-ODA -1.8754 0.4657 -4.027 5.65e-05 *** # Treatment9-HDA -2.1011 0.4792 -4.385 1.16e-05 *** # TreatmentHOB -0.8473 0.4661 -1.818 0.069091 . # TreatmentHVA -1.1649 0.4429 -2.630 0.008528 ** # TreatmentCut_esters -0.8179 0.4134 -1.978 0.047873 * # TreatmentCut_alkenes -0.8126 0.4143 -1.961 0.049827 * # TreatmentTG_esters -1.9417 0.4741 -4.095 4.21e-05 *** # TreatmentTG_acids -0.9188 0.4199 -2.188 0.028640 * # TreatmentQMP -1.6844 0.4580 -3.678 0.000235 *** Anova(fit, type="III") # Anova type III tests, to test overall effect of pheromone treatments # Analysis of Deviance Table (Type III Wald chisquare tests) # # Response: Ovaries_developed # Chisq Df Pr(>Chisq) # (Intercept) 3.6769 1 0.05517 . # Treatment 36.0360 9 3.907e-05 *** # posthoc tests against neg control using emmeans df_posthoc=data.frame(summary(contrast(emmeans(fit,~Treatment, type="response"), method="trt.vs.ctrl1", group="Acetone", adjust="none"))) # contrasts with acetone control df_posthoc$p.value=df_posthoc$p.value/2 # we use 1-sided FDR corrected p values df_posthoc$p.fdradj=p.adjust(df_posthoc$p.value, method="fdr") # FDR p value correction df_posthoc$sign=as.vector(stars.pval(df_posthoc$p.fdradj)) # significances shown in Fig. 2 df_posthoc$odds.ratio.lower95CB=df_posthoc$odds.ratio+1.96*df_posthoc$SE df_posthoc$odds.ratio.upper95CB=df_posthoc$odds.ratio-1.96*df_posthoc$SE df_posthoc$odds.ratio=1/df_posthoc$odds.ratio # to express odds ratio for workers to have developed ovaries in control relative to in each treatment, ie to reverse contrast compared to default df_posthoc$odds.ratio.lower95CB=1/df_posthoc$odds.ratio.lower95CB df_posthoc$odds.ratio.upper95CB=1/df_posthoc$odds.ratio.upper95CB df_posthoc=df_posthoc[c(9,1,2,3,4,5,6,7,8),c(1,2,9,10,5,6,7,8)] df_posthoc$contrast=gsub(" / Acetone", "", df_posthoc$contrast) colnames(df_posthoc)[1]="treatment" rownames(df_posthoc)=NULL write.csv(df_posthoc,"Table S2_posthoc tests.csv") df_posthoc # treatment odds.ratio odds.ratio.lower95CB odds.ratio.upper95CB z.ratio p.value p.fdradj sign # 1 QMP 5.388996 2.839523 52.75709 -3.676933 1.180275e-04 2.655620e-04 *** # 2 9-ODA 6.523348 3.410594 74.69971 -4.027469 2.819020e-05 8.457060e-05 *** # 3 9-HDA 8.175139 4.215498 134.69346 -4.384251 5.819276e-06 5.237349e-05 *** # 4 HOB 2.333373 1.219253 27.06060 -1.817452 3.457398e-02 3.457398e-02 * # 5 HVA 3.205680 1.715850 24.33624 -2.629639 4.273776e-03 7.692797e-03 ** # 6 Cut_esters 2.265720 1.251540 11.94663 -1.978252 2.395014e-02 2.805592e-02 * # 7 Cut_alkenes 2.253658 1.243649 11.99609 -1.961016 2.493859e-02 2.805592e-02 * # 8 TG_esters 6.970385 3.612851 98.63403 -4.095069 2.110211e-05 8.457060e-05 *** # 9 TG_acids 2.506324 1.374762 14.16780 -2.187933 1.433723e-02 2.150585e-02 * # effect plot of treatment effects on worker ovary development df=data.frame(effect("Treatment",fit,type="response", confidence.level=0.9)) # we use one-sided 95% confidence intervals, ie 95% confidence bounds to match our use of 1-sided p values df$label = c("",df_posthoc$sign[1:(nrow(df)-1)]) # for significance labels above bars df$Treatment=factor(df$Treatment,levels=c("Acetone","QMP","9-ODA","9-HDA","HOB","HVA","Cut_esters","Cut_alkenes","TG_esters","TG_acids"), labels=c("CONTROL","QMP","9-ODA","9-HDA","HOB","HVA","CUT_ESTERS","CUT_ALKENES","TG_ESTERS","TG_ACIDS")) df.label = df df.label$fit = df.label$upper df # Treatment fit se lower upper label # 1 CONTROL 0.31964120 0.08569767 0.19724531 0.4732153 # 2 9-ODA 0.06718174 0.02881851 0.03269767 0.1330322 *** # 3 9-HDA 0.05434532 0.02433015 0.02569992 0.1112727 *** # 4 HOB 0.16759956 0.06443942 0.08607801 0.3009074 * # 5 HVA 0.12782307 0.04869298 0.06668470 0.2311333 ** # 6 CUT_ESTERS 0.17174448 0.05768726 0.09618368 0.2877655 * # 7 CUT_ALKENES 0.17250513 0.05805590 0.09648149 0.2892547 * # 8 TG_ESTERS 0.06314519 0.02769172 0.03026458 0.1270677 *** # 9 TG_ACIDS 0.15785992 0.05483712 0.08684888 0.2697789 * # 10 QMP 0.08018911 0.03333837 0.03980131 0.1549465 *** colcontrol <<- "grey" colsQMP <<- colorRampPalette(c("firebrick2", "orange"))(5) colsCUT <<- c("steelblue2","steelblue") # c("springgreen3","springgreen4") colsTG <<- c("slateblue2","slateblue4") pl = qplot(x=Treatment, y=fit*100, fill=Treatment, data=df, geom="col", width=0.5, xlab="", ylab="WORKERS WITH DEVELOPED OVARIES (%)") + geom_linerange(aes(x=Treatment, ymin=lower*100, ymax=upper*100)) + theme_few() + theme(legend.position="none") + theme(axis.text.x=element_text(color="black", size=11, angle=-30, vjust=0.5, hjust=0.1)) + scale_fill_manual(values=c(colcontrol,colsQMP,colsCUT,colsTG)) + geom_text(data = df.label, aes(label=label), nudge_y=1.5, vjust=0.5, size=5) pl graph2ppt(pl, file="Fig2_effect plot_binomial glmm_treatments.pptx", width=7, height=5)