library(car) library(sjPlot) library(emmeans) library(lme4) library(tidyverse) #________________________________________________________----- #DISTANCE MODELS ----- #Read DPH BPH final dataset: ---- 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) #Model validation: par(mfrow=c(2,2)) #All ok plot(modist4.dph) plot(fitted(modist4.dph),residuals(modist4.dph)) qqnorm(resid(modist4.dph)) qqline(resid(modist4.dph)) # hist(resid(modist4.dph)) # acf(resid(modist4.dph)) # #Model results: summary(modist4.dph) Anova(modist4.dph) ##MODEL BPH Distance ---- dphdat3 = dphdat2 %>% filter(dph>0) modist4.bph<-glmer(cbind(BPH,nonBPH) ~ mindist +(1|location.id)+(1|juldayn), data=dphdat3,family=binomial(link="probit"),verbose=F, nAGQ=0) #Model validation: par(mfrow=c(2,2)) #All ok plot(modist4.bph) plot(fitted(modist4.bph),residuals(modist4.bph)) qqnorm(resid(modist4.bph)) qqline(resid(modist4.bph)) # hist(resid(modist4.bph)) # acf(resid(modist4.bph)) # #Model results: summary(modist4.bph) Anova(modist4.bph) #___________________________________________________________---- #STRUCTURE TYPE & DIEL MODELS ---- #Dataset including diel info: ---- dphdat.dn <- read.csv(here::here("Jacky Manuscript - Reef effect","Dryad","OFB_ReefEffect_HP-DPH-BPH_perDAYNIGHT_FinalDataset.csv"), header=TRUE) dphdat2.dn=dphdat.dn %>% filter(!is.na(closest.structure)) %>% #Removes CPODS > 1500m from structures filter(plus1pmaxreached==0) %>% #Removes days with >1% min max click mutate( date=as.Date(date), datef=as.factor(date), 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, juldayn=as.numeric(format(date,"%j"))) %>% ungroup() ##MODEL DPH StructureType: ---- mod2.dph.dn=glmer(cbind(dn.DPH,dn.nonDPH) ~ cpod.group*diel*closest.structure +(1|location.id)+(1|juldayn), data=dphdat2.dn,family=binomial(link="probit"),verbose=F, nAGQ=0) #logit, probit, cauchit Anova(mod2.dph.dn,type="III") #Interaction sig #Model validation: par(mfrow=c(2,2)) #All ok plot(mod2.dph.dn) plot(fitted(mod2.dph.dn),residuals(mod2.dph.dn)) qqnorm(resid(mod2.dph.dn)) qqline(resid(mod2.dph.dn)) # hist(resid(mod2.dph.dn)) # acf(resid(mod2.dph.dn)) # #Model results: Anova(mod2.dph.dn) summary(mod2.dph.dn) #PostHoc pairwise comparisons: emm.mod2.dph.dn = emmeans(mod2.dph.dn, c("cpod.group","closest.structure","diel"),by="closest.structure",type="response") pairs.emm.mod2.dph.dn=pairs(regrid(emm.mod2.dph.dn),type="response") #pvalues from here pairs.emm.mod2.dph.dn 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,diel,closest.structure,.group) ##MODEL BPH StructureType: ---- mod2.bph.dn=glmer(cbind(dn.BPH,dn.nonBPH) ~ cpod.group*diel*closest.structure +(1|location.id)+(1|juldayn), data=dphdat2.dn,family=binomial(link="probit"),verbose=F, nAGQ=0) Anova(mod2.bph.dn) #cpod.struc.group sig #Model validation: simulationOutput <- simulateResiduals(fittedModel = mod2.dph.dn,plot = T) #OK par(mfrow=c(2,2)) #All ok plot(mod2.bph.dn) plot(fitted(mod2.bph.dn),residuals(mod2.bph.dn)) qqnorm(resid(mod2.bph.dn)) qqline(resid(mod2.bph.dn)) # hist(resid(mod2.bph.dn)) # acf(resid(mod2.bph.dn)) # #Model results: Anova(mod2.bph.dn) summary(mod2.bph.dn) #PostHoc pairwise comparisons: emm.mod2.bph.dn = emmeans(mod2.bph.dn, c("cpod.group","closest.structure","diel"),by="closest.structure",type="response") pairs.emm.mod2.bph.dn=pairs(regrid(emm.mod2.bph.dn),type="response") #pvalues from here emm.letters.bph2=multcomp::cld(emm.mod2.bph.dn,Letters=letters) %>% select(cpod.struc.group,diel,.group) #___________________________________________________________---- #DIEL VARIATION - Pre-Post BeaDem installation ---- library(mgcv) library(DHARMa) library(ggplot2) library(sjPlot) library(tidyverse) #DATASET: 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") summary(gam1) plot.gam(gam1) #Check model: gam.check(gam1, rep=200) #k value OK simulationOutput <- simulateResiduals(fittedModel = gam1,plot = T) #OK ##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") summary(gam2) plot.gam(gam2) #Check model: gam.check(gam2, rep=200) #k value OK simulationOutput <- simulateResiduals(fittedModel = gam2,plot = T) #OK