# Statistical analyses of herbivorous fish abundance, herbivory rates and defecation rates outlined in # Zarco-Perello et al. Overwintering tropical herbivores accelerate detritus production on temperate reefs, # Proc Roy Soc B. #Packages ---- # Graphics install.packages("sjPlot", dependencies=TRUE) library(sjPlot) install.packages("sjmisc") library(sjmisc) install.packages("ggiraphExtra", dependencies=TRUE) library(ggiraphExtra) install.packages("reshape2") library(reshape2) install.packages("ggplot2") library("ggplot2") install.packages("ggpmisc") library(ggpmisc) install.packages("gridExtra") library(gridExtra) install.packages("grid") library(grid) # Data manipulation packages install.packages("plyr", dependencies= TRUE) library(plyr) install.packages("tidyr", dependencies= TRUE) library(tidyr) install.packages("dplyr") library(dplyr) install.packages("dlookr", dependencies= TRUE) library(dlookr) # Multipanel Scatterplots install.packages("PerformanceAnalytics", dependencies = TRUE) library("PerformanceAnalytics") ## Modeling Packages ---- install.packages("quantreg", dependencies = TRUE) # quantile regression library("quantreg") install.packages("lqmm") # mixed effects quantile regression library("lqmm") install.packages("nlme") # mixed effects linear regression library(nlme) install.packages("sjstats", dependencies = TRUE) # Summary and fit statistics for G-LMM library(sjstats) install.packages("glmmTMB") # generalized mixed effects regression library(glmmTMB) install.packages("DHARMa", dependencies = TRUE) # Diagnostics for GLMM library(DHARMa) #.---- # Main dataset ---- Herbivory <- read.csv("Full Dataset Weighted Variables.csv") str(Herbivory) Herbivory$Year <- as.factor(Herbivory$Year) # Fish Abundance ##### ---- #.---- ## Data Exploration ---- Fish.abund <-Herbivory%>% select(Siganus.DOV, Siganus.MaxN, Sydneyanus.DOV, Sydneyanus.MaxN, Temperate.DOV, Temperate.MaxN, Temperature) str(Fish.abund) chart.Correlation(Fish.abund, histogram=TRUE, pch=19) Sig.abund<-Herbivory%>% select(Temperature, Siganus.MaxN, Siganus.DOV) chart.Correlation(Sig.abund, histogram=TRUE, pch=19) Temp.abund<-Herbivory%>% select( Temperature, Temperate.MaxN, Temperate.DOV) chart.Correlation(Temp.abund, histogram=TRUE, pch=19) ggplot(data=Herbivory, aes(y=Siganus.DOV, x= Temperature, color=Year))+ geom_point()+ geom_smooth(method=lm) ggplot(data=Herbivory, aes(y=Temperate.DOV, x= Temperature, color=Year))+ geom_point()+ geom_smooth(method=lm) ## Fish data ### ---- Fish.abund.r <-Herbivory%>% select(Siganus.DOV, Siganus.MaxN, Sydneyanus.DOV, Sydneyanus.MaxN, Temperate.DOV, Temperate.MaxN, Temperature, Reef, Year) str(Fish.abund.r) # Siganus MaxN ---- # Generalized Linear Mixed Models GLMMtmb ---- ctrl= glmmTMBControl(optCtrl = list(iter.max = 1000, eval.max = 1000), profile = FALSE, collect = FALSE) glmm1.Siganus.MaxN=glmmTMB(Siganus.MaxN ~ Temperature + (1|Reef),family=nbinom1, control=ctrl, data=Fish.abund.r) glmm2.Siganus.MaxN=glmmTMB(Siganus.MaxN ~ Temperature + (1|Reef) + (1|Year), family=nbinom1, control=ctrl, data=Fish.abund.r) glmm3.Siganus.MaxN=glmmTMB(Siganus.MaxN ~ Temperature + (Temperature|Reef), family=poisson, control=ctrl, data=Fish.abund.r) glmm4.Siganus.MaxN=glmmTMB(Siganus.MaxN ~ Temperature + (Temperature|Reef/Year), family=nbinom2, control=ctrl, data=Fish.abund.r) AIC(glmm1.Siganus.MaxN, glmm2.Siganus.MaxN, glmm3.Siganus.MaxN, glmm4.Siganus.MaxN) # df AIC #glmm1.Siganus.MaxN 4 248.1241 #glmm2.Siganus.MaxN 5 249.1230 #glmm3.Siganus.MaxN 5 689.7900 #glmm4.Siganus.MaxN 9 NA anova(glmm1.Siganus.MaxN, glmm2.Siganus.MaxN) # No sig differences # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq) #glmm1.Siganus.MaxN 4 248.12 254.34 -120.06 240.12 #glmm2.Siganus.MaxN 5 249.12 256.90 -119.56 239.12 1.001 1 0.3171 simulationOutput <- simulateResiduals(fittedModel = glmm2.Siganus.MaxN, n=250 ) # Checking residuals plot(simulationOutput) # Good fit testDispersion(simulationOutput) summary(glmm2.Siganus.MaxN) # No sig differences # Estimate Std. Error z value Pr(>|z|) #(Intercept) 1.01196 1.64073 0.617 0.537 #Temperature 0.07737 0.07905 0.979 0.328 # Mixed-Effects Quantile Regression ---- fit.Siganus.MaxN$control$LP_tol_ll <- 1e-3 fit.Siganus.MaxN$control$LP_max_iter <- 1000 fit.Siganus.MaxN <- lqmm(Siganus.MaxN ~ Temperature, random = ~ 1, group = Reef, covariance = "pdDiag", tau = c(0.90), nK = 7, type = "normal", data = Fish.abund.r, control = lqmmControl(method = "df")) summary(fit.Siganus.MaxN, R = 200, seed = 52) # No sig differences #tau = 0.9 #Fixed effects: # Value Std. Error lower bound upper bound Pr(>|t|) #(Intercept) 2.3751 124.0134 -242.1740 246.924 0.9847 #Temperature 2.5000 7.2323 -11.7618 16.762 0.7300 #AIC: 369.2 (df = 4) # Plot ---- qs <- c(.90) Sig.MaxN <-ggplot(Fish.abund.r,aes(x=Temperature,y=Siganus.MaxN))+ geom_point(size=3,position=position_jitter(width=0.05,height=.05))+ geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.5, linetype="dashed")+ geom_smooth(method="lm", se=FALSE, colour="blue", size = 1, alpha = 0.3, linetype="dashed")+ xlab("Temperature")+ ylab(expression(Rabbitfish~Abundance~(MaxN)))+ #ylab(expression(Fish~Abundance~~"(count 125"*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ scale_x_continuous(breaks = seq(16, 22, by = 1))+ geom_hline(yintercept=0,linetype="dashed")+ theme( legend.title=element_blank(), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), #legend.direction="horizontal", legend.position = c(0.2, 0.85), legend.key.size = unit(.7, "cm"), axis.title.x = element_blank(), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1)) Sig.MaxN #----# # Siganus DOV ---- # Mixed-Effects Linear Regression ---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.Siganus.DOV=lme(Siganus.DOV ~ Temperature,random= ~ 1|Reef/Year,data=Fish.abund.r) lme2.Siganus.DOV=lme(Siganus.DOV ~ Temperature,random= ~ 1|Reef,data=Fish.abund.r) lme3.Siganus.DOV=lme(Siganus.DOV ~ Temperature,random= ~ 1+Temperature|Reef,data=Fish.abund.r) lme4.Siganus.DOV=lme(Siganus.DOV ~ Temperature,random= ~ Temperature|Reef,data=Fish.abund.r) AIC(lme1.Siganus.DOV, lme2.Siganus.DOV) # df AIC # lme1.Siganus.DOV 5 265.9784 # lme2.Siganus.DOV 4 267.3212 anova(lme1.Siganus.DOV, lme2.Siganus.DOV) summary(lme1.Siganus.DOV) # no sig differences # Fixed effects: Siganus.DOV ~ Temperature # Value Std.Error DF t-value p-value #(Intercept) -11.199136 15.052750 30 -0.7439927 0.4627 #Temperature 1.093905 0.716004 30 1.5277920 0.1370 # Residuals plot(lme1.Siganus.DOV) # Mixed-effects Quantile Regressions ---- fit.Siganus.DOV <- lqmm(Siganus.DOV ~ Temperature, random = ~ 1, group = Reef, covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Fish.abund.r, control = lqmmControl(method = "df")) fit.Siganus.MaxN$control$LP_tol_ll <- 1e-3 fit.Siganus.MaxN$control$LP_max_iter <- 1000 summary(fit.Siganus.DOV, R = 200, seed = 52) #no sig differences #tau = 0.9 #Fixed effects: # Value Std. Error lower bound upper bound Pr(>|t|) #(Intercept) 16.68647 22.23228 -27.15462 60.5276 0.4538 #Temperature 0.39214 1.31380 -2.19861 2.9829 0.7656 #AIC: 272.5 (df = 4) # Plot ---- qs <- c(0.9) Sig.DOV <-ggplot(Fish.abund.r,aes(x=Temperature,y=Siganus.DOV))+ geom_point(size=3,position=position_jitter(width=0.05,height=.05))+ geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.5, linetype="dashed")+ geom_smooth(method="lm", se=FALSE, colour="blue", size = 1, alpha = 0.3, linetype="dashed")+ xlab(bquote(" "*"Temperature"~degree*C*""))+ #ylab(expression(Fish~Abundance~(MaxN)))+ ylab(expression(Rabbitfish~Abundance~~"(count 125"*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ scale_x_continuous(breaks = seq(16, 22, by = 1))+ geom_hline(yintercept=0,linetype="dashed")+ theme( legend.title=element_blank(), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.position = "none", legend.key.size = unit(.7, "cm"), axis.title.x = element_text(face="bold", colour="153", size=15, margin = unit(c(0.2, 0, 0, 0), "cm")), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1)) Sig.DOV #.----- # Temperate MaxN ----- # Mixed-Effects Linear Regression ---- lme1.Temp.MaxN=lme(Temperate.MaxN ~ Temperature,random= ~ 1|Reef/Year,data=Fish.abund.r) lme2.Temp.MaxN=lme(Temperate.MaxN ~ Temperature,random= ~ 1|Reef,data=Fish.abund.r) lme3.Temp.MaxN=lme(Temperate.MaxN ~ Temperature,random= ~ 1+Temperature|Reef,data=Fish.abund.r) lme4.Temp.MaxN=lme(Temperate.MaxN ~ Temperature,random= ~ Temperature|Reef,data=Fish.abund.r) AIC(lme1.Temp.MaxN, lme2.Temp.MaxN, lme3.Temp.MaxN, lme4.Temp.MaxN) # df AIC # lme1.Temp.MaxN 5 296.2448 # lme2.Temp.MaxN 4 294.5209 - # lme3.Temp.MaxN 6 298.5209 # lme4.Temp.MaxN 6 298.5209 anova(lme2.Temp.MaxN, lme1.Temp.MaxN) summary(lme2.Temp.MaxN) # Fixed effects: Temperate.MaxN ~ Temperature # Value Std.Error DF t-value p-value # (Intercept) -9.217332 23.138115 30 -0.3983614 0.6932 # Temperature 0.984791 1.236099 30 0.7966929 0.4319 #Residuals plot(lme2.Temp.MaxN) # Mixed-effects Quantile Regression ---- fit.Temp.MaxN$control$LP_tol_ll <- 1e-3 fit.Temp.MaxN$control$LP_max_iter <- 1000 fit.Temp.MaxN <- lqmm(Temperate.MaxN ~ Temperature, random = ~ 1, group = Reef, covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Fish.abund.r, control = lqmmControl(method = "df")) summary(fit.Temp.MaxN, R = 500) #tau = 0.9 #Fixed effects: # Value Std. Error lower bound upper bound Pr(>|t|) #(Intercept) -31.6347 105.5361 -243.7174 180.448 0.7656 #Temperature 3.2692 5.8778 -8.5426 15.081 0.5806 # AIC: 343.7 (df = 4) # Plot ---- qs <- c(0.9) Temperate.MaxN <-ggplot(Fish.abund.r,aes(x=Temperature,y=Temperate.MaxN))+ geom_point(size=3,position=position_jitter(width=0.05,height=.05))+ geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.5, linetype="dashed")+ geom_smooth(method="lm", se=FALSE, colour="blue", size = 1, alpha = 0.3 ,linetype="dashed")+ xlab("Temperature")+ ylab(expression(Temperate~Fish~Abundance~(MaxN)))+ #ylab(expression(Fish~Abundance~~"(count 125"*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ scale_x_continuous(breaks = seq(16, 22, by = 1))+ geom_hline(yintercept=0,linetype="dashed")+ ylim(0, 150)+ theme( legend.title=element_blank(), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.position = "null", legend.key.size = unit(.7, "cm"), #axis.title.x = element_blank(), axis.title.x =element_blank(), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1)) Temperate.MaxN #.---- # # Temperate DOV ---- #.---- lme1.Temperate.DOV=lme(Temperate.DOV ~ Temperature,random= ~ 1|Reef/Year,data=Fish.abund.r) lme2.Temperate.DOV=lme(Temperate.DOV ~ Temperature,random= ~ 1|Reef,data=Fish.abund.r) lme3.Temperate.DOV=lme(Temperate.DOV ~ Temperature,random= ~ 1+Temperature|Reef,data=Fish.abund.r) lme4.Temperate.DOV=lme(Temperate.DOV ~ Temperature,random= ~ Temperature|Reef,data=Fish.abund.r) AIC(lme1.Temperate.DOV, lme2.Temperate.DOV, lme3.Temperate.DOV, lme4.Temperate.DOV) # df AIC #lme1.Temperate.DOV 5 232.7760 - #lme2.Temperate.DOV 4 233.5069 - #lme3.Temperate.DOV 6 237.1719 #lme4.Temperate.DOV 6 237.1719 anova(lme1.Temperate.DOV, lme2.Temperate.DOV) # no sig difference summary(lme2.Temperate.DOV) # Fixed effects: Temperate.DOV ~ Temperature # Value Std.Error DF t-value p-value #(Intercept) 5.228492 9.095268 30 0.5748585 0.5697 #Temperature -0.072658 0.484935 30 -0.1498302 0.8819 #Residuals plot(lme4.Temperate.DOV) # Mixed-Effects Quantile Regression ---- fit.Temp.DOV <- lqmm(Temperate.DOV ~ Temperature, random = ~ 1, group = Reef, covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Fish.abund.r, control = lqmmControl(method = "df")) fit.Temp.DOV$control$LP_tol_ll <- 1e-3 fit.Temp.DOV$control$LP_max_iter <- 1000 summary(fit.Temp.DOV, R = 500) #tau = 0.9 #Fixed effects: # Value Std. Error lower bound upper bound Pr(>|t|) #(Intercept) -4.27507 20.92391 -45.53611 36.9860 0.8383 #Temperature 0.91667 1.14795 -1.34704 3.1804 0.4255 #AIC: 258.4 (df = 4) # Plot ---- qs <- c(0.9) Temperate.DOV <-ggplot(Fish.abund.r,aes(x=Temperature,y=Temperate.DOV))+ geom_point(size=3,position=position_jitter(width=0.02,height=.02))+ geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.5, linetype="dashed")+ geom_smooth(method="lm", se=FALSE,colour="blue", size = 1, alpha = 0.3, linetype="dashed")+ xlab(bquote(" "*"Temperature"~degree*C*""))+ #ylab(expression(Temperate~Fish~Abundance~(MaxN)))+ ylab(expression(Temperate~Fish~Abundance~~"(count 125 "*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ scale_x_continuous(breaks = seq(16, 22, by = 1))+ geom_hline(yintercept=0,linetype="dashed")+ ylim(0, 60)+ theme( legend.title=element_blank(), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.position = "Null", legend.key.size = unit(.7, "cm"), #axis.title.x = element_blank(), axis.title.x =element_text(face="bold", colour="153", size=15), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1)) Temperate.DOV #.---- # Join Graphs ---- #.---- # This just shows you on screen what plot will look like grid.arrange(Sig.MaxN,Temperate.MaxN, Sig.DOV, Temperate.DOV, ncol=2) # Use arrangeGrob to pass this to ggsave! Fish.abundances<-arrangeGrob(Sig.MaxN,Temperate.MaxN, Sig.DOV, Temperate.DOV, ncol=2) ggsave(Fish.abundances, file="Fig 1 All Fish Abundances.png",width = 40, height = 30,units = "cm") ####.---- ### Herbivory ##### ---- ####.---- #Data Exploration ----- Herb.sig<-Herbivory%>% select(Temperature, Siganus.MaxN, Siganus.drift.cm.hr, Siganus.browsing.cm.hr, Siganus.kelp.cm.hr, Siganus.kelp.g.hr, Siganus.total.cm.hr) chart.Correlation(Herb.sig, histogram=TRUE, pch=19) Herb.syd<-Herbivory%>% select(Temperature, Sydneyanus.MaxN, Sydneyanus.drift.cm.hr, Sydneyanus.browsing.cm.hr, Sydneyanus.kelp.cm.hr, Sydneyanus.kelp.g.hr, Sydneyanus.total.cm.hr) chart.Correlation(Herb.syd, histogram=TRUE, pch=19) ggplot(data=Herbivory, aes(y=Siganus.kelp.g.hr, x= Siganus.MaxN, color=Year))+ geom_point()+ stat_smooth(method="lm") ggplot(data=Herbivory, aes(y=Siganus.kelp.g.hr, x= Siganus.MaxN, color=Reef))+ geom_point()+ stat_smooth(method="lm") ggplot(data=Herbivory, aes(y=Sydneyanus.kelp.g.hr, x= Sydneyanus.MaxN, color=Year))+ geom_point()+ geom_smooth(method=lm) ggplot(data=Herbivory, aes(y=Sydneyanus.kelp.g.hr, x= Sydneyanus.MaxN, color=Reef))+ geom_point()+ geom_smooth(method=lm) # Kelp Herbivory ##### ---- # Siganus ### ---- # Mixed-Effects Linear Regression ---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim') lme1.kelp.R=lme(Siganus.kelp.cm.hr ~ Siganus.MaxN*Temperature,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.kelp.R=lme(Siganus.kelp.cm.hr ~ Siganus.MaxN*Temperature,random= ~ 1|Reef/Year, control=ctrl,data=Herbivory) lme3.kelp.R=lme(Siganus.kelp.cm.hr ~ Siganus.MaxN*Temperature,random= ~ Siganus.MaxN|Reef, control=ctrl,data=Herbivory) lme4.kelp.R=lme(Siganus.kelp.cm.hr ~ Siganus.MaxN*Temperature,random= ~ Siganus.MaxN|Reef/Year, control=ctrl,data=Herbivory) AIC(lme1.kelp.R,lme2.kelp.R, lme3.kelp.R, lme4.kelp.R) # df AIC #lme1.kelp.R 6 415.3427- #lme2.kelp.R 7 417.3438 #lme3.kelp.R 8 419.1991 #lme4.kelp.R 11 422.8617 anova(lme1.kelp.R, lme4.kelp.R) # no sig difference summary(lme1.kelp.R) # Value Std.Error DF t-value p-value #(Intercept) -2.7868987 10.652389 28 -0.261622 0.7955 #Siganus.MaxN -1.3034676 0.311331 28 -4.186763 0.0003 #Temperature 0.1427152 0.571622 28 0.249667 0.8047 #Siganus.MaxN:Temperature 0.0848805 0.016362 28 5.187712 0.0001 *** # interaction graphic package sjplot plot_model(lme1.kelp.R, type = "int", mdrt.values = "meansd")+ geom_point(mapping = aes(x=Siganus.MaxN, y=Siganus.kelp.cm.hr), inherit.aes = FALSE, data=Herbivory) # Residuals plot p<-plot_model(lme1.kelp.R, type = "diag") plot_grid(p) # Mixed-effects Quantile Regressions ---- fit.kelpg.R.q2$control$LP_tol_ll <- 1e-3 fit.kelpg.R.q2$control$LP_max_iter <- 1000 fit.kelpg.R.q <- lqmm(Siganus.kelp.g.hr ~ Siganus.MaxN*Temperature, random = ~ 1, group = Reef, covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Herbivory, control = lqmmControl(method = "gs")) summary(fit.kelpg.R.q, R = 200, seed = 52) #tau = 0.9 # Value Std. Error lower bound upper bound Pr(>|t|) # (Intercept) -2.7632948 8.0164385 -18.5713632 13.0448 0.7307 # Siganus.MaxN -1.3024392 0.1918059 -1.6806722 -0.9242 1.26e-10 *** # Temperature 0.1425507 0.4651334 -0.7746722 1.0598 0.7596 # Siganus.MaxN:Temperature 0.1010824 0.0099631 0.0814355 0.1207 < 2.2e-16 *** # AIC: [1] 261.8 (df = 6) # Plot kelp.herb.Sig <- plot_model(lme1.kelp.R, type = "int", mdrt.values = "meansd", se=TRUE, colors = c("blue", "orange","firebrick"))+ geom_point(mapping = aes(x=Siganus.MaxN, y=Siganus.kelp.cm.hr), inherit.aes = FALSE, data=Herbivory, size=3)+ #scale_colour_gradient(low = "blue", high = "red")+ #ylab(expression(Kelp~~Consumption~~Rates~~"("*g *""*hr^"-1"*")"))+ ylab(expression(Kelp~~Herbivory~~Rates~~"("*cm^"2"*""*hr^"-1"*")"))+ #xlab(expression(Rabbitfish~Abundance~(MaxN)))+ geom_hline(yintercept=0,linetype="dashed")+ #xlim(0, 150)+ #ylim(0, 1300)+ theme( plot.title = element_blank(), legend.title=element_text(size = 15), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.direction="horizontal", legend.position = c(0.4, 0.95), legend.key.size = unit(1, "cm"), #axis.title.x = element_text(face="bold", colour="153", size=13, # margin=unit(c(0.2, 0, 0, 0), "cm")), axis.title.x = element_blank(), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1))+ annotate("text", size = 5, x = 27, y = 1130, label = "Interaction P < 0.0001", fontface =3) kelp.herb.Sig # Kyphosus ### ---- # Mixed-Effects Linear Regression ---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.kelp.K=lme(Sydneyanus.kelp.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.kelp.K=lme(Sydneyanus.kelp.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ 1|Year/Reef, control=ctrl,data=Herbivory) lme3.kelp.K=lme(Sydneyanus.kelp.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ Sydneyanus.MaxN|Reef, control=ctrl,data=Herbivory) lme4.kelp.K=lme(Sydneyanus.kelp.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ Sydneyanus.MaxN|Year/Reef, control=ctrl,data=Herbivory) AIC(lme1.kelp.K, lme2.kelp.K, lme3.kelp.K, lme4.kelp.K) # df AIC #lme1.kelp.K 6 388.7777 #lme2.kelp.K 7 392.5417 #lme3.kelp.K 8 385.1589 #lme4.kelp.K 11 375.6001 - anova(lme4.kelp.K, lme2.kelp.K) # sig difference! summary(lme4.kelp.K) # Value Std.Error DF t-value p-value #(Intercept) 6.092794 66.39387 20 0.0917674 0.9278 #Sydneyanus.MaxN -14.287698 6.95403 20 -2.0545928 0.0532 #Temperature -0.096073 3.53735 20 -0.0271596 0.9786 #Sydneyanus.MaxN:Temperature 0.942303 0.32047 20 2.9403794 0.0081 # interaction graphic package sjplot plot_model(lme4.kelp.K, type = "int", mdrt.values = "meansd")+ geom_point(mapping = aes(x=Sydneyanus.MaxN, y=Sydneyanus.kelp.cm.hr), inherit.aes = FALSE, data=Herbivory) # Model does not fit points, not a good fit summary(lme1.kelp.K) # Value Std.Error DF t-value p-value #(Intercept) -148.25624 106.21528 28 -1.3958089 0.1737 #Sydneyanus.MaxN -1.90534 6.88549 28 -0.2767183 0.7840 #Temperature 9.09841 5.59620 28 1.6258194 0.1152 #Sydneyanus.MaxN:Temperature 0.11135 0.34845 28 0.3195626 0.7517 # Residuals plot p<-plot_model(lme1.kelp.K, type = "diag") # Seems like non-constant residual variance: heteroscedasticity plot_grid(p) # Quantile Regression Mixed-Effects ---- fit.kelpg.R.q$control$LP_tol_ll <- 1e-3 fit.kelpg.R.q$control$LP_max_iter <- 1000 fit.kelpg.K.q <- lqmm(Sydneyanus.kelp.g.hr ~ Sydneyanus.MaxN*Temperature, random = ~ 1, group = Reef, covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Herbivory, control = lqmmControl(method = "gs")) summary(fit.kelpg.K.q, R = 200, seed = 52) #tau = 0.9 #Fixed effects: # Value Std. Error lower bound upper bound Pr(>|t|) #(Intercept) -10.934609 7.127442 -24.989614 3.1204 0.1266 #Sydneyanus.MaxN 0.062197 0.805334 -1.525886 1.6503 0.9385 #Temperature 0.661703 0.429758 -0.185761 1.5092 0.1252 #Sydneyanus.MaxN:Temperature 0.030796 0.051654 -0.071064 0.1327 0.5517 # AIC: 235.3 (df = 6) # Plot ---- Kyph.herb.kelp <-ggplot(Herbivory,aes(Sydneyanus.MaxN, y=Sydneyanus.kelp.cm.hr))+ geom_point(size=3,position=position_jitter(width=0.05,height=.05), aes(colour=Temperature))+ scale_colour_gradient2(low = "blue", mid= "orange", high = "red", midpoint = 18.5)+ #geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.3)+ #geom_smooth(method="lm", se=FALSE, colour="blue", size = 1, alpha = 0.3)+ ylab(expression(Kelp~~Consumption~~Rates~~"("*g*""*hr^"-1"*")"))+ xlab(expression(K.~Sydneyanus~Abundance~(MaxN)))+ geom_hline(yintercept=0,linetype="dashed")+ #ylab(expression(Fish~Abundance~~"(count 125"*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ #scale_x_continuous(breaks = seq(16, 22, by = 1))+ #xlim(0, 150)+ ylim(0, 1200)+ theme( #legend.title=element_blank(), legend.title=element_text(size = 15), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), #legend.direction="horizontal", legend.position = c(0.8, 0.7), legend.justification='left', legend.key.size = unit(1, "cm"), axis.title.x = element_blank(), axis.title.y = element_blank(), #axis.title.x = element_text(face="bold", colour="153", size=13, # margin=unit(c(0.2, 0, 0, 0), "cm")), #axis.title.y = element_text(face="bold", colour="153", size=13, # margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1)) Kyph.herb.kelp #.---- # Drifting Herbivory ##### ---- #.---- # Siganus ### ---- # Data Exploration ---- ggplot(data=Herbivory, aes(y=Siganus.drift.cm.hr, x= Siganus.MaxN, color=Year))+ geom_point()+ stat_smooth(method="lm") ggplot(data=Herbivory, aes(y=Siganus.drift.cm.hr, x= Siganus.MaxN, color=Reef))+ geom_point()+ stat_smooth(method="lm") # Mixed-Effects Linear Regression ---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.drift.R=lme(Siganus.drift.cm.hr ~ Siganus.MaxN*Temperature,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.drift.R=lme(Siganus.drift.cm.hr ~ Siganus.MaxN*Temperature,random= ~ 1|Year/Reef, control=ctrl,data=Herbivory) lme3.drift.R=lme(Siganus.drift.cm.hr ~ Siganus.MaxN*Temperature,random= ~ Siganus.MaxN|Reef, control=ctrl,data=Herbivory) lme4.drift.R=lme(Siganus.drift.cm.hr ~ Siganus.MaxN*Temperature,random= ~ Siganus.MaxN|Year/Reef, control=ctrl,data=Herbivory) AIC(lme1.drift.R, lme2.drift.R, lme3.drift.R, lme4.drift.R) # df AIC #lme1.drift.R 6 328.8702- #lme2.drift.R 7 330.9581 #lme3.drift.R 8 332.8366 #lme4.drift.R 11 338.9560 anova(lme1.drift.R, lme4.drift.R) # no sig difference summary(lme1.drift.R) # Value Std.Error DF t-value p-value #(Intercept) -57.63233 41.31734 28 -1.3948701 0.1740 #Siganus.MaxN 0.47549 1.21058 28 0.3927809 0.6975 #Temperature 3.67787 2.21312 28 1.6618435 0.1077 #Siganus.MaxN:Temperature -0.02113 0.06351 28 -0.3327472 0.7418 # Quantile Regression Mixed-Effects ---- fit.drift.R.q2$control$LP_tol_ll <- 1e-3 fit.drift.R.q2$control$LP_max_iter <- 1000 fit.drift.R.q2 <- lqmm(Siganus.drift.cm.hr ~ Siganus.MaxN*Temperature, random = ~ Siganus.MaxN, group = Reef, #random intercepts and slopes covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Herbivory, control = lqmmControl(method = "df")) summary(fit.drift.R.q2, R = 200, seed = 52) #tau = 0.9 #Fixed effects: # Value Std. Error lower bound upper bound Pr(>|t|) #(Intercept) -3.1534e+01 4.4882e+01 -1.2004e+02 56.9713 0.4831 #Siganus.MaxN -1.6530e-01 7.3889e+00 -1.4736e+01 14.4053 0.9822 #Temperature 4.3601e+00 3.2467e+00 -2.0423e+00 10.7624 0.1808 #Siganus.MaxN:Temperature -9.5895e-03 4.0435e-01 -8.0694e-01 0.7878 0.9811 # AIC: 349.5 (df = 8) # Plot ---- Sig.herb.drift <-ggplot(Herbivory,aes(x=Siganus.MaxN,y=Siganus.drift.cm.hr))+ geom_point(size=3,position=position_jitter(width=0.05,height=.05), aes(colour=Temperature))+ scale_colour_gradient2(low = "blue", mid= "orange", high = "red", midpoint = 18.5)+ #geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.3)+ #geom_smooth(method="lm", se=FALSE, colour="blue", size = 1, alpha = 0.3)+ ylab(expression(Drifting~~Herbivory~~Rates~~"("*cm^"2"*""*hr^"-1"*")"))+ xlab(expression(Rabbitfish~Abundance~(MaxN)))+ geom_hline(yintercept=0,linetype="dashed")+ #ylab(expression(Fish~Abundance~~"(count 125"*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ #scale_x_continuous(breaks = seq(16, 22, by = 1))+ ylim(0, 100)+ theme( #legend.title=element_blank(), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 11), #legend.direction="horizontal", legend.position = "null", #legend.position = c(0.8, 0.85), legend.justification='left', legend.key.size = unit(.5, "cm"), axis.title.x = element_text(face="bold", colour="153", size=15, margin=unit(c(0.2, 0, 0, 0), "cm")), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1)) Sig.herb.drift # Kyphosus ### ---- # Data Exploration ---- ggplot(data=Herbivory, aes(y=Sydneyanus.drift.cm.hr, x= Sydneyanus.MaxN, color=Reef))+ geom_point()+ stat_smooth(method="lm") ggplot(data=Herbivory, aes(y=Sydneyanus.drift.cm.hr, x= Sydneyanus.MaxN, color=Year))+ geom_point()+ stat_smooth(method="lm") # Mixed-Effects Linear Regression ---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.drift.K=lme(Sydneyanus.drift.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.drift.K=lme(Sydneyanus.drift.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ 1|Year/Reef, control=ctrl,data=Herbivory) lme3.drift.K=lme(Sydneyanus.drift.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ Sydneyanus.MaxN|Reef, control=ctrl,data=Herbivory) lme4.drift.K=lme(Sydneyanus.drift.cm.hr ~ Sydneyanus.MaxN*Temperature,random= ~ Sydneyanus.MaxN|Year/Reef, control=ctrl,data=Herbivory) AIC(lme1.drift.K, lme2.drift.K, lme3.drift.K, lme4.drift.K) # df AIC #lme1.drift.K 6 219.5555 #lme2.drift.K 7 220.4116 #lme3.drift.K 8 212.1436 - #lme4.drift.K 11 215.8004 anova(lme1.drift.K, lme3.drift.K) # sig dif summary(lme3.drift.K) # no sig difference # Value Std.Error DF t-value p-value #(Intercept) -3.0118300 5.973175 28 -0.5042260 0.6180 #Sydneyanus.MaxN -0.8123009 0.550326 28 -1.4760354 0.1511 #Temperature 0.1722853 0.319694 28 0.5389072 0.5942 #Sydneyanus.MaxN:Temperature 0.0541979 0.027070 28 2.0021271 0.0550 # Quantile Regression Mixed-Effects ---- fit.drift.K.q$control$LP_tol_ll <- 1e-3 fit.drift.K.q$control$LP_max_iter <- 1000 fit.drift.K.q2 <- lqmm(Sydneyanus.drift.cm.hr ~ Sydneyanus.MaxN*Temperature, random = ~ Sydneyanus.MaxN, group = Reef, #random intercepts and slopes covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Herbivory, control = lqmmControl(method = "gs")) summary(fit.drift.K.q2, R = 500, seed = 52) #tau = 0.9 # Value Std. Error lower bound upper bound Pr(>|t|) #(Intercept) -2.116500 4.788889 -11.559981 7.3270 0.6590 #Sydneyanus.MaxN -0.793787 0.516498 -1.812298 0.2247 0.1259 #Temperature 0.124323 0.304029 -0.475209 0.7239 0.6830 #Sydneyanus.MaxN:Temperature 0.052649 0.033089 -0.012601 0.1179 0.1132 #AIC: 224.1 (df = 6) # Plot ---- Kyph.herb.drift <-ggplot(Herbivory,aes(x=Sydneyanus.MaxN,y=Sydneyanus.drift.cm.hr))+ geom_point(size=3,position=position_jitter(width=0.05,height=.05), aes(colour=Temperature))+ scale_colour_gradient2(low = "blue", mid= "orange", high = "red", midpoint = 18.5)+ #geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.3)+ #geom_smooth(method="lm", se=FALSE, colour="blue", size = 1, alpha = 0.3)+ ylab(expression(Drifting~~Herbivory~~Rates~~"("*cm^"2"*""*hr^"-1"*")"))+ xlab(expression(K.~sydneyanus~Abundance~(MaxN)))+ geom_hline(yintercept=0,linetype="dashed")+ #ylab(expression(Fish~Abundance~~"(count 125"*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ #scale_x_continuous(breaks = seq(16, 22, by = 1))+ ylim(0, 100)+ theme( #legend.title=element_blank(), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 10), legend.direction="horizontal", legend.position = "null", #legend.position = c(0.05, 0.9), legend.justification='left', legend.key.size = unit(.5, "cm"), axis.title.x = element_text(face="bold", colour="153", size=15, margin=unit(c(0.2, 0, 0, 0), "cm")), #axis.title.y = element_text(face="bold", colour="153", size=13, # margin=unit(c(0, 0.2, 0, 0), "cm")), axis.title.y = element_blank(), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1)) Kyph.herb.drift ## Join Plots #####---- # This just shows you on screen what plot will look like grid.arrange(kelp.herb.Sig, Kyph.herb.kelp, Sig.herb.drift, Kyph.herb.drift, ncol=2) # Use arrangeGrob to pass this to ggsave! Fish.herbivory<-arrangeGrob(kelp.herb.Sig, Kyph.herb.kelp, Sig.herb.drift, Kyph.herb.drift, ncol=2) ggsave(Fish.herbivory, file="FIG 3. All FISH HERBIVORY.png",width = 40, height = 30,units = "cm") #.---- # Excretion Rates ##### ---- #.---- # Siganus ### ---- # Data Exploration ---- Excr.sig<-Herbivory%>% select(Temperature, Siganus.MaxN, Siganus.total.cm.hr, Rabbit.faeces.kg.hr) chart.Correlation(Excr.sig, histogram=TRUE, pch=19) ggplot(data=Herbivory, aes(y=Rabbit.faeces.kg.hr, x= Siganus.MaxN, color=Reef))+ geom_point()+ stat_smooth(method="lm") ggplot(data=Herbivory, aes(y=Rabbit.faeces.kg.hr, x= Siganus.MaxN, color=Year))+ geom_point()+ stat_smooth(method="lm") ggplot(data=Herbivory, aes(y=Rabbit.faeces.kg.hr, x= Siganus.MaxN, color=Reef))+ geom_point()+ stat_smooth(method="lm") # Transformation # Defecation Herbivory$Rabbit.faeces.kg.hr.log <- transform(Herbivory$Rabbit.faeces.kg.hr, method = "log+1") plot(Herbivory$Rabbit.faeces.kg.hr.log) Herbivory$Rabbit.faeces.kg.hr.log = as.numeric(Herbivory$Rabbit.faeces.kg.hr.log) # Abundance Herbivory$Siganus.MaxN.log <- transform(Herbivory$Siganus.MaxN, method = "log+1") plot(Herbivory$Siganus.MaxN.log) Herbivory$Siganus.MaxN.log = as.numeric(Herbivory$Siganus.MaxN.log) # Herbivory Herbivory$Siganus.total.cm.hr.log <- transform(Herbivory$Siganus.total.cm.hr, method = "log+1") plot(Herbivory$Siganus.total.cm.hr.log) Herbivory$Siganus.total.cm.hr.log = as.numeric(Herbivory$Siganus.total.cm.hr.log) # Mixed-Effects Linear Regression---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.faeces.R.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.MaxN.log*Temperature,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.faeces.R.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.MaxN.log*Temperature,random= ~ 1|Year/Reef, control=ctrl,data=Herbivory) lme3.faeces.R.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.MaxN.log*Temperature,random= ~ Siganus.MaxN|Reef, control=ctrl,data=Herbivory) lme4.faeces.R.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.MaxN.log*Temperature,random= ~ Siganus.MaxN|Year/Reef, control=ctrl,data=Herbivory) AIC(lme1.faeces.R.l, lme2.faeces.R.l, lme3.faeces.R.l, lme4.faeces.R.l) # df AIC #lme1.faeces.R.l 6 65.82223 - #lme2.faeces.R.l 7 67.82591 #lme3.faeces.R.l 8 69.82294 #lme4.faeces.R.l 11 75.82486 anova(lme1.faeces.R.l, lme3.faeces.R.l) # no sig difference summary(lme1.faeces.R.l) # Value Std.Error DF t-value p-value #(Intercept) -0.0428889 0.9740133 28 -0.044033 0.9652 #Siganus.MaxN.log -1.2243932 0.3822240 28 -3.203339 0.0034 #Temperature 0.0005672 0.0542511 28 0.010456 0.9917 #Siganus.MaxN.log:Temperature 0.0802298 0.0208175 28 3.853957 0.0006 plot_model(lme1.faeces.R.l, type = "int", mdrt.values = "meansd")+ geom_point(mapping = aes(x=Siganus.MaxN.log, y=Rabbit.faeces.kg.hr.log), inherit.aes = FALSE, data=Herbivory) # Residuals plot p<-plot_model(lme1.faeces.R.l, type = "diag") # Seems like non-constant residual variance: heteroscedasticity plot_grid(p) # Quantile Regression Mixed-Effects ---- fit.faeces.R.q$LP_tol_ll <- 1e-3 fit.faeces.R.q$LP_max_iter <- 1000 fit.faeces.R.q <- lqmm(Rabbit.faeces.kg.hr.log ~ Siganus.MaxN.log*Temperature, random = ~ Siganus.MaxN.log, group = Reef, #random intercepts covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Herbivory, control = lqmmControl(method = "gs")) fit.boot.faeces.R.q <- boot(fit.faeces.R.q) summary(fit.boot.faeces.R.q, R = 500, seed = 52) #tau = 0.9 # Value Std. Error lower bound upper bound Pr(>|t|) # (Intercept) -0.025093 0.822287 -1.656688 1.6065 0.9757 # Siganus.MaxN.log -1.230551 0.299861 -1.825541 -0.6356 8.360e-05 *** # Temperature 0.015661 0.051261 -0.086052 0.1174 0.7606 # Siganus.MaxN.log:Temperature 0.088579 0.015164 0.058491 0.1187 6.613e-08 *** # Plot ---- Sig.faeces.int <- plot_model(lme1.faeces.R.l, type = "int", mdrt.values = "meansd", colors = c("blue", "orange","firebrick"))+ geom_point(mapping = aes(x=Siganus.MaxN.log, y=Rabbit.faeces.kg.hr.log), size=3, inherit.aes = FALSE, data=Herbivory)+ ylab(expression(Excretion~~Rates~~"log("*kg *" "*hr^"-1"*")"))+ xlab(expression(Rabbitfish~Abundance~log(MaxN)))+ geom_hline(yintercept=0,linetype="dashed")+ #xlim(0, 85)+ #ylim(0, 3)+ theme( plot.title = element_blank(), legend.title=element_text(size = 15), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.direction="horizontal", legend.position = c(0.3, 0.95), legend.key.size = unit(1, "cm"), #axis.title.x = element_blank(), axis.title.x = element_text(face="bold", colour="153", size=15, margin=unit(c(0.2, 0, 0, 0), "cm")), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1))+ annotate("text", size = 5, x = 0.5, y = 2, label = "P > 0.001", fontface =3) Sig.faeces.int # Herbivory rates ---- # Mixed-Effects Linear Regression---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.faeces.R.h.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.total.cm.hr.log,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.faeces.R.h.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.total.cm.hr.log,random= ~ 1|Year/Reef, control=ctrl,data=Herbivory) lme3.faeces.R.h.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.total.cm.hr.log,random= ~ Siganus.total.cm.hr.log|Reef, control=ctrl,data=Herbivory) lme4.faeces.R.h.l=lme(Rabbit.faeces.kg.hr.log ~ Siganus.total.cm.hr.log,random= ~ Siganus.total.cm.hr.log|Year/Reef, control=ctrl,data=Herbivory) AIC(lme1.faeces.R.h.l, lme2.faeces.R.h.l, lme3.faeces.R.h.l, lme4.faeces.R.h.l) # df AIC #lme1.faeces.R.h.l 4 6.471884 #lme2.faeces.R.h.l 5 5.906192 #lme3.faeces.R.h.l 6 9.854969 #lme4.faeces.R.h.l 9 -14.978218 anova(lme1.faeces.R.l, lme4.faeces.R.l) # no sig difference summary(lme2.faeces.R.h.l) # Value Std.Error DF t-value p-value #(Intercept) 0.4078919 0.3422015 30 1.191964 0.2426 #Siganus.total.cm.hr 0.0083138 0.0010237 30 8.121328 0.0001 plot_model(lme2.faeces.R.h.l, type = "pred", terms = c("Siganus.total.cm.hr.log"))+ geom_point(mapping = aes(x=Siganus.total.cm.hr.log, y=Rabbit.faeces.kg.hr.log), inherit.aes = FALSE, data=Herbivory) # Residuals plot p<-plot_model(lme2.faeces.R.h, type = "diag") plot_grid(p) # Plot ---- Sig.faeces.herb.l <-plot_model(lme2.faeces.R.h.l, type = "pred", terms = c("Siganus.total.cm.hr.log"))+ geom_point(mapping = aes(x=Siganus.total.cm.hr.log, y=Rabbit.faeces.kg.hr.log, colour=Temperature), size=3,position=position_jitter(width=0.05,height=.05), inherit.aes = FALSE, data=Herbivory)+ scale_colour_gradient2(low = "blue", mid="orange",high = "red", midpoint = 18.5)+ ylab(expression(Excretion~~Rates~~"log("*kg *" "*hr^"-1"*")"))+ xlab(expression(Total~~Herbivory~~Rates~~"log("*cm^"2"*""*hr^"-1"*")"))+ geom_hline(yintercept=0,linetype="dashed")+ #xlim(0, 85)+ ylim(0, 2.5)+ theme( plot.title = element_blank(), legend.title=element_text(size = 15), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.direction="horizontal", #legend.position = c(0.25, 0.90), legend.position = "null", legend.key.size = unit(1, "cm"), #axis.title.x = element_blank(), axis.title.x = element_text(face="bold", colour="153", size=15, margin=unit(c(0.2, 0, 0, 0), "cm")), axis.title.y = element_text(face="bold", colour="153", size=15, margin=unit(c(0, 0.2, 0, 0), "cm")), #strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1))+ annotate("text", size = 5, x = 1, y = 2.1, label = "P < 0.0001", fontface =3) Sig.faeces.herb.l # Kyphosus ### ---- # Data Exploration ---- Excr.Kyph<-Herbivory%>% select(Temperature, Sydneyanus.MaxN, Sydneyanus.total.cm.hr, Sydneyanus.faeces.kg.hr ) chart.Correlation(Excr.Kyph, histogram=TRUE, pch=19) ggplot(data=Herbivory, aes(y=Sydneyanus.faeces.kg.hr, x= Sydneyanus.MaxN, color=Year))+ geom_point()+ stat_smooth(method="lm") ggplot(data=Herbivory, aes(y=Sydneyanus.faeces.kg.hr, x= Sydneyanus.MaxN, color=Reef))+ geom_point()+ stat_smooth(method="lm") # Transformation # Defecation Herbivory$Sydneyanus.faeces.kg.hr.log <- transform(Herbivory$Sydneyanus.faeces.kg.hr, method = "log+1") plot(Herbivory$Sydneyanus.faeces.kg.hr.log) Herbivory$Sydneyanus.faeces.kg.hr.log = as.numeric(Herbivory$Sydneyanus.faeces.kg.hr.log) # Abundance Herbivory$Sydneyanus.MaxN.log <- transform(Herbivory$Sydneyanus.MaxN, method = "log+1") plot(Herbivory$Sydneyanus.MaxN.log) Herbivory$Sydneyanus.MaxN.log = as.numeric(Herbivory$Sydneyanus.MaxN.log) #Herbivory Herbivory$Sydneyanus.total.cm.hr.log <- transform(Herbivory$Sydneyanus.total.cm.hr, method = "log+1") plot(Herbivory$Sydneyanus.total.cm.hr.log) Herbivory$Sydneyanus.total.cm.hr.log = as.numeric(Herbivory$Sydneyanus.total.cm.hr.log) # Mixed-Effects Linear Regression ---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.faeces.K=lme(Sydneyanus.faeces.kg.hr ~ Sydneyanus.MaxN*Temperature,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.faeces.K=lme(Sydneyanus.faeces.kg.hr ~ Sydneyanus.MaxN*Temperature,random= ~ 1|Year/Reef, control=ctrl,data=Herbivory) lme3.faeces.K=lme(Sydneyanus.faeces.kg.hr ~ Sydneyanus.MaxN*Temperature,random= ~ Sydneyanus.MaxN|Reef, control=ctrl,data=Herbivory) lme4.faeces.K=lme(Sydneyanus.faeces.kg.hr ~ Sydneyanus.MaxN*Temperature,random= ~ Sydneyanus.MaxN|Year/Reef, control=ctrl,data=Herbivory) AIC(lme1.faeces.K, lme2.faeces.K, lme3.faeces.K, lme4.faeces.K) # df AIC #lme1.faeces.K 6 91.62214 #lme2.faeces.K 7 93.26152 #lme3.faeces.K 8 83.18236 #lme4.faeces.K 11 80.57276 - anova(lme1.faeces.K, lme4.faeces.K) # Sig difference summary(lme4.faeces.K) # Value Std.Error DF t-value p-value #(Intercept) -0.7002369 0.5358436 20 -1.3067935 0.2061 #Sydneyanus.MaxN -0.0449324 0.0615229 20 -0.7303365 0.4737 #Temperature 0.0391583 0.0285113 20 1.3734322 0.1848 #Sydneyanus.MaxN:Temperature 0.0041932 0.0025333 20 1.6552597 0.1135 # Residuals plot p<-plot_model(lme4.faeces.K, type = "diag") plot_grid(p) # Quantile Regression Mixed-Effects ---- fit.drift.K.q$control$LP_tol_ll <- 1e-3 fit.drift.K.q$control$LP_max_iter <- 1000 fit.faeces.K.q <- lqmm(Sydneyanus.faeces.kg.hr ~ Sydneyanus.MaxN*Temperature, random = ~ 1, group = Reef, #random intercepts covariance = "pdSymm", tau = c(0.90), nK = 7, type = "normal", data = Herbivory, control = lqmmControl(method = "gs")) summary(fit.faeces.K.q, R = 200, seed = 52) #tau = 0.9 # Value Std. Error lower bound upper bound Pr(>|t|) # (Intercept) -1.2517179 0.4401308 -2.1196368 -0.3838 0.004921 ** # Sydneyanus.MaxN 0.0016497 0.0788522 -0.1538434 0.1571 0.983329 # Temperature 0.0762082 0.0267478 0.0234628 0.1290 0.004845 ** # Sydneyanus.MaxN:Temperature 0.0024137 0.0061342 -0.0096827 0.0145 0.694385 #AIC: 97.56 (df = 6) # Plot ---- Kyph.faeces <-ggplot(Herbivory,aes(x=Sydneyanus.MaxN.log,y=Sydneyanus.faeces.kg.hr.log))+ geom_point(size=3, aes(colour=Temperature))+ scale_colour_gradient2(low = "blue", mid="orange", high = "red", midpoint = 18.5)+ #geom_quantile(quantiles=qs, colour="red", size = 1, alpha = 0.3)+ #geom_smooth(method="lm", se=FALSE,colour="blue", size = 1, alpha = 0.3)+ #ylab(expression(k.~Sydneyanus~~Excretion~~Rates~~"("*kg *" "*hr^"-1"*")"))+ xlab(expression(K.~Sydneyanus~Abundance~~log(MaxN)))+ geom_hline(yintercept=0,linetype="dashed")+ #ylab(expression(Fish~Abundance~~"(count 125"*m^"-2"*")"))+ #scale_fill_manual(values=c("#66FF00", "#33FFCC", "#66FF00", "#33FFCC", "#66FF00"))+ #scale_x_continuous(breaks = seq(0, 150, by = 25))+ ylim(0, 2.5) + theme( #legend.title=element_blank(), legend.title = element_text(size = 15), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.direction="vertical", legend.position = c(0.8, 0.7), #legend.position = "none", legend.justification='left', legend.key.size = unit(1, "cm"), axis.title.x = element_text(face="bold", colour="153", size=15, margin=unit(c(0.2, 0, 0, 0), "cm")), axis.title.y = element_blank(), #axis.title.y = element_text(face="bold", colour="153", size=15, # margin=unit(c(0, 0.2, 0, 0), "cm")), strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1))+ annotate("text", size = 5, x = 1, y = 2, label = "Temperature P = 0.005", fontface =3) Kyph.faeces # Mixed-Effects Linear Regression BITE RATES WITH TRANSFORMATION---- ctrl <- lmeControl(maxIter= 1000, msMaxIter= 1000, msMaxEval=1000, tolerance= 1e-3, opt='optim'); lme1.faeces.k.h.l=lme(Sydneyanus.faeces.kg.hr.log ~ Sydneyanus.total.cm.hr.log,random= ~ 1|Reef, control=ctrl,data=Herbivory) lme2.faeces.k.h.l=lme(Sydneyanus.faeces.kg.hr.log ~ Sydneyanus.total.cm.hr.log,random= ~ 1|Year/Reef, control=ctrl,data=Herbivory) lme3.faeces.k.h.l=lme(Sydneyanus.faeces.kg.hr.log ~ Sydneyanus.total.cm.hr.log,random= ~ Sydneyanus.MaxN|Reef, control=ctrl,data=Herbivory) lme4.faeces.k.h.l=lme(Sydneyanus.faeces.kg.hr.log ~ Sydneyanus.total.cm.hr.log,random= ~ Sydneyanus.MaxN|Year/Reef, control=ctrl,data=Herbivory) AIC(lme1.faeces.k.h.l, lme2.faeces.k.h.l, lme3.faeces.k.h.l, lme4.faeces.k.h.l) # df AIC #lme1.faeces.k.h 4 7.523099 - #lme2.faeces.k.h 5 9.810609 #lme3.faeces.k.h 6 11.489396 #lme4.faeces.k.h 9 18.059047 anova(lme1.faeces.k.h.l, lme4.faeces.k.h.l) # no sig difference summary(lme1.faeces.k.h.l) # Value Std.Error DF t-value p-value #(Intercept) 0.00890964 0.04822945 30 0.184734 0.8547 #Sydneyanus.total.cm.hr.log 0.12706031 0.02082304 30 6.101908 0.0001 plot_model(lme1.faeces.k.h.l, type = "pred", terms = c("Sydneyanus.total.cm.hr.log"))+ geom_point(mapping = aes(x=Sydneyanus.total.cm.hr.log, y=Sydneyanus.faeces.kg.hr.log), inherit.aes = FALSE, data=Herbivory) # Residuals p<-plot_model(lme1.faeces.k.h.l, type = "diag") plot_grid(p) # Plot --- Syd.faeces.herb.l <-plot_model(lme1.faeces.k.h.l, type = "pred", terms = c("Sydneyanus.total.cm.hr.log"))+ geom_point(mapping = aes(x=Sydneyanus.total.cm.hr.log, y=Sydneyanus.faeces.kg.hr.log, colour=Temperature), size=3,position=position_jitter(width=0.01,height=.01), inherit.aes = FALSE, data=Herbivory)+ scale_colour_gradient2(low = "blue", mid = "orange", high = "red", midpoint = 18.5)+ ylab(expression(Excretion~~Rates~~"log("*kg *" "*hr^"-1"*")"))+ xlab(expression(Total~~Herbivory~~Rates~~"log("*cm^"2"*" "*hr^"-1"*")"))+ geom_hline(yintercept=0,linetype="dashed")+ #xlim(0, 85)+ ylim(0, 2.3)+ theme( plot.title = element_blank(), legend.title=element_text(size = 15), axis.text.x = element_text(colour="153", size=15), axis.text.y = element_text(colour="153", size=15), panel.background = element_rect(fill="white"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.text = element_text(size = 15), legend.direction="horizontal", #legend.position = c(0.25, 0.90), legend.position = "null", legend.key.size = unit(1, "cm"), axis.title.y = element_blank(), axis.title.x = element_text(face="bold", colour="153", size=15, margin=unit(c(0.2, 0, 0, 0), "cm")), #axis.title.y = element_text(face="bold", colour="153", size=15, # margin=unit(c(0, 0.2, 0, 0), "cm")), #strip.text.x = element_text(size=15, face="bold"), strip.background = element_rect(fill="#FFFFFF"))+ theme(axis.line.x = element_line(size = 1, colour = "black"), axis.line.y = element_line(color="black", size = 1))+ annotate("text", size = 5, x = 0.5, y = 1.3, label = "P < 0.001", fontface =3) Syd.faeces.herb.l #.---- ## Join Plots #####---- # This just shows you on screen what plot will look like grid.arrange(Sig.faeces.int, Kyph.faeces, Sig.faeces.herb.l, Syd.faeces.herb.l, ncol=2) # Use arrangeGrob to pass this to ggsave! Fish.excretion <-arrangeGrob(Sig.faeces.int, Kyph.faeces, Sig.faeces.herb.l, Syd.faeces.herb.l, ncol=2) ggsave(Fish.excretion, file="Fig. 4. Defecation Rates.png",width = 40, height = 30,units = "cm")