---
title: "Continuous bubble streams for controlling marine biofouling on static artificial structures"
author: "Grant A. Hopkins, Fletcher Gilbertson, Oliver Floerl, Paula Casanovas, Matt Pine, Patrick Cahill"
date: "2 March 2021"
output: 
  word_document:
    toc: true 
    toc_depth: 5
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(knitr)
```

```{r warning=FALSE, results=TRUE, echo=TRUE, message=FALSE}
rm(list=ls())

library(tidyverse)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(reshape2)
library(nortest)
library(lmerTest)
library(emmeans)
library(MASS)
library(sjmisc)
library(sjPlot)
library(flextable)

#Lab data
bubble_lab <- read.csv("S6 Lab data.csv")
summary(bubble_lab)

bubble_lab_experiment <- bubble_lab %>% 
  group_by(Mechanism, Organism,  SurfaceType, Settlementperiodhrs,Flowrate) %>%
  summarise(repetitions = length(COUNTresponse)) %>%
  data.frame()

bubble_lab$Settlementperiodhrs <- as.factor(bubble_lab$Settlementperiodhrs)
bubble_lab$Flowrate <- factor(bubble_lab$Flowrate, levels = c("Nil",  "Low",  "Med", "High"))

bubble_lab$Settlementperiodhrs  <- sub(3, "3 hours", bubble_lab$Settlementperiodhrs)
bubble_lab$Settlementperiodhrs  <- sub(120, "120 hours", bubble_lab$Settlementperiodhrs)
bubble_lab$Settlementperiodhrs  <- factor(bubble_lab$Settlementperiodhrs, levels = c("3 hours", "120 hours"))

#Field data

##Vertical surfaces
field_vertical <- read.csv("S7a Field data (vertical surfaces).csv")

##Horizontal and angled surfaces 
field_horizontal_angled <- read.csv("S7b Field data (horizontal angled surfaces).csv")
```


## Methods

### Statistical analyses

#### Laboratory trials

We performed generalized linear models (Dobson 1990) to test how flow rate, settlement periond and surface type affect the efficiency of the bubble treatment. For the scouring and disruption trials, we modeled individual counts of Ciona and Crassostrea using a negative binomial generalized model, because the ratio between the residual deviance and the degrees of freedom showed that Poisson models were overdispersed. For each trial, we did a different model for each species. When appropiate, we model a subset of the data to undesrtand the differences better. To test if there were differences between the different factor combinations of the models, we calculated the predicted marginal means from the models and compared them with a *z*-test. 


```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
##### Scouring trial
## scouring trial plot

bubble_lab_scouring <- bubble_lab[bubble_lab$Mechanism=="Scouring",]
bubble_lab_scouring$SurfaceType <- ifelse(bubble_lab_scouring$SurfaceType=="FR", "FR-IS1100",bubble_lab_scouring$SurfaceType)


ggplot_bubble_scouring <- ggplot(aes(y = COUNTresponse, x = Flowrate, col=SurfaceType, fill=SurfaceType), data = bubble_lab_scouring) +
  geom_boxplot(size=1, alpha=0.8) +
  scale_colour_manual("", values = c("#999999","#000000ff")) +
  scale_fill_manual("", values = c( "#999999","#000000ff")) +
  labs(y="Number of individuals", x="Flow rate")  +
  theme_bw(base_size = 18) +
  theme(legend.key.size = unit(1, "cm"), axis.ticks = element_blank(),strip.text.y = element_text(face = "italic"), legend.position = "bottom") +
  facet_grid(cols = vars(Settlementperiodhrs), rows = vars(Organism), scales = "free_y", switch = 'x') 

##analyses 

#Ciona
#Poisson model - overdispersed
ciona_scouring_glm_poisson <- glm(COUNTresponse ~   Flowrate * Settlementperiodhrs, data=bubble_lab_scouring[bubble_lab_scouring$Organism=="Ciona",], family="poisson")

summary(ciona_scouring_glm_poisson)

#Quasipoisoon model - better fit, still not that good
ciona_scouring_glm_quasipoisson <- glm(COUNTresponse ~   Flowrate * Settlementperiodhrs, data=bubble_lab_scouring[bubble_lab_scouring$Organism=="Ciona",], family="quasipoisson")

summary(ciona_scouring_glm_quasipoisson)

#Negative binomial model - best we could do?
ciona_scouring_glm_NB <- glm.nb(COUNTresponse ~   Flowrate * Settlementperiodhrs, data=bubble_lab_scouring[bubble_lab_scouring$Organism=="Ciona",], link="log")

summary(ciona_scouring_glm_NB)

#Plot the interaction
#emmip(ciona_scouring_glm_NB, Flowrate ~ Settlementperiodhrs)

#plot_model(ciona_scouring_glm_NB, type = "pred", terms = c("Flowrate", "Settlementperiodhrs"))

#Get the pairwise comparison from the marginal means
ciona_scouring_emm <- emmeans(ciona_scouring_glm_NB, ~Flowrate * Settlementperiodhrs)
ciona_scouring_emm_pairs <- data.frame(pairs(ciona_scouring_emm))

ciona_scouring_emm_pairs <- ciona_scouring_emm_pairs[order(ciona_scouring_emm_pairs$p.value),]
ciona_scouring_emm_pairs$p.value <- round(ciona_scouring_emm_pairs$p.value,5)

ciona_scouring_emm_pairs

#Crassostrea
Crassostrea_scouring_glm_poisson <- glm(COUNTresponse ~   Flowrate * Settlementperiodhrs * SurfaceType, data=bubble_lab_scouring[bubble_lab_scouring$Organism=="Crassostrea",], family="poisson")

summary(Crassostrea_scouring_glm_poisson)

Crassostrea_scouring_glm_quasipoisson <- glm(COUNTresponse ~   Flowrate * Settlementperiodhrs * SurfaceType, data=bubble_lab_scouring[bubble_lab_scouring$Organism=="Crassostrea",], family="quasipoisson")

summary(Crassostrea_scouring_glm_quasipoisson)

Crassostrea_scouring_glm_NB <- glm.nb(COUNTresponse ~   Flowrate * Settlementperiodhrs * SurfaceType, data=bubble_lab_scouring[bubble_lab_scouring$Organism=="Crassostrea",], link="log")

summary(Crassostrea_scouring_glm_NB)

#plot the three way interaction
#emmip(Crassostrea_scouring_glm_NB, Flowrate ~ Settlementperiodhrs | SurfaceType)
#emmip(Crassostrea_scouring_glm_NB, Settlementperiodhrs ~ Flowrate | SurfaceType)

#plot_model(Crassostrea_scouring_glm_NB, type = "pred", terms = c("Flowrate", "Settlementperiodhrs", "SurfaceType"))

#pairwise comparison
Crassostrea_scouring_emm <- emmeans(Crassostrea_scouring_glm_NB, ~Flowrate * Settlementperiodhrs * SurfaceType)
Crassostrea_scouring_emm_pairs <- data.frame(pairs(Crassostrea_scouring_emm))

Crassostrea_scouring_emm_pairs <- Crassostrea_scouring_emm_pairs[order(Crassostrea_scouring_emm_pairs$p.value),]
Crassostrea_scouring_emm_pairs$p.value <- round(Crassostrea_scouring_emm_pairs$p.value,5)


Crassostrea_scouring_glm_NB_simple <- emmeans::contrast(Crassostrea_scouring_emm,
                     method = "revpairwise",
                     simple = "each",
                     combine = TRUE,
                     adjust = "none") %>%
  summary(infer = TRUE)

# #model only for FR surface type and settlement period 120 hrs
Crassostrea_scouring_glm_NB_120_FR <- glm.nb(COUNTresponse ~   Flowrate, data=bubble_lab_scouring[bubble_lab_scouring$Organism=="Crassostrea" & bubble_lab_scouring$SurfaceType %in% c("FR-IS1100") & bubble_lab_scouring$Settlementperiodhrs=="120 hours",], link="log")
```


```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
##### Disruption trial

### disruption trial plot

bubble_lab_disruption <- bubble_lab[bubble_lab$Mechanism=="Disruption",]

bubble_lab_disruption$Flowrate <- factor(bubble_lab_disruption$Flowrate, levels=c("Nil","Med"))

bubble_lab_disruption$SurfaceType <- ifelse(bubble_lab_disruption$SurfaceType=="FR", "FR-IS1100",bubble_lab_disruption$SurfaceType)

ggplot_bubble_disruption <- ggplot(aes(y = COUNTresponse, x = Flowrate, col=SurfaceType, fill=SurfaceType), data = bubble_lab_disruption) +
  geom_boxplot(size=1.2, alpha=0.8) +
  scale_colour_manual("", values = c("#999999","#000000ff")) +
  scale_fill_manual("", values = c( "#999999","#000000ff")) +
  labs(y="Number of individuals", x="Flow rate")  +
  theme_bw(base_size = 18) +
  theme(legend.key.size = unit(1, "cm"), axis.ticks = element_blank(),strip.text.y = element_text(face = "italic"), legend.position = "bottom") +
  facet_grid(rows = vars(Organism), scales = "free_y", switch = 'x') 

#Ciona
ciona_disruption_glm_poisson <- glm(COUNTresponse ~ 0 +  Flowrate * SurfaceType, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Ciona",], family="poisson")

summary(ciona_disruption_glm_poisson)

ciona_disruption_glm_quasipoisson <- glm(COUNTresponse ~  0 + Flowrate * SurfaceType, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Ciona",], family="quasipoisson")

summary(ciona_disruption_glm_quasipoisson)

ciona_disruption_glm_NB <- glm.nb(COUNTresponse ~  0 + Flowrate * SurfaceType, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Ciona",], link = "log")

summary(ciona_disruption_glm_NB)

#emmip(ciona_disruption_glm_NB, Flowrate ~  SurfaceType)
#plot_model(ciona_disruption_glm_NB, type = "pred", terms = c("Flowrate", "SurfaceType"))

ciona_disruption_glm_NB_FR <- glm.nb(COUNTresponse ~   0 + Flowrate, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Ciona" & bubble_lab_disruption$SurfaceType=="FR-IS1100",], link="log")

ciona_disruption_glm_NB_ACR <- glm.nb(COUNTresponse ~  0 +  Flowrate, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Ciona" & bubble_lab_disruption$SurfaceType=="ACR",], link="log")

summary(ciona_disruption_glm_NB_FR)
summary(ciona_disruption_glm_NB_ACR)


#Crassostrea
Crassostrea_disruption_glm_poisson <- glm(COUNTresponse ~ 0 +  Flowrate * SurfaceType, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Crassostrea",], family="poisson")

summary(Crassostrea_disruption_glm_poisson)

Crassostrea_disruption_glm_quasipoisson <- glm(COUNTresponse ~  0 + Flowrate * SurfaceType, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Crassostrea",], family="quasipoisson")

summary(Crassostrea_disruption_glm_quasipoisson)

Crassostrea_disruption_glm_NB <- glm.nb(COUNTresponse ~  0 + Flowrate * SurfaceType, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Crassostrea",], link = "log")

summary(Crassostrea_disruption_glm_NB)

#plot the three way interaction
#emmip(Crassostrea_scouring_glm_NB, Flowrate ~ Settlementperiodhrs | SurfaceType)
#emmip(Crassostrea_scouring_glm_NB, Settlementperiodhrs ~ Flowrate | SurfaceType)

#plot_model(Crassostrea_scouring_glm_NB, type = "pred", terms = c("Flowrate", "Settlementperiodhrs", "SurfaceType"))

Crassostrea_disruption_glm_NB_FR <- glm.nb(COUNTresponse ~   0 + Flowrate, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Crassostrea" & bubble_lab_disruption$SurfaceType=="FR-IS1100",], link="log")

Crassostrea_disruption_glm_NB_ACR <- glm.nb(COUNTresponse ~  0 +  Flowrate, data=bubble_lab_disruption[bubble_lab_disruption$Organism=="Crassostrea" & bubble_lab_disruption$SurfaceType=="ACR",], link="log")

summary(Crassostrea_disruption_glm_NB_FR)
summary(Crassostrea_disruption_glm_NB_ACR)

```

#### Field trials

##### Vertical surfaces 

We performed a ordinal mixed effect model (Agresti 2002) to test how surface type (CONC or FR-IS1100) and the bubble treatment affect the levels of fouling (LoF) over the experimental period. We included a interaction term in the models between surface type and bubble treatment to test for a combined effect fo these factors. Because of the repetitions of the observations over 13 weeks, we added the different weeks at which the observations were taken as a random effect on the models to account for the effect of time in the results. To test if there were differences between the different factor combinations of the models, we calculated the predicted marginal means from the models and compared them with a *z*-test. 

We performed a GLM (with Gamma error distribution) to test how surface type and bubble treatment affect the end biomass of fouling present after the experimental period (assessed as dry weight). We included a interaction term in the models between surface type and bubble treatment to test for a combined effect fo these factors.

We performed a beta regresion to test how surface type and bubble treatment affect the end percent cover of macro-fouling present after the experimental period (assessed by point counts of the plate photographs). We included a interaction term in the models between surface type and bubble treatment to test for a combined effect fo these factors.

Agresti, A. (2002) Categorical Data Analysis. Second edition. Wiley.


```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
#########################
# Levels of fouling (LOF)

#summary(field_vertical)

names(field_vertical)

##########
# LOF plot

bubble_pilot_LOF_melt <- melt(field_vertical[,c("Treatment", "Frame", "Panel", "LoF_week0", "LoF_week1", "LoF_week4", "LoF_week6", "LoF_week9", "LoF_week13" )], id.vars = c("Treatment", "Frame", "Panel"))

bubble_pilot_LOF_melt$Treatment <- ifelse(bubble_pilot_LOF_melt$Treatment=="CONTROL","Control","Bubbled")
bubble_pilot_LOF_melt$Surface <- ifelse(grepl("Concrete",bubble_pilot_LOF_melt$Panel), "CONC", "FR-IS1100")
bubble_pilot_LOF_melt$Week <- as.numeric(gsub("LoF_week","",bubble_pilot_LOF_melt$variable))

bubble_pilot_LOF_plot <- ggplot(bubble_pilot_LOF_melt, aes(y=value, x=as.factor(Week),colour=Treatment)) +     geom_point()+
    geom_line(aes(group=Treatment)) +
    #geom_boxplot(position="fill", stat="identity", alpha=0.75)  +
    scale_colour_manual("", values = c( "#999999","red")) +
    labs(y="Level of fouling", x="Week")  +
    facet_grid(cols = vars(Surface))  +
  theme_bw(base_size = 18) +
  theme(legend.key.size = unit(1, "cm"), axis.ticks = element_blank(),legend.position="bottom") 

#######################################
## Ordinal logistic regression for LOF

## Ordinal mixed model with week as random effect
library(ordinal)
bubble_LOF_clmm <- clmm(as.factor(value) ~ Surface * Treatment + (1|Week), data = bubble_pilot_LOF_melt[bubble_pilot_LOF_melt$Week!="0",])
summary(bubble_LOF_clmm)

#pairwise comparison
bubble_LOF_clmm_emm <- emmeans(bubble_LOF_clmm, ~Surface * Treatment)
bubble_LOF_clmm_pairs <- data.frame(pairs(bubble_LOF_clmm_emm))

## Ordinal regression model for repeated scores, including an autoregressive component (AR1)
 
# Fitted using a generalized estimating equation (GEE) methodology. The algorithm estimates the correlation parameter by minimizing the generalized variance of the regression parameters at each step of the fitting algorithm.

library(repolr)

bubble_pilot_LOF_melt$value <- as.integer(bubble_pilot_LOF_melt$value)
bubble_pilot_LOF_melt$Week <- as.integer(bubble_pilot_LOF_melt$Week)
bubble_pilot_LOF_melt$Treatment <- as.factor(bubble_pilot_LOF_melt$Treatment)
bubble_pilot_LOF_melt$Surface <- as.factor(bubble_pilot_LOF_melt$Surface)

unique_units <- unique(bubble_pilot_LOF_melt[,c("Frame","Panel","Treatment")])

bubble_pilot_LOF_melt$unit <- NA
for(i in 1:nrow(unique_units))
{
bubble_pilot_LOF_melt$unit <- ifelse(bubble_pilot_LOF_melt$Frame == unique_units$Frame[i] & bubble_pilot_LOF_melt$Panel == unique_units$Panel[i] & bubble_pilot_LOF_melt$Treatment == unique_units$Treatment[i], i, bubble_pilot_LOF_melt$unit)  
}  

bubble_pilot_LOF_melt$value <- bubble_pilot_LOF_melt$value+1
bubble_pilot_LOF_melt$Week <- bubble_pilot_LOF_melt$Week+1

bubble_LOF_w0 <- bubble_pilot_LOF_melt[,c("unit", "Surface", "Treatment", "Week","value")]
bubble_LOF_w0 <- bubble_LOF_w0[order(bubble_LOF_w0$unit, bubble_LOF_w0$Week),]

unique(bubble_LOF_w0$unit)
unique(bubble_LOF_w0$value)
unique(bubble_LOF_w0$Week)

bubble_LOF_w0$Treatment_Surface <- as.factor(paste(bubble_LOF_w0$Treatment,bubble_LOF_w0$Surface))

bubble_LOF_treatment <- repolr(value~Treatment, data=bubble_LOF_w0, categories=6, subjects="unit",times=c(1,2,5,7,10,14), corr.mod="ar1", alpha=0.5)

bubble_LOF_surface <- repolr(value~Surface, data=bubble_LOF_w0, categories=6, subjects="unit",times=c(1,2,5,7,10,14), corr.mod="ar1", alpha=0.5)

#The interaction model doesn't run because it generate a singular matrix 
#bubble_LOF_Treatment_Surface <- repolr(value~Treatment*Surface, data=bubble_LOF_w0, categories=6, subjects="unit",times=c(1,2,5,7,10,14), corr.mod="ar1", alpha=0.5)  

#summary(bubble_LOF_treatment)
#summary(bubble_LOF_surface)

###############
# PERCENT COVER

bubble_pilot_percent <- field_vertical[,c("Treatment", "Frame", "Panel", "Notes", "Point_count_Number_of_points_classified_in_image", "Point_count_Biofilm", "Point_count_Bare_space")]

bubble_pilot_percent$Point_count_macrofouling <- 100 - (bubble_pilot_percent$Point_count_Biofilm + bubble_pilot_percent$Point_count_Bare_space)

###############
## Percent cover plot

bubble_pilot_percent <- bubble_pilot_percent[bubble_pilot_percent$Notes == "",]

bubble_pilot_percent$Treatment <- ifelse(bubble_pilot_percent$Treatment=="CONTROL","Control","Bubbled")
bubble_pilot_percent$Surface <- ifelse(grepl("Concrete",bubble_pilot_percent$Panel), "CONC", "FR-IS1100")

bubble_pilot_percent_plot <- ggplot(bubble_pilot_percent, aes(y=Point_count_macrofouling, x=Treatment, colour=Treatment)) + 
    geom_boxplot()+
    #geom_boxplot(position="fill", stat="identity", alpha=0.75)  +
    scale_colour_manual("", values = c( "#999999","red")) +
    labs(y="Macrofouling cover (%)", x="")  +
    facet_grid(cols = vars(Surface))  +
  theme_bw(base_size = 18) +
  theme(legend.position="none") 


bubble_pilot_percent_melt <- melt(bubble_pilot_percent[,c("Treatment", "Frame", "Surface", "Point_count_Biofilm", "Point_count_Bare_space", "Point_count_macrofouling" )], id.vars = c("Treatment", "Frame", "Surface"))

bubble_pilot_percent_plot2 <- ggplot(bubble_pilot_percent_melt, aes(y=value, x=Treatment, colour=variable)) + 
    geom_boxplot()+
    scale_colour_manual("", values = c("#000000ff","#999999","red"),
                       labels=c("Biofilm", "Bare space", "Macrofouling")) +
    labs(y="% cover of macrofouling", x="")  +
    facet_grid(cols = vars(Surface))  +
  theme_bw(base_size = 18) +
  theme(legend.position="bottom") 

#########################################
## Beta regresion for percent cover of Macrofouling

# Using beta regression
library(betareg)

bubble_pilot_percent$macrofouling <- (bubble_pilot_percent$Point_count_macrofouling/100)+0.001
summary(bubble_pilot_percent)

percent_macrofouling_beta <- betareg(macrofouling ~ Surface * Treatment, data = bubble_pilot_percent)

summary(percent_macrofouling_beta)

#pairwise comparison
percent_macrofouling_beta_emm <- emmeans(percent_macrofouling_beta, ~Surface * Treatment)
percent_macrofouling_emm_pairs <- data.frame(pairs(percent_macrofouling_beta_emm))

#########
# BIOMASS

bubble_pilot_biomass <- field_vertical[,c("Treatment", "Frame", "Panel", "Notes", "Biomass_dry_weight")]
bubble_pilot_biomass <- bubble_pilot_biomass[bubble_pilot_biomass$Notes == "",]

###############
## Biomass plot

bubble_pilot_biomass$Treatment <- ifelse(bubble_pilot_biomass$Treatment=="CONTROL","Control","Bubbled")
bubble_pilot_biomass$Surface <- ifelse(grepl("Concrete",bubble_pilot_biomass$Panel), "CONC", "FR-IS1100")

bubble_pilot_biomass_plot <- ggplot(bubble_pilot_biomass, aes(y=Biomass_dry_weight, x=Treatment, colour=Treatment)) + 
    geom_boxplot()+
    #geom_boxplot(position="fill", stat="identity", alpha=0.75)  +
    scale_colour_manual("", values = c( "#999999","red")) +
    labs(y="Dry weight (g)", x="")  +
    facet_grid(cols = vars(Surface))  +
  theme_bw(base_size = 18) +
  theme(legend.position="none") 

##################
## GLM for biomass

summary(bubble_pilot_biomass)

shapiro.test(bubble_pilot_biomass$Biomass_dry_weight)
ggdensity(bubble_pilot_biomass$Biomass_dry_weight)
gghistogram(bubble_pilot_biomass$Biomass_dry_weight)

# Using Gamma family, useful for highly positively skewed data

Biomass_glm <- glm(Biomass_dry_weight ~ Surface * Treatment, data = bubble_pilot_biomass, family = "Gamma")
summary(Biomass_glm)

#pairwise comparison
Biomass_glm_emm <- emmeans(Biomass_glm, ~Surface * Treatment)
Biomass_glm_emm_pairs <- data.frame(pairs(Biomass_glm_emm))
```

##### Horizontal and angled surfaces 

For the first three "rounds" of the trials, we performed generalized linear models to test how surface type and orientation angle affect the efficiency of the bubble treatment on reducing the percent cover of biofouling present after the experimental period. We included a interaction term in the models between surface type and orientation angle to test for a combined effect fo these factors. Because the repetitions of the trial had to be carried out at different dates, we added the different "rounds" as another fixed effect on the generalized linear models to account for the effect of seasonality in the results. We performed a differnt model for each of the components of percent cover measured: barespace, biofilm and  macrofouling. For the final "round" of the trial, when the effect of bubbles was tested against a control for only two of the surface types, we inspected the results graphically. We separately model the percent cover of barespace, biofilm and macrofouling using a binomial generalized model. To test if there were differences between the different factor combinations of the models, we calculated the predicted marginal means from the models and compared them with a *z*-test. 
 
All statistical analyses were performed within the 'R' statistical and programming environment (RCoreTeam 2019). We use the package emmeans (Lenth 2019) to estimate the marginal means (least-squares means) for the factor combinations from the generalized linear models, and compare these means. 

Dobson, A. J. (1990) An Introduction to Generalized Linear Models. London: Chapman and Hall.

R Core Team (2019). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL
  https://www.R-project.org/.
  
Russell Lenth (2019). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.3.5.1.
  https://CRAN.R-project.org/package=emmeans


```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
## field trial for horizontal and angled surfaces

summary(field_horizontal_angled)

field_horizontal_angled_melt <- melt(field_horizontal_angled[,2:8], id.vars = c("Round", "Treatment", "Surface", "Angle"))

field_horizontal_angled_melt$Surface <- ifelse(field_horizontal_angled_melt$Surface=="IS1000", "FR-IS1000",field_horizontal_angled_melt$Surface )
field_horizontal_angled_melt$Surface <- ifelse(field_horizontal_angled_melt$Surface=="IS1100", "FR-IS1100",field_horizontal_angled_melt$Surface )

field_horizontal_angled_melt$variable <- as.character(field_horizontal_angled_melt$variable)
field_horizontal_angled_melt$variable <- ifelse(field_horizontal_angled_melt$variable=="Barespace", "Bare space",
field_horizontal_angled_melt$variable)

# Percent cover for R1, R2 and R3 

FirstRounds_Percent_cover_plot_v1 <- ggplot(field_horizontal_angled_melt[field_horizontal_angled_melt$Round!="Final",], aes(fill=variable, y=value, x=Round)) + 
    geom_bar(position="fill", stat="identity", alpha=0.75)  +
    scale_fill_manual("", values = c( "#999999","#000000ff","red")) +
    labs(y="Percent cover", x="Round")  +
    facet_grid(rows = vars(Angle), cols = vars(Surface))  +
  theme_bw(base_size = 18) +
  theme(legend.key.size = unit(1, "cm"), axis.ticks = element_blank(),legend.position="bottom") 

FirstRounds_Percent_cover_plot_v2 <- ggplot(field_horizontal_angled_melt[field_horizontal_angled_melt$Round!="Final",], aes(fill=variable, y=value, x=Round)) + 
    geom_bar(position="fill", stat="identity", alpha=0.75)  +
    scale_fill_manual("", values = c( "#999999","#000000ff","red")) +
    labs(y="Percent cover", x="Round")  +
    facet_grid(rows = vars(Surface), cols = vars(Angle))  +
  theme_bw(base_size = 18) +
  theme(legend.key.size = unit(1, "cm"), axis.ticks = element_blank(),legend.position="bottom") 

FirstRounds_Percent_cover_plot_v3 <- ggplot(field_horizontal_angled_melt[field_horizontal_angled_melt$Round!="Final",], aes(fill=variable, y=value, x=Angle)) + 
    geom_bar(position="fill", stat="identity", alpha=0.75)  +
    scale_fill_manual("", values = c( "#999999","#000000ff","red")) +
    labs(y="Percent cover", x="Angle")  +
    facet_grid(rows = vars(Surface), cols = vars(Round))  +
  theme_bw(base_size = 18) +
  theme(legend.key.size = unit(1, "cm"), axis.ticks = element_blank(),legend.position="bottom") 

# Percent cover for Final

FinalRound_Percent_cover_plot <- ggplot(field_horizontal_angled_melt[field_horizontal_angled_melt$Round=="Final",], aes(fill=variable, y=value, x=Treatment)) + 
    geom_bar(position="fill", stat="identity", alpha=0.75)  +
    scale_fill_manual("", values = c( "#999999","#000000ff","red")) +
    labs(y="Percent cover", x="Treatment")  +
    facet_grid(cols = vars(Surface))  +
  theme_bw(base_size = 18) +
  theme(legend.key.size = unit(1, "cm"), axis.ticks = element_blank(),legend.position="bottom") 


################
##R1, R2 and R3

#Barespace plot

ggplot_Barespace_field <- ggplot(aes(y = Barespace, x = Surface), data = field_horizontal_angled[field_horizontal_angled$Round!="Final",]) +
  geom_boxplot(size=1.2, alpha=0.4, fill="#999999") +
  labs(y="Percent of barespace", x="Surface type")  +
  theme_bw(base_size = 22) +
  theme(legend.key.size = unit(2, "cm"), axis.ticks = element_blank()) +
  facet_grid(rows = vars(Angle), scales = "free_y", switch = 'x') 

#Barespace analysis

field_horizontal_angled$BarespaceInteger <- field_horizontal_angled$Barespace *100

Barespace_glm_binomial <- glm(cbind(BarespaceInteger, 100-BarespaceInteger) ~ Surface * Angle + Round, data = field_horizontal_angled[field_horizontal_angled$Round!="Final",], family = binomial(logit))

summary(Barespace_glm_binomial)


#pairwise comparison
Barespace_glm_binomial_emm <- emmeans(Barespace_glm_binomial, ~Surface * Angle)
Barespace_glm_binomial_emm_pairs <- data.frame(pairs(Barespace_glm_binomial_emm))

Barespace_glm_binomial_emm_pairs <- Barespace_glm_binomial_emm_pairs[order(Barespace_glm_binomial_emm_pairs$p.value),]
Barespace_glm_binomial_emm_pairs$p.value <- round(Barespace_glm_binomial_emm_pairs$p.value,5)

#Biofilm plot

ggplot_Biofilm_field <- ggplot(aes(y = Biofilm, x = Surface), data = field_horizontal_angled[field_horizontal_angled$Round!="Final",]) +
  geom_boxplot(size=1.2, alpha=0.4, fill="#999999") +
  labs(y="Percent of biofilm", x="Surface type")  +
  theme_bw(base_size = 22) +
  theme(legend.key.size = unit(2, "cm"), axis.ticks = element_blank()) +
  facet_grid(rows = vars(Angle), scales = "free_y", switch = 'x') 

#Biofilm analysis

field_horizontal_angled$BiofilmInteger <- field_horizontal_angled$Biofilm *100

Biofilm_glm_binomial <- glm(cbind(BiofilmInteger, 100-BiofilmInteger) ~ Surface * Angle  + Round, data = field_horizontal_angled[field_horizontal_angled$Round!="Final",], family = binomial(logit))

summary(Biofilm_glm_binomial)

#pairwise comparison
Biofilm_glm_binomial_emm <- emmeans(Biofilm_glm_binomial, ~Surface * Angle)
Biofilm_glm_binomial_emm_pairs <- data.frame(pairs(Biofilm_glm_binomial_emm))

Biofilm_glm_binomial_emm_pairs <- Biofilm_glm_binomial_emm_pairs[order(Biofilm_glm_binomial_emm_pairs$p.value),]
Biofilm_glm_binomial_emm_pairs$p.value <- round(Biofilm_glm_binomial_emm_pairs$p.value,5)

#Macrofouling plot

ggplot_Macrofouling_field <- ggplot(aes(y = Macrofouling, x = Surface), data = field_horizontal_angled[field_horizontal_angled$Round=="Final",]) +
  geom_point(size=4, alpha=0.8, fill="#999999",position = position_dodge2(width = 1)) +
  labs(y="percent of Macrofouling", x="surface type")  +
  theme_bw(base_size = 22) +
  theme(legend.key.size = unit(2, "cm"), axis.ticks = element_blank()) +
  facet_grid(cols = vars(Treatment), scales = "free_y", switch = 'x') 

#Macrofouling analysis

field_horizontal_angled$MacrofoulingInteger <- field_horizontal_angled$Macrofouling *100

Macrofouling_glm_binomial <- glm(cbind(MacrofoulingInteger, 100-MacrofoulingInteger) ~ Surface * Angle + Round, data = field_horizontal_angled[field_horizontal_angled$Round!="Final",], family = binomial(logit))

summary(Macrofouling_glm_binomial)

#pairwise comparison
Macrofouling_glm_binomial_emm <- emmeans(Macrofouling_glm_binomial, ~Surface * Angle)
Macrofouling_glm_binomial_emm_pairs <- data.frame(pairs(Macrofouling_glm_binomial_emm))

Macrofouling_glm_binomial_emm_pairs <- Macrofouling_glm_binomial_emm_pairs[order(Macrofouling_glm_binomial_emm_pairs$p.value),]
Macrofouling_glm_binomial_emm_pairs$p.value <- round(Macrofouling_glm_binomial_emm_pairs$p.value,5)

#######
##Final

#Barespace plot

ggplot_Barespace_field2 <- ggplot(aes(y = Barespace, x = Surface), data = field_horizontal_angled[field_horizontal_angled$Round=="Final",]) +
  geom_boxplot(size=1.2, alpha=0.4, fill="#999999") +
  labs(y="percent of barespace", x="surface type")  +
  theme_bw(base_size = 22) +
  theme(legend.key.size = unit(2, "cm"), axis.ticks = element_blank()) +
  facet_grid(cols = vars(Treatment), scales = "free_y", switch = 'x') 

#Barespace analysis

field_horizontal_angled$BarespaceInteger <- field_horizontal_angled$Barespace *100

Barespace_glm_binomial2 <- glm(cbind(BarespaceInteger, 100-BarespaceInteger) ~ Surface * Treatment, data = field_horizontal_angled[field_horizontal_angled$Round=="Final",], family = binomial(logit))

summary(Barespace_glm_binomial2)

#pairwise comparison
Barespace_glm_binomial_emm2 <- emmeans(Barespace_glm_binomial2, ~ Surface * Treatment)
Barespace_glm_binomial_emm_pairs2 <- data.frame(pairs(Barespace_glm_binomial_emm2))

Barespace_glm_binomial_emm_pairs2 <- Barespace_glm_binomial_emm_pairs2[order(Barespace_glm_binomial_emm_pairs2$p.value),]
Barespace_glm_binomial_emm_pairs2$p.value <- round(Barespace_glm_binomial_emm_pairs2$p.value,5)

#Biofilm plot

ggplot_Biofilm_field2 <- ggplot(aes(y = Biofilm, x = Surface), data = field_horizontal_angled[field_horizontal_angled$Round=="Final",]) +
  geom_boxplot(size=1.2, alpha=0.4, fill="#999999") +
  labs(y="percent of Biofilm", x="surface type")  +
  theme_bw(base_size = 22) +
  theme(legend.key.size = unit(2, "cm"), axis.ticks = element_blank()) +
  facet_grid(cols = vars(Treatment), scales = "free_y", switch = 'x') 

#Biofilm analysis

field_horizontal_angled$BiofilmInteger <- field_horizontal_angled$Biofilm *100

Biofilm_glm_binomial2 <- glm(cbind(BiofilmInteger, 100-BiofilmInteger) ~ Surface * Treatment, data = field_horizontal_angled[field_horizontal_angled$Round=="Final",], family = binomial(logit))

summary(Biofilm_glm_binomial2)

#pairwise comparison
Biofilm_glm_binomial_emm2 <- emmeans(Biofilm_glm_binomial2, ~Surface * Treatment)
Biofilm_glm_binomial_emm_pairs2 <- data.frame(pairs(Biofilm_glm_binomial_emm2))

Biofilm_glm_binomial_emm_pairs2 <- Biofilm_glm_binomial_emm_pairs2[order(Biofilm_glm_binomial_emm_pairs2$p.value),]
Biofilm_glm_binomial_emm_pairs2$p.value <- round(Biofilm_glm_binomial_emm_pairs2$p.value,5)

#Macrofouling plot

ggplot_Macrofouling_field2 <- ggplot(aes(y = Macrofouling, x = Surface), data = field_horizontal_angled[field_horizontal_angled$Round=="Final",]) +
  geom_boxplot(size=1.2, alpha=0.4, fill="#999999") +
   labs(y="Percent of macrofouling", x="surface type")  +
  theme_bw(base_size = 22) +
  theme(legend.key.size = unit(2, "cm"), axis.ticks = element_blank()) +
  facet_grid(cols = vars(Treatment), scales = "free_y", switch = 'x') 

#Macrofouling analysis

field_horizontal_angled$MacrofoulingInteger <- field_horizontal_angled$Macrofouling *100

Macrofouling_glm_binomial2 <- glm(cbind(MacrofoulingInteger, 100-MacrofoulingInteger) ~ Surface * Treatment, data = field_horizontal_angled[field_horizontal_angled$Round=="Final",], family = binomial(logit))

summary(Macrofouling_glm_binomial2)

#pairwise comparison
Macrofouling_glm_binomial_emm2 <- emmeans(Macrofouling_glm_binomial2, ~Surface * Treatment)
Macrofouling_glm_binomial_emm_pairs2 <- data.frame(pairs(Macrofouling_glm_binomial_emm))

Macrofouling_glm_binomial_emm_pairs2 <- Macrofouling_glm_binomial_emm_pairs2[order(Macrofouling_glm_binomial_emm_pairs2$p.value),]
Macrofouling_glm_binomial_emm_pairs2$p.value <- round(Macrofouling_glm_binomial_emm_pairs2$p.value,5)
```

### Results

#### Laboratory trials

The results of the scouring trial showed that bubble high flow rate reduced the number of Ciona individuals at both 3 and 120 hours settled on the acrylic surface (Figure 1, Table 1 and 2). In contrast, there were no differences between the number of Crassostrea individuals settled on the acrylic surface at any flow rate, regardless of settlement time. However, all bubble flow rates reduced the number of Crassostrea individuals compared with the Nil control when individuals were settled for 120 hrs on the FR surface (Figure 1, Table 1 and 2). The disruption trial showed that a medium flow rate decreased the number of Crassostrea individuals to almost zero in both surface types tested when compared with the no bubble treatment. For Ciona, there was a significant decrease in individual numbers in the medium flow treatment compared with the no bubble treatment, but this difference was more obvious in the acrilic surface (Figure 2, Table 3). 

##### Figure 1. Ciona and Crassostrea scouring plot

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}

#Ciona and Crassostrea scouring plot
ggplot_bubble_scouring

```

<br><br><br><br>

##### Table 1. Results of the negative binomial models for the scouring trial

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
#Ciona scouring results

ciona_scouring_glm_NB_summary <- summary(ciona_scouring_glm_NB)

#Crassostrea scouring results

Crassostrea_scouring_glm_NB_summary <- summary(Crassostrea_scouring_glm_NB)

Crassostrea_scouring_glm_NB_120_FR_summary <- summary(Crassostrea_scouring_glm_NB_120_FR)

ciona_scouring_glm_NB_summary_coef <- data.frame(ciona_scouring_glm_NB_summary$coefficients)
Crassostrea_scouring_glm_NB_summary_coef <- data.frame(Crassostrea_scouring_glm_NB_summary$coefficients)
Crassostrea_scouring_glm_NB_120_FR_summary_coef <- data.frame(Crassostrea_scouring_glm_NB_120_FR_summary$coefficients)

scouring_results <- rbind(ciona_scouring_glm_NB_summary_coef,
                       Crassostrea_scouring_glm_NB_summary_coef,
                       Crassostrea_scouring_glm_NB_120_FR_summary_coef)

scouring_results$model <- c(rep(paste(as.character(ciona_scouring_glm_NB_summary$call$formula)[2:3], collapse = " ~ "),8),
                                         rep(paste(as.character(Crassostrea_scouring_glm_NB_summary$call$formula)[2:3], collapse = " ~ "),16),
                                         rep(paste(as.character(Crassostrea_scouring_glm_NB_120_FR_summary$call$formula)[2:3], collapse = " ~ "),4))

scouring_results$data <- c(rep("Ciona, acrylic surface, 3 and 120 settelment hours",8),
                           rep("Crassostrea, acrylic and FR surfaces, 3 and 120 settelment hours",16),
                           rep("Crassostrea, FR surface, 120 settelment hours",4))

scouring_results$factor <- row.names(scouring_results)


scouring_results <- scouring_results[c("data","model","factor","Estimate","z.value","Pr...z..")]

scouring_results$Estimate <- round(scouring_results$Estimate,3)
scouring_results$z.value <- round(scouring_results$z.value,3)

scouring_results$Pr...z.. <- round(scouring_results$Pr...z..,5)
names(scouring_results)[6] <- "p.value"

scouring_results$model <- gsub("COUNTresponse ", "", scouring_results$model)

scouring_results_gropued <- as_grouped_data(x = scouring_results, groups = c("data","model"))
#head(scouring_results_gropued)

Table_1_prety <-  as_flextable(scouring_results_gropued) %>% 
   bold(j = 1, i = ~ !is.na(data), bold = TRUE, part = "body" ) %>%
   bold(j = 1, i = ~ !is.na(model), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 1. Results of the negative binomial models for the scouring trial.") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(4,1,1,1)) %>%
   theme_booktabs()

#loop to colorise the columns

columns <- names(scouring_results_gropued)[3:6]

#anything with a p.value of more less than 0.05 is formatted light red
for (i in seq_along(columns)){
eval(parse(text= paste0("Table_1_prety<-", "bg (Table_1_prety,i=~", "p.value", "<0.05, j=~", columns[i], " ,bg=\"#ffd5d5ff\")")))
}

Table_1_prety <- fontsize(Table_1_prety, part = "all", size = 8) 

Table_1_prety
```

<br><br><br><br>

##### Table 2. Pairwise comparisons of the treatment combinations for the scouring trial.

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}

ciona_scouring_emm_pairs$species <- "Ciona"
Crassostrea_scouring_emm_pairs$species <- "Crassostrea"

scouring_emm_pairs <- rbind(ciona_scouring_emm_pairs, Crassostrea_scouring_emm_pairs)

scouring_emm_pairs_gropued <- as_grouped_data(x = scouring_emm_pairs[scouring_emm_pairs$p.value<0.05,-4], groups = c("species"))
#head(scouring_results_gropued)

Table_2_prety <-  as_flextable(scouring_emm_pairs_gropued) %>% 
   bold(j = 1, i = ~ !is.na(species), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 2. Pairwise comparisons of the treatment combinations for the scouring trial (only showing significant differences).") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(2,1,1,1,1)) %>%
   theme_booktabs()

#loop to colorise the columns

#columns <- names(scouring_emm_pairs_gropued)[2:6]

#anything with a p.value of more less than 0.05 is formatted light red

# for (i in seq_along(columns)){
# eval(parse(text= paste0("Table_2_prety<-", "bg (Table_2_prety,i=~", "p.value", "<0.05, j=~", columns[i], " ,bg=\"#ffd5d5ff\")")))
# }

Table_2_prety <- fontsize(Table_2_prety, part = "all", size = 8) 

Table_2_prety
```

<br><br><br><br>

##### Figure 2. Ciona and Crassostrea disruption plot

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}

#Ciona and Crassostrea disruption plot
ggplot_bubble_disruption

```

<br><br><br><br>

##### Table 1. Results of the negative binomial models for the scouring trial

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
#Ciona disruption results

ciona_disruption_glm_NB_summary <- summary(ciona_disruption_glm_NB)

ciona_disruption_glm_NB_FR_summary <- summary(ciona_disruption_glm_NB_FR)
ciona_disruption_glm_NB_ACR_summary <- summary(ciona_disruption_glm_NB_ACR)

#Crassostrea disruption results

Crassostrea_disruption_glm_NB_summary <- summary(Crassostrea_disruption_glm_NB)

Crassostrea_disruption_glm_NB_FR_summary <- summary(Crassostrea_disruption_glm_NB_FR)
Crassostrea_disruption_glm_NB_ACR_summary <- summary(Crassostrea_disruption_glm_NB_ACR)

#model coefficients

ciona_disruption_glm_NB_summary_coef <- data.frame(ciona_disruption_glm_NB_summary$coefficients)
ciona_disruption_glm_NB_FR_summary_coef <- data.frame(ciona_disruption_glm_NB_FR_summary$coefficients)
ciona_disruption_glm_NB_ACR_summary_coef <- data.frame(ciona_disruption_glm_NB_ACR_summary$coefficients)

Crassostrea_disruption_glm_NB_summary_coef <- data.frame(Crassostrea_disruption_glm_NB_summary$coefficients)
Crassostrea_disruption_glm_NB_FR_summary_coef <- data.frame(Crassostrea_disruption_glm_NB_FR_summary$coefficients)
Crassostrea_disruption_glm_NB_ACR_summary_coef <- data.frame(Crassostrea_disruption_glm_NB_ACR_summary$coefficients)

disruption_results <- rbind(ciona_disruption_glm_NB_summary_coef,
                          ciona_disruption_glm_NB_FR_summary_coef,
                          ciona_disruption_glm_NB_ACR_summary_coef,
                          Crassostrea_disruption_glm_NB_summary_coef,
                          Crassostrea_disruption_glm_NB_FR_summary_coef,
                          Crassostrea_disruption_glm_NB_ACR_summary_coef)

disruption_results$model <- c(rep(paste(as.character(ciona_disruption_glm_NB_summary$call$formula)[2:3], collapse = " ~ "),4),
                              rep(paste(as.character(ciona_disruption_glm_NB_FR_summary$call$formula)[2:3], collapse = " ~ "),2),
                              rep(paste(as.character(ciona_disruption_glm_NB_ACR_summary$call$formula)[2:3], collapse = " ~ "),2),
                              
                              rep(paste(as.character(Crassostrea_disruption_glm_NB_summary$call$formula)[2:3], collapse = " ~ "),4),
                              rep(paste(as.character(Crassostrea_disruption_glm_NB_FR_summary$call$formula)[2:3], collapse = " ~ "),2),
                              rep(paste(as.character(Crassostrea_disruption_glm_NB_ACR_summary$call$formula)[2:3], collapse = " ~ "),2))

disruption_results$species <- c(rep("Ciona",8),
                              rep("Crassostrea",8))

disruption_results$factor <- row.names(disruption_results)

disruption_results <- disruption_results[c("species","model","factor","Estimate","z.value","Pr...z..")]

disruption_results$Estimate <- round(disruption_results$Estimate,3)
disruption_results$z.value <- round(disruption_results$z.value,3)

disruption_results$Pr...z.. <- round(disruption_results$Pr...z..,5)
names(disruption_results)[6] <- "p.value"

disruption_results$model <- gsub("COUNTresponse ~ 0 +", "~", disruption_results$model)
disruption_results$model <- gsub("~+", "", disruption_results$model)

disruption_results_gropued <- as_grouped_data(x = disruption_results, groups = c("species","model"))
#head(scouring_results_gropued)

Table_3_prety <-  as_flextable(disruption_results_gropued) %>% 
   bold(j = 1, i = ~ !is.na(species), bold = TRUE, part = "body" ) %>%
   bold(j = 1, i = ~ !is.na(model), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 1. Results of the negative binomial models for the disruption trial.") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(4,1,1,1)) %>%
   theme_booktabs()

#loop to colorise the columns

columns <- names(disruption_results_gropued)[3:6]

#anything with a p.value of more less than 0.05 is formatted light red
for (i in seq_along(columns)){
eval(parse(text= paste0("Table_3_prety<-", "bg (Table_3_prety,i=~", "p.value", "<0.05, j=~", columns[i], " ,bg=\"#ffd5d5ff\")")))
}

Table_3_prety <- fontsize(Table_3_prety, part = "all", size = 8) 

Table_3_prety
```

<br><br><br><br>

#### Field trials

##### Vertical surfaces

The model results looking at the effect of CONC and FR-IS1100 surfaces and bubble treatment on the level of fouling (LOF) showed significant differences among the surfaces and bubble treatment, including a significant interaction between the different surfaces and bubble treatment and control (all p-values < 0.05). Also, the model analysing the dry weight of fouling in the panels at the end of the experiment showed significant differences among the CONC and FR-IS1100 surfaces and bubble treatment (all p-values < 0.05 including a significant interation between surfaces and bubble treatment).

###### Figure 3. Level of fouling and biomass accumulation for the vertical surfaces field trial.

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
library(cowplot)
plot_grid(bubble_pilot_biomass_plot, bubble_pilot_percent_plot, bubble_pilot_LOF_plot, nrow = 3, labels = "AUTO", rel_heights = c(0.75, 0.75,1), ncol = 1, align = 'v')
```

<br><br><br><br>

###### Table 4. Level of fouling model results

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
LOF_model_summary <- summary(bubble_LOF_clmm)

#model coefficients

LOF_model_summary_coef <- data.frame(LOF_model_summary$coefficients[6:8,])

field_results <- rbind(LOF_model_summary_coef)

field_results$model <- c(rep(paste(as.character(LOF_model_summary$call$formula)[2:3], collapse = " ~ "),3))

field_results$factor <- row.names(field_results)

field_results <-field_results[c("model","factor","Estimate","z.value","Pr...z..")]

field_results$Estimate <- round(field_results$Estimate,3)
field_results$z.value <- round(field_results$z.value,3)

field_results$Pr...z.. <- round(field_results$Pr...z..,5)
names(field_results)[5] <- "p.value"

field_results$p.value <- ifelse(field_results$p.value<=0.001, "< 0.001", field_results$p.value)

field_results[field_results$model=="as.factor(value) ~ Surface * Treatment + (1 | Week)",]$model <- "LOF ~ Surface * Treatment + (1 | Week)"

field_results_gropued <- as_grouped_data(x = field_results, groups = c("model"))

Table_4_prety <-  as_flextable(field_results_gropued) %>% 
   bold(j = 1, i = ~ !is.na(model), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 4. Results of the generalized linear mixed model for the level of fouling.") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(4,1,1,1)) %>%
   theme_booktabs()

#loop to colorise the columns

columns <- names(field_results_gropued)[2:5]

#anything with a p.value of less than 0.05 is formatted light red
for (i in seq_along(columns)){
eval(parse(text= paste0("Table_4_prety<-", "bg (Table_4_prety,i=~", "p.value", "<0.05, j=~", columns[i], " ,bg=\"#ffd5d5ff\")")))
}

Table_4_prety <- fontsize(Table_4_prety, part = "all", size = 8) 

Table_4_prety
```

<br><br><br><br>

###### Table 5. Pairwise comparisons of the treatment combinations of the level of fouling results.

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
#pairwise comparison
bubble_LOF_clmm_pairs$estimate <- round(bubble_LOF_clmm_pairs$estimate,3)
bubble_LOF_clmm_pairs$SE <- round(bubble_LOF_clmm_pairs$SE,3)
bubble_LOF_clmm_pairs$z.ratio <- round(bubble_LOF_clmm_pairs$z.ratio,3)
bubble_LOF_clmm_pairs$p.value <- ifelse(bubble_LOF_clmm_pairs$p.value<0.001,"< 0.001", bubble_LOF_clmm_pairs$p.value)

Table_5_prety <-  flextable(bubble_LOF_clmm_pairs) %>% 
   align(part = "all") %>% # left align
   set_caption(caption = "Table 5. Pairwise comparisons of the treatment combinations for the levels of fouling (only showing significant differences).") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(3,1,1,1,1,1)) %>%
   theme_booktabs()

Table_5_prety <- fontsize(Table_5_prety, part = "all", size = 8) 

Table_5_prety
```

<br><br><br><br>

###### Table 6. Percent cover model results

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
percent_macrofouling_model_summary <- summary(percent_macrofouling_beta)

#model coefficients

percent_macrofouling_model_summary_coef <- data.frame(percent_macrofouling_model_summary$coefficients$mean[1:3,])

field_results$model <- c(rep(paste(as.character(percent_macrofouling_model_summary$call$formula)[2:3], collapse = " ~ "),3))

field_results$factor <- row.names(field_results)

field_results <-field_results[c("model","factor","Estimate","z.value","p.value")]

field_results$Estimate <- round(field_results$Estimate,3)
field_results$z.value <- round(field_results$z.value,3)

field_results[field_results$model=="macrofouling ~ Surface * Treatment",]$model <- "percent macrofouling cover ~ Surface * Treatment"

field_results_gropued <- as_grouped_data(x = field_results, groups = c("model"))

Table_6_prety <-  as_flextable(field_results_gropued) %>% 
   bold(j = 1, i = ~ !is.na(model), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 6. Results of the beta regresion for the percent macrofouling cover.") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(4,1,1,1)) %>%
   theme_booktabs()

#loop to colorise the columns

columns <- names(field_results_gropued)[2:5]

#anything with a p.value of less than 0.05 is formatted light red
for (i in seq_along(columns)){
eval(parse(text= paste0("Table_6_prety<-", "bg (Table_6_prety,i=~", "p.value", "<0.05, j=~", columns[i], " ,bg=\"#ffd5d5ff\")")))
}

Table_6_prety <- fontsize(Table_6_prety, part = "all", size = 8) 

Table_6_prety
```

<br><br><br><br>

###### Table 7. Pairwise comparisons of the treatment combinations of the percent cover.

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
#pairwise comparison

percent_macrofouling_emm_pairs$estimate <- round(percent_macrofouling_emm_pairs$estimate,3)
percent_macrofouling_emm_pairs$SE <- round(percent_macrofouling_emm_pairs$SE,3)
percent_macrofouling_emm_pairs$z.ratio <- round(percent_macrofouling_emm_pairs$z.ratio,3)
percent_macrofouling_emm_pairs$p.value <- ifelse(percent_macrofouling_emm_pairs$p.value<0.001,"< 0.001", as.character(round(as.numeric(percent_macrofouling_emm_pairs$p.value),4)))

Table_7_prety <-  flextable(percent_macrofouling_emm_pairs) %>% 
   align(part = "all") %>% # left align
   set_caption(caption = "Table 7. Pairwise comparisons of the treatment combinations for the macrofouling percent cover.") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(3,1,1,1,1,1)) %>%
   theme_booktabs()

Table_7_prety <- fontsize(Table_7_prety, part = "all", size = 8) 

Table_7_prety
```

###### Table 8. Biomass accumulation

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
biomass_model_summary <- summary(Biomass_glm)

#model coefficients

biomass_model_summary_coef <- data.frame(biomass_model_summary$coefficients)

field_results <- rbind(biomass_model_summary_coef)

field_results$model <- c(rep(paste(as.character(biomass_model_summary$call$formula)[2:3], collapse = " ~ "),4))

field_results$factor <- row.names(field_results)

field_results <-field_results[c("model","factor","Estimate","Std..Error","t.value","Pr...t..")]

field_results$Estimate <- round(field_results$Estimate,3)
field_results$t.value <- round(field_results$t.value,3)

field_results$Pr...z.. <- round(field_results$Pr...t..,5)
names(field_results)[5] <- "p.value"

field_results$Std..Error <- round(field_results$Std..Error,3)
names(field_results)[4] <- "Std.Error"

field_results$p.value <- ifelse(field_results$p.value<=0.001, "< 0.001", field_results$p.value)

field_results[field_results$model=="Biomass_dry_weight ~ Surface * Treatment",]$model <- "dry weight ~ Surface * Treatment"

field_results_gropued <- as_grouped_data(x = field_results[,c(1:5)], groups = c("model"))

Table_8_prety <-  as_flextable(field_results_gropued) %>% 
   bold(j = 1, i = ~ !is.na(model), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 8. Results of the generalized linear model for the biomass accumulation.") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(4,1,1,1)) %>%
   theme_booktabs()

#loop to colorise the columns

columns <- names(field_results_gropued)[2:5]

#anything with a p.value of less than 0.05 is formatted light red
for (i in seq_along(columns)){
eval(parse(text= paste0("Table_8_prety<-", "bg (Table_8_prety,i=~", "p.value", "<0.05, j=~", columns[i], " ,bg=\"#ffd5d5ff\")")))
}

Table_8_prety <- fontsize(Table_8_prety, part = "all", size = 8) 

Table_8_prety
```


<br><br><br><br>

###### Table 9. Pairwise comparisons of the treatment combinations of the biomass accumulation results.

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
#pairwise comparison

Biomass_glm_emm_pairs$estimate <- round(Biomass_glm_emm_pairs$estimate,3)
Biomass_glm_emm_pairs$SE <- round(Biomass_glm_emm_pairs$SE,3)
Biomass_glm_emm_pairs$z.ratio <- round(Biomass_glm_emm_pairs$z.ratio,3)
Biomass_glm_emm_pairs$p.value <- ifelse(Biomass_glm_emm_pairs$p.value<0.001,"< 0.001", Biomass_glm_emm_pairs$p.value)

Table_9_prety <-  flextable(Biomass_glm_emm_pairs) %>% 
   align(part = "all") %>% # left align
   set_caption(caption = "Table 9. Pairwise comparisons of the treatment combinations for the biommass model (only showing significant differences).") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(3,1,1,1,1,1)) %>%
   theme_booktabs()

Table_9_prety <- fontsize(Table_9_prety, part = "all", size = 8) 

Table_9_prety
```

##### Horizontal and angled surfaces

For the first three rounds of the trials, we found that there was an important temporal component on the results, which was more important on determining the percent cover of barespace, biofilm and macrofouling than the surface type or the position angle (Figure 3 and tables 4 and 5). However, the barespace and biofilm cover did show significant differences between some surface types and angles (Table 5). Figure 4 shows how the bubbled treatment reduced the percent cover of macrofouling and biofilm in the IS1000 surface, and reduced the macrofouling in the Poly sourface. In both cases, macrofouling was not detected in the bubbled treatment.

###### Figure 4. First Rounds Percent cover plot

There are three versions of this plot because I'm not sure which one would be best for ilustrating the story you are telling.

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}
#Percent cover for R1, R2 and R3 

FirstRounds_Percent_cover_plot_v1

FirstRounds_Percent_cover_plot_v2

FirstRounds_Percent_cover_plot_v3

```

<br><br><br><br>

###### Figure 5. Final Rounds Percent cover plot

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}

#Percent cover for Final

FinalRound_Percent_cover_plot

```

<br><br><br><br>

###### Table 10. First round percent cover model results

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}

Barespace_glm_binomial_summary <- summary(Barespace_glm_binomial)
Biofilm_glm_binomial_summary <- summary(Biofilm_glm_binomial)
Macrofouling_glm_binomial_summary <- summary(Macrofouling_glm_binomial)

#model coefficients

Barespace_glm_summary_coef <- data.frame(Barespace_glm_binomial_summary$coefficients)
Biofilm_glm_summary_coef <- data.frame(Biofilm_glm_binomial_summary$coefficients)
Macrofouling_glm_summary_coef <- data.frame(Macrofouling_glm_binomial_summary$coefficients)

field_results <- rbind(Barespace_glm_summary_coef,
                          Biofilm_glm_summary_coef,
                          Macrofouling_glm_summary_coef)

field_results$model <- c(rep(paste(as.character(Barespace_glm_binomial_summary$call$formula)[2:3], collapse = " ~ "),10),
                              rep(paste(as.character(Biofilm_glm_binomial_summary$call$formula)[2:3], collapse = " ~ "),10),
                              rep(paste(as.character(Macrofouling_glm_binomial_summary$call$formula)[2:3], collapse = " ~ "),10))

field_results$percent <- c(rep("Barespace",10),
                           rep("Biofilm",10),
                           rep("Macrofouling",10))

field_results$factor <- row.names(field_results)

field_results <-field_results[c("percent","model","factor","Estimate","z.value","Pr...z..")]

field_results$Estimate <- round(field_results$Estimate,3)
field_results$z.value <- round(field_results$z.value,3)

field_results$Pr...z.. <- round(field_results$Pr...z..,5)
names(field_results)[6] <- "p.value"

field_results[field_results$model=="cbind(BarespaceInteger, 100 - BarespaceInteger) ~ Surface * Angle + Round",]$model <- "Barespace ~ Surface * Angle + Round"
field_results[field_results$model=="cbind(BiofilmInteger, 100 - BiofilmInteger) ~ Surface * Angle + Round",]$model <- "Biofilm ~ Surface * Angle + Round"
field_results[field_results$model=="cbind(MacrofoulingInteger, 100 - MacrofoulingInteger) ~ Surface * Angle + Round",]$model <- "Macrofouling ~ Surface * Angle + Round"

field_results_gropued <- as_grouped_data(x = field_results, groups = c("percent","model"))

Table_10_prety <-  as_flextable(field_results_gropued) %>% 
   bold(j = 1, i = ~ !is.na(percent), bold = TRUE, part = "body" ) %>%
   bold(j = 1, i = ~ !is.na(model), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 10. Results of the first round percent cover model results.") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(4,1,1,1)) %>%
   theme_booktabs()

#loop to colorise the columns

columns <- names(field_results_gropued)[3:6]

#anything with a p.value of more less than 0.05 is formatted light red
for (i in seq_along(columns)){
eval(parse(text= paste0("Table_10_prety<-", "bg (Table_10_prety,i=~", "p.value", "<0.05, j=~", columns[i], " ,bg=\"#ffd5d5ff\")")))
}

Table_10_prety <- fontsize(Table_10_prety, part = "all", size = 8) 

Table_10_prety
```


<br><br><br><br>

###### Table 11. First rounds pairwise comparisons of the treatment combinations for the field trial.

<br><br><br><br>

```{r warning=FALSE, echo=TRUE, results=TRUE, fig.width=16, fig.height=12}

Barespace_glm_binomial_emm_pairs$percent <- "Barespace"
Biofilm_glm_binomial_emm_pairs$percent <- "Biofilm"
Macrofouling_glm_binomial_emm_pairs$percent <- "Macrofouling"

percent_emm_pairs <- rbind(Barespace_glm_binomial_emm_pairs, Biofilm_glm_binomial_emm_pairs,Macrofouling_glm_binomial_emm_pairs)

percent_emm_pairs_gropued <- as_grouped_data(x = percent_emm_pairs[percent_emm_pairs$p.value<0.05,-4], groups = c("percent"))

Table_11_prety <-  as_flextable(percent_emm_pairs_gropued) %>% 
   bold(j = 1, i = ~ !is.na(percent), bold = TRUE, part = "body" ) %>%
   align(part = "all") %>% # left align
   set_caption(caption = "Table 11. Pairwise comparisons of the treatment combinations for the first round of the field trial (only showing significant differences).") %>%
   font(fontname = "Times New Roman", part = "all") %>% 
   width(width = c(2,1,1,1,1)) %>%
   theme_booktabs()

Table_11_prety <- fontsize(Table_11_prety, part = "all", size = 8) 

Table_11_prety
```