#___________________________________________________________---- #Figure2: GLMM Distance ---- library(car) library(sjPlot) library(emmeans) library(lme4) library(tidyverse) library(cowplot) dphdat <- read.csv(here::here("Jacky Manuscript - Reef effect","Dryad","OFB_ReefEffect_HP-DPH-BPH_perDAY_FinalDataset.csv") ,header=TRUE) dphdat2=dphdat %>% filter(plus1pmaxreached==0) %>% #Removes days with >1% min max click mutate( date=as.Date(date), datef=as.factor(date), DPH=as.integer(dph), nonDPH=total.hours-dph, BPH=as.integer(bph), nonBPH=dph-bph, propDPH=DPH/total.hours, propBPH=BPH/DPH, total.hours=as.integer(total.hours), julday=format(date,"%j"), juldayn=as.numeric(format(date,"%j"))) ##MODEL DPH Distance: ---- modist4.dph=glmer(cbind(DPH,nonDPH) ~ mindist +(1|juldayn) + (1|location.id), data=dphdat2,family=binomial(link="probit"),verbose=F, nAGQ=0) Anova(modist4.dph,type = "II") #*** #Plot model: predmindist.modist4=jtools::make_predictions(modist4.dph,pred="mindist",data=jackydata) %>% mutate(label="***") colnames(predmindist.modist4)[6]="propDPH" plot.midist4=ggplot() + geom_point(data=dphdat2,aes(x=mindist,y=propDPH), position=position_jitter(width = 20,height = 0.025),alpha=0.4, shape=21, col="#6a7987",fill="#9dadbc" )+ geom_ribbon(data=predmindist.modist4,aes(x=mindist,ymin=ymin,ymax=ymax),fill="#005c86",alpha=0.4) + geom_line(data=predmindist.modist4,aes(x=mindist,y=propDPH),colour="#005c86",size=1) + geom_text(data=predmindist.modist4, aes(x=1400,y=0.9,label=label),size=10)+ labs(title = "", x="Distance to offshore structures (m)", y="p(occurrence)") + scale_y_continuous(limit=c(0,1),breaks=seq(0.2,0.8,by=0.2),expand = c(0,0))+ scale_x_continuous(limit=c(0,2750),breaks=seq(0,2500,by=500),expand = c(0,0))+ theme_bw()+ theme(axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), axis.title.x = element_blank(), axis.title.y = element_text(size=12), plot.title = element_text(face = "bold",size=11, hjust = 0.5), strip.text = element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "right") plot.midist4 ##MODEL BPH Distance to offshore structures ---- dphdat3 = dphdat2 %>% filter(dph>0) modist4.bph<-glmer(cbind(BPH,nonBPH) ~ mindist +(1|location.id)+(1|juldayn), data=dphdat3,family=binomial(link="log"),verbose=F, nAGQ=0) Anova(modist4.bph, type="II") #*** #Plot model: predmindist.modist4buzz=jtools::make_predictions(modist4.bph,pred="mindist",data=dphdat3) %>% mutate(label="***") colnames(predmindist.modist4buzz)[6]="propBPH" per.increase=(0.51-0.35)/0.35*100 per.increase #22% higher closer to structures plot.midist4.buzz=ggplot() + geom_point(data=dphdat3,aes(x=mindist,y=propBPH), position=position_jitter(width = 20,height = 0.025),alpha=0.4, shape=21, col="#6a7987",fill="#9dadbc" )+ geom_ribbon(data=predmindist.modist4buzz,aes(x=mindist,ymin=ymin,ymax=ymax),fill="#7f4362",alpha=0.4) + geom_line(data=predmindist.modist4buzz,aes(x=mindist,y=propBPH),colour="#7f4362",size=1) + geom_text(data=predmindist.modist4buzz, aes(x=1400,y=0.9,label=label),size=10)+ labs(title = "", x="Distance to offshore structures (m)", y="p(foraging)") + scale_y_continuous(limit=c(0,1),breaks=seq(0.2,0.8,by=0.2),expand = c(0,0))+ scale_x_continuous(limit=c(0,2750),breaks=seq(0,2500,by=500),expand = c(0,0))+ theme_bw()+ theme(axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), axis.title.x = element_blank(), axis.title.y = element_text(size=12), plot.title = element_text(face = "bold",size=11, hjust = 0.5), strip.text = element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "right") plot.midist4.buzz ##Save Distance Models: ---- medplots= ggdraw() + draw_plot(plot.midist4, x = 0, y = 0, width = .5, height = 1) + draw_plot(plot.midist4.buzz, x = .5, y = 0, width = .5, height = 1) + draw_plot_label(label = c("A", "B"), size = 15, x = c(0, 0.5), y = c(1, 1)) medplots newmed=ggpubr::annotate_figure(medplots, bottom=grid::textGrob("Distance to offshore structures (m)", gp = grid::gpar(cex = 1.1))) ggsave(newmed, file=here::here("Jacky Manuscript - Reef effect","MS Plots","Fig2.png"), width = 18, height = 9, units = "cm") #___________________________________________________________________---- #Figure3: Diel/StructureType ---- library(car) library(sjPlot) library(emmeans) library(lme4) library(tidyverse) library(cowplot) all2021PODdata.dn <- read.csv(here::here("Jacky Manuscript - Reef effect","Dryad","OFB_ReefEffect_OFB_ReefEffect_HP-DPH-BPH_perDAYNIGHT_FinalDataset.csv"),header=TRUE) jackydata.dn=all2021PODdata.dn %>% filter(plus1pmaxreached==0) %>% #Removes days with >1% min max click filter(!is.na(closest.structure)) %>% #Removes CPODS > 1500m from structures mutate( date=as.Date(date), closest.structure=ifelse(closest.structure=="WTA/B","Beatrice Demonstrator",closest.structure), dn.DPH=dn.dph, dn.nonDPH=dn.toth-dn.dph, dn.BPH=dn.bph, dn.nonBPH=dn.dph-dn.bph, dn.propDPH=dn.DPH/dn.toth, dn.propBPH=dn.BPH/dn.DPH, julday=format(date,"%j"), juldayn=as.numeric(format(date,"%j"))) mindist.min=min(jackydata.dn$mindist) mindist.max=max(jackydata.dn$mindist) #Data subset: closest CPODs to the platforms (<1500m): cpodsub.dn=jackydata.dn %>% filter(cpod.group=="Non-structure\nCPOD" | cpod.group =="Structure\nCPOD") #Platform CPODs (<200m) and nearest CPODs (<1500m) ##MODEL DPH PlatformType: ---- mod2.dph.dn=glmer(cbind(dn.DPH,dn.nonDPH) ~ cpod.group*diel*closest.structure +(1|location.id)+(1|juldayn), data=cpodsub.dn,family=binomial(link="probit"),verbose=F, nAGQ=0) Anova(mod2.dph.dn,type="III") #Matrix PostHoc pairwise comparisons: mod2B.dph.dn=glmer(cbind(dn.DPH,dn.nonDPH) ~ cpod.struc.group*diel +(1|location.id)+(1|juldayn), data=cpodsub.dn,family=binomial(link="probit"),verbose=F, nAGQ=0) emm.mod2B.dph.dn = emmeans(mod2B.dph.dn, ~ cpod.struc.group|diel,type="response") pairs.emm.mod2.dph.dn=pairs(regrid(emm.mod2B.dph.dn),type="response") #pvalues from here matrix.emm.dph=pwpm(regrid(emm.mod2B.dph.dn), by=NULL,digits = 3) matrix.emm.dph write.csv(matrix.emm.dph,here::here("Jacky Manuscript - Reef effect","MS Plots","TableS3_DPH_NoFormated.csv"),row.names = T) #Plot model: #New data: cpod.group = c("Non-structure\nCPOD","Structure\nCPOD") closest.structure=c("Jacky","Beatrice Bravo","Beatrice Demonstrator") diel = c("day", "night") predData2=as.data.frame(expand.grid(cpod.group,closest.structure,diel)) %>% setNames(c("cpod.group","closest.structure","diel")) #Predictions on new data: predmod2.dn=jtools::make_predictions(mod2.dph.dn,new_data=predData2) colnames(predmod2.dn)[4]="propDPH" #Add posthoc comparisons: emm.mod2.dph.dn.l = emmeans(mod2.dph.dn, c("cpod.group","closest.structure","diel"),type="response") emm.letters.dph2=multcomp::cld(emm.mod2.dph.dn.l,Letters=letters) %>% select(cpod.group,closest.structure,diel,.group) %>% mutate(let.group=gsub("[[:space:]]", "", .group), let.group= case_when(nchar(let.group)<5 ~ let.group, nchar(let.group)>=5 ~ paste0(substr(let.group,1,4),"\n",substr(let.group,5,nchar(let.group))))) predmod2.dn=left_join(predmod2.dn,emm.letters.dph2) %>% mutate( closest.structure=factor(closest.structure,levels=c("Jacky","Beatrice Bravo","Beatrice Demonstrator"),ordered = T), cpod.group=factor(cpod.group,levels=c("Non-structure\nCPOD","Structure\nCPOD"),ordered = T)) #Order raw data: cpodsub.dn.ord=cpodsub.dn %>% mutate( closest.structure=factor(closest.structure,levels=c("Jacky","Beatrice Bravo","Beatrice Demonstrator"),ordered = T), cpod.group=factor(cpod.group,levels=c("Non-structure\nCPOD","Structure\nCPOD"),ordered = T)) plot.dph.mod2=ggplot() + geom_point(data=cpodsub.dn.ord,aes(x=cpod.group,y=dn.propDPH,fill=diel), position=position_jitterdodge(jitter.width = 0.1,jitter.height = 0.02,dodge.width = 0.7), alpha=0.05, shape=21, col="#6a7987" ) + geom_errorbar(data=predmod2.dn,aes(x=cpod.group,ymin=ymin,ymax=ymax,colour=diel), lwd=1,width=0.2, position=position_dodge(width = 0.7)) + geom_point(data=predmod2.dn,aes(x=cpod.group,y=propDPH,fill=diel,shape=diel), size=2, position=position_dodge(width = 0.7)) + geom_text(data=predmod2.dn,aes(x=cpod.group,y=0,label=let.group,colour=diel,group=diel), size=3, colour="black", position=position_dodge(width = 0.7)) + labs(title = "", x="", y="p(occurrence)") + scale_shape_manual(values=c(21,23))+ scale_y_continuous(limit=c(-0.05,1),breaks=seq(0,1,by=0.2),expand = c(0,0.05))+ scale_colour_manual(values = c("#e69c24","#005c86")) + scale_fill_manual(values = c("#e69c24","#005c86")) + facet_wrap(.~closest.structure,scales="free_x") + theme_bw()+ theme(axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), axis.title.x = element_blank(), axis.title.y = element_text(size=12), plot.title = element_text(face = "bold",size=11, hjust = 0.5), strip.text = element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") plot.dph.mod2 ##Model BPH PlatformType: ---- cpodsub.dn.buzz=cpodsub.dn %>% filter(dn.dph>0) mod2.bph.dn=glmer(cbind(dn.BPH,dn.nonBPH) ~ cpod.group*diel*closest.structure +(1|location.id)+(1|juldayn), data=cpodsub.dn.buzz,family=binomial(link="probit"),verbose=F, nAGQ=0) Anova(mod2.bph.dn,type="III") #Matrix PostHoc pairwise comparisons: mod2B.bph.dn=glmer(cbind(dn.BPH,dn.nonBPH) ~ cpod.struc.group*diel +(1|location.id)+(1|juldayn), data=cpodsub.dn.buzz,family=binomial(link="probit"),verbose=F, nAGQ=0) emm.mod2B.bph.dn = emmeans(mod2B.bph.dn, ~ cpod.struc.group|diel,type="response") pairs.emm.mod2B.bph.dn=pairs(regrid(emm.mod2B.bph.dn),type="response") #pvalues from here matrix.emm=pwpm(regrid(emm.mod2B.bph.dn), by=NULL,digits = 3) matrix.emm write.csv(matrix.emm,here::here("Jacky Manuscript - Reef effect","MS Plots","TableS4_BPH_NoFormated.csv"),row.names = T) #Plot model: #New data: cpod.group = c("Non-structure\nCPOD","Structure\nCPOD") closest.structure=c("Jacky","Beatrice Bravo","Beatrice Demonstrator") diel = c("day", "night") predData2=as.data.frame(expand.grid(cpod.group,closest.structure,diel)) %>% setNames(c("cpod.group","closest.structure","diel")) #Predict on new data: predmod2buzz.dn=jtools::make_predictions(mod2.bph.dn,new_data=predData2) #Add posthoc comparison: emm2.mod2.bph.dn.l = emmeans::emmeans(mod2.bph.dn,c("cpod.group","diel","closest.structure"),type="response") emm.letters.bph2.l=multcomp::cld(emm2.mod2.bph.dn.l,Letters=letters) %>% select(cpod.group,diel,closest.structure,.group) colnames(predmod2buzz.dn)[4]="propBPH" predmod2buzz.dn=left_join(predmod2buzz.dn,emm.letters.bph2.l) %>% mutate( let.group=gsub("[[:space:]]", "", .group), closest.structure=factor(closest.structure,levels=c("Jacky","Beatrice Bravo","Beatrice Demonstrator"),ordered = T), cpod.group=factor(cpod.group,levels=c("Non-structure\nCPOD","Structure\nCPOD"),ordered = T)) #Order raw data: cpodsub.dn.buzz.ord=cpodsub.dn.buzz %>% mutate( closest.structure=factor(closest.structure,levels=c("Jacky","Beatrice Bravo","Beatrice Demonstrator"),ordered = T), cpod.group=factor(cpod.group,levels=c("Non-structure\nCPOD","Structure\nCPOD"),ordered = T)) #Plot: plot.bph.mod2=ggplot() + geom_point(data=cpodsub.dn.buzz.ord,aes(x=cpod.group,y=dn.propBPH,fill=diel), position=position_jitterdodge(jitter.width = 0.1,jitter.height = 0.02,dodge.width = 0.7), alpha=0.05, shape=21, col="#6a7987" ) + geom_errorbar(data=predmod2buzz.dn,aes(x=cpod.group,ymin=ymin,ymax=ymax,colour=diel), lwd=1,width=0.2, position=position_dodge(width = 0.7)) + geom_point(data=predmod2buzz.dn,aes(x=cpod.group,y=propBPH,fill=diel,shape=diel), size=2, position=position_dodge(width = 0.7)) + geom_text(data=predmod2buzz.dn,aes(x=cpod.group,y=0,label=let.group,group=diel), size=3, colour="black", position=position_dodge(width = 0.7)) + labs(title = "", x="", y="p(foraging)") + scale_shape_manual(values=c(21,23))+ scale_y_continuous(limit=c(-0.05,1),breaks=seq(0,1,by=0.2),expand = c(0,0.05))+ scale_colour_manual(values = c("#e69c24","#7f4362")) + scale_fill_manual(values = c("#e69c24","#7f4362")) + facet_wrap(.~closest.structure,scales="free_x") + theme_bw()+ theme(axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), axis.title.x = element_text(size=12), axis.title.y = element_text(size=12), plot.title = element_text(face = "bold",size=11, hjust = 0.5), strip.text = element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") plot.bph.mod2 ##Save Distance plots: ---- displots= ggdraw() + draw_plot(plot.dph.mod2, x =0, y = 0.5, width = 1, height =0.5) + draw_plot(plot.bph.mod2, x = 0, y = 0, width = 1, height = 0.5) + draw_plot_label(label = c("A", "B"), size = 15, x = c(0, 0), y = c(1, 0.5)) displots ggsave(displots, file=here::here("Jacky Manuscript - Reef effect","MS Plots","Fig3.png"), width = 18, height = 15, units = "cm") #_______________________________________________---- #Figure4: Diel comparison Pre-Post installation ---- library(mgcv) library(sjPlot) dicomp=read.csv(here::here("Jacky Manuscript - Reef effect","Dryad","OFB_ReefEffect_HP-PPM-BPM_Pre-Post-FinalDataset.csv"), header = T) %>% mutate( por.pres=as.factor(por.pres), depl=as.factor(depl), year=as.factor(year), structure=as.factor(structure), buzz.pres=as.factor(ifelse(bpm>0,1,0))) ##GAM MODEL: Hourly Presence: ---- gam1=gam(por.pres ~ s(hour,k=4, bs="cc",by=structure) + structure, data=dicomp, family=binomial(link = "probit"), method = "REML") #Plot model: hour = seq(0,23.5,by=0.5) structure=as.factor(c("Present","Absent")) predData=as.data.frame(expand.grid(hour,structure)) %>% setNames(c("hour","structure")) predgam1=jtools::make_predictions(gam1, new_data = predData, outcome.scale="response") dicomp2=dicomp %>% mutate( por.pres=as.numeric(as.character(por.pres)), por.pres=case_when(structure=="Absent" ~ por.pres, structure=="Present" ~ por.pres-0.03 ) ) plot.gam1=ggplot() + geom_point(data=dicomp2,aes(x=hour,y=por.pres,shape=structure,colour=structure), size=1,position = position_jitter(width=0.4,height = 0)) + geom_ribbon(data=predgam1,aes(x=hour,ymin=ymin,ymax=ymax,fill=structure), alpha=0.2) + geom_line(data=predgam1,aes(x=hour,y=por.pres,colour=structure,linetype=structure), size=1) + labs(title = "", x="Time of day", y="p(occurrence)") + scale_linetype_manual(values=c("dashed","solid"))+ scale_shape_manual(values=c(1,2))+ scale_y_continuous(limits=c(-0.05,1.02),breaks=seq(0,1,by=0.2),expand = c(0,0))+ scale_x_continuous(limits = c(0,23.5), breaks = seq(0,21,by=3),labels = paste0(seq(0,21,by=3),"h"),expand = c(0.01,0)) + scale_colour_manual(values = c("#9dadbc","#005c86")) + scale_fill_manual(values = c("#9dadbc","#005c86")) + theme_bw()+ theme(axis.text.x = element_text(size=10), axis.text.y = element_text(size=10), axis.title.y = element_text(size=12), axis.title.x = element_blank(), plot.title = element_text(face = "bold",size=11, hjust = 0.5), strip.text = element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="none") plot.gam1 ##GAM MODEL: Hourly Foraging: ---- dicompbuzz=dicomp %>% filter(por.pres!="0") gam2=gam(buzz.pres ~ s(hour,k=4, bs="cc",by=structure) + structure, data=dicompbuzz, family=binomial(link = "probit"), method = "REML") #Plot model: hour = seq(0,23.5,by=0.5) structure=as.factor(c("Present","Absent")) predData=as.data.frame(expand.grid(hour,structure)) %>% setNames(c("hour","structure")) predgam2=jtools::make_predictions(gam2, new_data = predData, outcome.scale="response") dicomp2b=dicompbuzz %>% mutate( por.pres=as.numeric(as.character(por.pres)), buzz.pres=as.numeric(as.character(buzz.pres)), buzz.pres=case_when(structure=="Absent" ~ buzz.pres, structure=="Present" ~buzz.pres-0.03 ) ) plot.gam2=ggplot() + geom_ribbon(data=predgam2,aes(x=hour,ymin=ymin,ymax=ymax,fill=structure), alpha=0.2) + geom_line(data=predgam2,aes(x=hour,y=buzz.pres,colour=structure,linetype=structure), size=1) + geom_point(data=dicomp2b,aes(x=hour,y=buzz.pres,shape=structure,colour=structure), size=1,position = position_jitter(width=0.4,height = 0)) + labs(title = "", x="Time of day", y="p(foraging)") + scale_linetype_manual(values=c("dashed","solid"))+ scale_shape_manual(values=c(1,2))+ scale_y_continuous(limits=c(-0.05,1.02),breaks=seq(0,1,by=0.2),expand = c(0,0))+ scale_x_continuous(limits = c(0,23.5), breaks = seq(0,21,by=3),labels = paste0(seq(0,21,by=3),"h"),expand = c(0.01,0)) + scale_colour_manual(values = c("#9dadbc","#7F4362")) + scale_fill_manual(values = c("#9dadbc","#7F4362")) + theme_bw()+ theme(axis.text.x = element_text(size=10), axis.text.y = element_text(size=10), axis.title.y = element_text(size=12), axis.title.x = element_blank(), plot.title = element_text(face = "bold",size=11, hjust = 0.5), strip.text = element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="none") plot.gam2 ##Save GAM plots: ---- gamplots= ggdraw() + draw_plot(plot.gam1, x =0, y = 0, width = 0.5, height =1) + draw_plot(plot.gam2, x = 0.5, y = 0, width = 0.5, height = 1) + draw_plot_label(label = c("A", "B"), size = 15, x = c(0, 0.5), y = c(1, 1)) gamplots ggsave(gamplots, file=here::here("Jacky Manuscript - Reef effect","MS Plots","Fig4.png"), width = 18, height = 9, units = "cm")