Guam_Data_Analysis

set up

# load packages
library(readr)
library(ggplot2)
library(dplyr)
library(corrplot)
library(DataExplorer)
library(reshape2)
library(tidyverse) 
library(MASS)
library(car)
library(ggpubr)
library(sjPlot)
library(ggord)
library(ggordiplots)
library(Hmisc)
library(rstatix)
library(qqplotr)

# load data
Guam <- read_csv("Guam_Final.csv")

scale predictor variables

# Scale only the predictor variables
Guam <- Guam %>%
  mutate(across(`>1m`:Abundance, scale))

Exploratory Analysis

###### Distribution ######
# Morphology
Guam %>%
  select(Mounding, `Laminar Columnar`, `Encrusting Flat`, Branching) %>%
  plot_density()

# Genera
Guam %>%
   dplyr::select(Porites, Montipora, Dipsastrea, Leptastrea) %>%
  plot_density()

# Sizes
Guam %>%
  select(`0-5cm`:`>1m`) %>%
  plot_density()

# Community
Guam %>%
  select(RICHNESS, SHANNON,Abundance) %>%
  plot_density()

# Complexity metrics
Guam %>%
  dplyr::select(surf_complex, mean_curvature, mean_slope, mean_vrm) %>%
  plot_density()

# Environmental
Guam %>%
  dplyr::select(wave_exposure) %>%
  plot_density()

###### Normality #######
plot_qq(Guam$Mounding)
plot_qq(Guam$`Laminar Columnar`)
plot_qq(Guam$`Encrusting Flat`)
plot_qq(Guam$surf_complex)
plot_qq(Guam$mean_profile_curvature)
plot_qq(Guam$mean_vrm)
plot_qq(Guam$s11_20)
plot_qq(Guam$Porites)
plot_qq(Guam$Leptastrea)
plot_qq(Guam$Montipora)

Correlation Plots

covGuam <- Hmisc::rcorr(cor(Guam[, -c(1:3)]), type='spearman')

# Sizes with Genera - three plots for better visibility
corrplot(corr=covGuam$r[3:14, 15:25], p.mat=covGuam$P[3:14, 15:25], sig.level=0.05, method='number') 
corrplot(corr=covGuam$r[3:14, 26:36], p.mat=covGuam$P[3:14, 26:36], sig.level=0.05,  method='number') 
corrplot(corr=covGuam$r[3:14, 37:46], p.mat=covGuam$P[3:14, 37:46], sig.level=0.05, method='number') 

# Sizes with Morphology
corrplot(corr=covGuam$r[3:14,49:59], p.mat=covGuam$P[3:14,49:59], sig.level=0.05, method='number')

# Genera with Morphology - three plots for better visibility
corrplot(corr=covGuam$r[15:25, 49:59], p.mat=covGuam$P[15:25, 49:59], sig.level=0.05, method='number')
corrplot(corr=covGuam$r[26:36, 49:59], p.mat=covGuam$P[26:36, 49:59], sig.level=0.05, method='number')
corrplot(corr=covGuam$r[37:46, 49:59], p.mat=covGuam$P[37:46, 49:59], sig.level=0.05,  method='number')

# Complexity with Morphology
corrplot(corr=covGuam$r[49:59, 60:65], p.mat=covGuam$P[49:59, 60:65], sig.level=0.05, method='number')

# Complexity with Size
corrplot(corr=covGuam$r[3:14, 60:65], p.mat=covGuam$P[3:14, 60:65], sig.level=0.05, method='number')

# Complexity with Community
corrplot(corr=covGuam$r[47:48, 60:65], p.mat=covGuam$P[47:48, 60:65], sig.level=0.05, method='number')

simple linear regressions

####### vrm & community assemblage #####
vrm_diver <- lm(mean_vrm ~ SHANNON, data = Guam)
summary(vrm_diver)

summary(lm(surf_complex ~ SHANNON, data = Guam))

# plot
diverplot <- ggplot(Guam, aes(x = SHANNON, y = mean_vrm)) +
  geom_point(size=0.6) + geom_smooth(method = "lm", se = TRUE, color="#4669F6", fill = "#69B4F6") +
  theme_classic() +
  labs(tag = "A") +
  xlab("Coral Genus Diversity") + ylab("VRM") +
  theme(axis.text = element_text(size=6, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 
# complexity ~ abundance 
vrm_abund <- lm(surf_complex ~ Abundance, data = Guam)
summary(vrm_abund)
# plot
abundplot <- ggplot(Guam, aes(x = Abundance, y = surf_complex)) +
  geom_point(size=0.6) + geom_smooth(method = "lm", se = TRUE, color="#4669F6", fill = "#69B4F6") +
  theme_classic() +
  labs(tag = "B") +
  xlab("Coral Abundance") + ylab("Surface Complexity") +
  theme(axis.text = element_text(size=6, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 
# complexity ~ richness 
vrm_rich <- lm(mean_slope ~ Abundance, data = Guam)
summary(vrm_rich)
ggplot(Guam, aes(x = RICHNESS, y = mean_profile_curvature)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)
# plot
richplot <- ggplot(Guam, aes(x = RICHNESS, y = surf_complex)) +
  geom_point(size=0.6) + geom_smooth(method = "lm", se = TRUE, color="#4669F6", fill = "#69B4F6") +
  theme_classic() +
  labs(tag = "C") +
  xlab("Coral Genus Richness") + ylab("Surface Complexity") +
  theme(axis.text = element_text(size=6, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 


gridExtra::grid.arrange(diverplot, abundplot, richplot, ncol=3)

##### surface complexity & morphologies ########
summary(lm(surf_complex ~ Columnar, data = Guam))
summary(lm(surf_complex ~ `Encrusting Columnar`, data = Guam))
summary(lm(surf_complex ~ `Encrusting Flat`, data = Guam))
summary(lm(surf_complex ~ Laminar, data = Guam))
summary(lm(surf_complex ~ `Laminar Columnar`, data = Guam))
summary(lm(surf_complex ~ Mounding, data = Guam))


##### Profile Curvature & morphologies ########
summary(lm(mean_profile_curvature ~ Branching, data = Guam))
summary(lm(mean_profile_curvature ~ `Encrusting Columnar`, data = Guam))

##### Planform Curvature & morphologies ########
summary(lm(mean_plan_curvature ~ Mounding, data = Guam))
summary(lm(mean_plan_curvature ~ Columnar, data = Guam))

##### surface complexity & genera ########
summary(lm(surf_complex ~ Astreopora, data = Guam))
ggplot(Guam, aes(x = Astreopora, y = surf_complex)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)

summary(lm(surf_complex ~ Dipsastrea, data = Guam))
ggplot(Guam, aes(x = Dipsastrea, y = surf_complex)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)

summary(lm(mean_slope ~ Millepora, data = Guam))
ggplot(Guam, aes(x = Sinularia, y = surf_complex)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)


##### Profile Curvature & genera ########
summary(lm(mean_profile_curvature ~ Porites, data = Guam))


##### surface complexity & sizes ########
summary(lm(mean_plan_curvature ~ `>1m`, data = Guam))
ggplot(Guam, aes(x = `>1m`, y = surf_complex)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)

summary(lm(surf_complex ~ s31-40, data = Guam))
summary(lm(surf_complex ~ s41_50, data = Guam))
summary(lm(surf_complex ~ s51_60, data = Guam))
summary(lm(surf_complex ~ s61_70, data = Guam))
summary(lm(surf_complex ~ s1m, data = Guam))

ggplot(Guam, aes(x = `0-5cm`, y = surf_complex)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)

#####  Profile Curvature & sizes ########
summary(lm(mean_profile_curvature ~ s41_50, data = Guam))


summary(lm(SHANNON ~ Branching, data = Guam))
ggplot(Guam, aes(x = SHANNON, y = surf_complex)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)






ggplot(Guam, aes(x = Branching, y = mean_profile_curvature)) +
  geom_point(size=0.6) + geom_smooth(method = "lm", se = TRUE, color="#4669F6", fill = "#69B4F6") +
  theme_classic() +
  labs(tag = "A") +
  xlab("Branching") + ylab("mean_profile_curvature") +
  theme(axis.text = element_text(size=6, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 

ggplot(Guam, aes(x = `Encrusting Columnar`, y = mean_profile_curvature)) +
  geom_point(size=0.6) + geom_smooth(method = "lm", se = TRUE, color="#4669F6", fill = "#69B4F6") +
  theme_classic() +
  labs(tag = "A") +
  xlab("Encrsuting Columnar") + ylab("mean_profile_curvature") +
  theme(axis.text = element_text(size=6, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 

percentages of traits bases on all data

RawBenthicData %>% 
  group_by(morphology) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
    ggplot(mapping = aes(x = reorder(morphology, -Freq), Freq)) + 
  geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
    theme_classic() +
  xlab("Count") + ylab("Morphology Percentage") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold"))  +
     theme(axis.text.x = element_text(angle = 55, vjust = 0.8, hjust=0.7))

RawBenthicData %>% 
  group_by(genus) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = (n/sum(n)) * 100) %>%
  ggplot(mapping = aes(x = reorder(genus, -Freq), Freq)) + 
  geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
  theme_classic() +
  xlab("Count") + ylab("Genus Percentage") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
   theme(axis.text.x = element_text(angle = 55, vjust = 0.8, hjust=0.5))

RawBenthicData %>% 
  group_by(size) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = (n/sum(n)) * 100) %>%
    ggplot(mapping = aes(x = reorder(size,-Freq), Freq)) + 
  geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
    theme_classic() +
  xlab("Count") + ylab("Size Percentage") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
   theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.5))



RawBenthicData %>% 
  group_by(size, wave_exposure) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  ggplot(mapping = aes(x = size, y = Freq, fill = as.factor(wave_exposure))) + 
   xlab("Size Category") + ylab("Proportion") +
  theme_classic() +
  geom_bar(stat = "identity", position = "stack") +
    scale_fill_brewer(palette = "Accent") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
   theme(axis.text.x = element_text(angle = 30, vjust = 0.8, hjust=0.7)) +
   theme(legend.position="none") 

mille <- RawBenthicData %>% 
  group_by(size, genus) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  filter(genus == "Millepora") 

mille$size <- factor(mille$size, levels = c("6-10cm", "11-20cm", "21-30cm", 
              "41-50cm",  "61-70cm", "81-90cm","91-100cm"))

mille %>%
  ggplot(mapping = aes(x = size, y = n)) + 
   xlab("Size Category") + ylab("n") +
  theme_classic() +
  geom_bar(stat = "identity", position = "dodge", fill = "#CFE5B3") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) 

##### sizes among morphology #######
encflat <- RawBenthicData %>% 
  group_by(morphology, size) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  filter(morphology == "Encrusting Flat") 
# reorder
encflat$size <- factor(encflat$size, levels = c("0-5cm", "6-10cm", "11-20cm", "21-30cm", "31-40cm",
                                            "41-50cm", "51-60cm", "61-70cm", "71-80cm", "81-90cm", ">1m"))

## plot
enflat <- ggplot(encflat, mapping = aes(x = size, y =n)) + 
          geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
  labs(title = "Encrusting Flat") +
  theme_classic() +
  labs(tag = "A") +
  xlab("Size Category") + ylab("n") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.5))

##### sizes among morphology #######
RawBenthicData$morphology <- as.factor(RawBenthicData$morphology)
levels(RawBenthicData$morphology)[levels(RawBenthicData$morphology)=='Encrusting Mounding'] <- 'Mounding'
                                            
mound <- RawBenthicData %>% 
  group_by(morphology, size) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  filter(morphology == "Mounding") 
# reorder
mound$size <- factor(mound$size, levels = c("0-5cm", "6-10cm", "11-20cm", "21-30cm", "31-40cm",
                                    "41-50cm", "51-60cm", "61-70cm", "71-80cm", "81-90cm","91-100cm", ">1m"))

## plot
mounding <- ggplot(mound, mapping = aes(x = size, y =n)) + 
          geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
  labs(title = "Mounding") +
  theme_classic() +
  labs(tag = "B") +
  xlab("Size Category") + ylab("n") +
  theme(axis.text = element_text(size=7, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.5))


# laminar columnar
lamcol <- RawBenthicData %>% 
  group_by(morphology, size) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  filter(morphology == "Laminar Columnar") 
# reorder
lamcol$size <- factor(lamcol$size, levels = c("0-5cm", "6-10cm", "11-20cm", "21-30cm", "31-40cm",
              "41-50cm", "51-60cm", "61-70cm", "71-80cm", "81-90cm","91-100cm" ,">1m"))
## plot
lamcol <- ggplot(lamcol, mapping = aes(x = size, y =n)) + 
          geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
  labs(title = "Laminar Columnar") +
  theme_classic() +
  labs(tag = "C") +
  xlab("Size Category") + ylab("n") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.5))

#  columnar
col <- RawBenthicData %>% 
  group_by(morphology, size) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  filter(morphology == "Columnar") 
# reorder
col$size <- factor(col$size, levels = c("0-5cm", "6-10cm", "11-20cm", "21-30cm", "31-40cm",
              "41-50cm", "51-60cm", "61-70cm", "71-80cm", "81-90cm", ">1m"))
## plot
col <- ggplot(col, mapping = aes(x = size, y =n)) + 
          geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
  labs(title = "Columnar") +
  theme_classic() +
  labs(tag = "D") +
  xlab("Size Category") + ylab("n") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.5))

#  laminar
lam <- RawBenthicData %>% 
  group_by(morphology, size) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  filter(morphology == "Laminar") 
# reorder
lam$size <- factor(lam$size, levels = c("0-5cm", "6-10cm", "11-20cm", "21-30cm", "31-40cm",
              "41-50cm", "51-60cm", "61-70cm", "71-80cm", "81-90cm","91-100cm", ">1m"))
## plot
lam <- ggplot(lam, mapping = aes(x = size, y =n)) + 
          geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
  labs(title = "Laminar") +
  theme_classic() +
  labs(tag = "E") +
  xlab("Size Category") + ylab("n") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.5))

#  branching
bran <- RawBenthicData %>% 
  group_by(morphology, size) %>%
  dplyr::summarise(n = n()) %>%
  mutate(Freq = n/sum(n)) %>%
  filter(morphology == "Branching") 
# reorder
bran$size <- factor(bran$size, levels = c("0-5cm", "6-10cm", "11-20cm", "21-30cm", "31-40cm",
              "41-50cm", "51-60cm", "61-70cm", "71-80cm", ">1m"))
## plot
bran <- ggplot(bran, mapping = aes(x = size, y =n)) + 
          geom_bar(stat = "identity", position = "dodge", fill = "#8BBFFA") +
  theme(legend.position="none") +
  labs(title = "Branching") +
  theme_classic() +
  labs(tag = "F") +
  xlab("Size Category") + ylab("n") +
  theme(axis.text = element_text(size=9, color="black")) +
  theme(axis.title.x = element_text(size = 10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.5))



gridExtra::grid.arrange(enflat, mounding,lamcol, col, lam, bran, ncol=2)

multiple linear regressions

########### Surf Complex VS SIZE CATEGORIES #############
# Fit the full model 
surf_sizes <- lm(mean_slope ~., data=Guam[c(64,2:13)])   
# Stepwise regression model
surf_sizes_step <- stepAIC(surf_sizes, direction = "both", trace = FALSE)
# Results
summary(surf_sizes_step)
plot(surf_sizes_step)
avPlots(surf_sizes_step)
vif_values <- vif(surf_sizes_step)
print(vif_values)
barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue")
abline(v = 5, lwd = 2, lty = 2)
C <- plot_model(surf_sizes,  sort.est = TRUE, show.values = TRUE, vline.color = "grey",
    value.size = 2.5, dot.size = 1.5, line.size = 0.8, value.offset = .35, title = " ") + 
   ylim(-.2, .4) +
   theme_bw() +
   labs(tag="A") +
  theme(axis.text = element_text(size=9.5, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 


bla <- lm(surf_complex ~., data=Guam[c(63,2:13)])

##### plot #####
Guam %>%
  subset(select=c(`>1m`, `0-5cm`, `11-20cm`, `41-50cm`, `51-60cm`, surf_complex)) %>%
  pivot_longer(!surf_complex, names_to = "Trait", values_to = "Proportion") %>%
  ggplot(aes(Proportion, surf_complex, fill = Trait)) +
  geom_smooth(method="lm", size=1.2, aes(color = Trait)) +
  #scale_fill_brewer(palette="Set2") +
 # scale_colour_brewer(palette = "Set2") +
  theme(legend.position = "right", legend.text = element_text(size=11), 
         legend.title =   element_text(size=13.5, face = "bold"),
          axis.text = element_text(size=12, color="black"), legend.key.size = unit(9,"mm")) +
  theme(axis.title.x = element_text(size = 13.5, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =13.5, color = "black", face = "bold")) +
  bgcolor("white") +
  labs(title = "Multiple Linear Regression for Surface Complexity") + 
  theme(plot.title = element_text(size=13.5,face = "bold", color = "black", hjust = 0.5)) +
  xlab("Traits (scaled)") + ylab("Surface Complexity") +
  theme(panel.background = element_rect(fill=NA, colour = "black", size=1, linetype = "solid")) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
         panel.grid.minor = element_blank()) 

## once colonies get passed 40cm, there's a positive relationship between surfcomplex & size

########### Surf Complex VS MORPHOLOGY CATEGORIES #############
# Fit the full model 
surf_morph <- lm(surf_complex ~., data=Guam[c(63,46:56)])   
# Stepwise regression model
surf_morph_step <- stepAIC(surf_morph, direction = "both", trace = FALSE)
# Results
summary(surf_morph_step)
plot(surf_morph_step)
avPlots(surf_morph_step)
vif_values <- vif(surf_morph_step)
print(vif_values)
barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue")
abline(v = 5, lwd = 2, lty = 2)
B <- plot_model(surf_morph,  sort.est = TRUE, show.values = TRUE, vline.color = "grey",
    value.size = 2.5, dot.size = 1.5, line.size = 0.8, value.offset = .35, title = " ") + 
   ylim(-.25,.5) +
   theme_bw() +
   labs(tag="B") +
  theme(axis.text = element_text(size=9.5, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 

summary(lm(mean_plan_curvature ~., data=Guam[c(68,14:45)]))  

########### Surf Complex VS GENERA CATEGORIES #############
# Fit the full model 
surf_gen <- lm(surf_complex ~., data=Guam[c(63,14:45)])   
# Stepwise regression model
surf_gen_step <- stepAIC(surf_gen, direction = "both", trace = FALSE)
# Results
summary(surf_gen_step)
plot(surf_gen_step)
avPlots(surf_gen_step)
vif_values <- vif(surf_gen_step)
print(vif_values)
barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue")
abline(v = 5, lwd = 2, lty = 2)
C <- plot_model(surf_gen,  sort.est = TRUE, show.values = TRUE, vline.color = "grey",
    value.size = 2.5, dot.size = 1.5, line.size = 0.8, value.offset = .35, title = " ") + 
   ylim(-.3, .2) +
   theme_bw() +
   labs(tag="C") +
  theme(axis.text = element_text(size=9.5, color="black")) +
  theme(axis.title.x = element_text(size = 9, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =9, color = "black", face = "bold")) 


###### sizes ######
A <- lm(surf_complex ~., data=Guam[c(63,2:13)])   
B <- lm(mean_slope ~., data=Guam[c(64,2:13)])   
C <- lm(mean_vrm ~., data=Guam[c(65,2:13)])
D <- lm(mean_profile_curvature ~., data=Guam[c(67,2:13)])   
FF <- lm(mean_plan_curvature ~., data=Guam[c(68,2:13)])   

sjPlot::plot_models( A, B, C, D, FF,  dot.size = 3.9, line.size = 0.9,
   legend.title = "Complextiy Metric", pacing = 0.6, p.shape = TRUE, vline.color = "grey") +
   ylim(-3.7, 6.1) +
   theme_classic() +
   theme(axis.text = element_text(size=16, color="black")) +
   labs(tag="A") 

###### morphologies ######
A <- lm(surf_complex ~., data=Guam[c(63,46:56)])  
B <- lm(mean_slope ~., data=Guam[c(64,46:56)])   
C <- lm(mean_vrm ~., data=Guam[c(65,46:56)])
D <- lm(mean_profile_curvature ~., data=Guam[c(67,46:56)])   
FF <- lm(mean_plan_curvature ~., data=Guam[c(68,46:56)])   

sjPlot::plot_models( A, B, C, D, FF,  dot.size = 3.9, line.size = 0.9,
    legend.title = "Complextiy Metric", spacing = 0.6, p.shape = TRUE, vline.color = "grey") +
  ylim(-3.7, 3) +
  theme_classic() +
  theme(axis.text = element_text(size=16, color="black")) +
     labs(tag="B") 

###### morphologies ######
A <- lm(surf_complex ~., data=Guam[c(63,14:45)])  
A <- stepAIC(A, direction = "both", trace = FALSE)
B <- lm(mean_slope ~., data=Guam[c(64,14:45)]) 
B <- stepAIC(B, direction = "both", trace = FALSE)
C <- lm(mean_vrm ~., data=Guam[c(65,14:45)])
C <- stepAIC(C, direction = "both", trace = FALSE)
D <- lm(mean_profile_curvature ~., data=Guam[c(67,14:45)])   
D <- stepAIC(D, direction = "both", trace = FALSE)
FF <- lm(mean_plan_curvature ~., data=Guam[c(68,14:45)])  
FF <- stepAIC(FF, direction = "both", trace = FALSE)

sjPlot::plot_models( A, B, C, D, FF,  dot.size = 3.9, line.size = 0.9,
   legend.title = "Complextiy Metric", spacing = 0.6, p.shape = TRUE, vline.color = "grey") +
  ylim(-4.5, 2.5) +
   theme_classic() +
  theme(axis.text = element_text(size=16, color="black")) +
   labs(tag="C") 

community assemblage regressions

###### VRM - taxonomy #####
# model
vrm_tax <- lm(mean_vrm ~  RICHNESS + Abundance + SHANNON, data = Guam)
# results
summary(vrm_tax)
plot(vrm_tax)
avPlots(vrm_tax)
vif_values <- vif(vrm_tax)
print(vrm_tax)
barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue")
abline(v = 5, lwd = 2, lty = 2)

# plot estimates #
plot_model(vrm_tax,  sort.est = TRUE, show.values = TRUE, vline.color = "grey",
    value.size = 4, dot.size = 3, line.size = 1.5, value.offset = .3) + 
   ylim(-.25, .25) +
   theme_bw(base_size=13) +
   labs(title = "Estimates from VRM Regression Model") 

RDA

# subset traits
Traits <- Guam %>%
    subset(select=c(plot, `>1m`:`Mounding Lobate`)) %>%
    column_to_rownames(var = "plot")

# subset structure
Structure <-  Guam %>%
   subset(select=c(plot, surf_complex:mean_plan_curvature)) %>%
   subset(select=-c(mean_curvature)) %>%
  column_to_rownames(var = "plot")

# Run Structure ~ Trait RDA 
rda1 <- vegan::rda(Structure ~ ., data = Traits, scale=TRUE)
# Results
anova(rda1) # is the model significant?
anova.cca(rda1) # is the model significant?
head (summary (rda1))
anova(rda1, permutations = 1000, by="terms") 
anova(rda1, by="axis") 
RsquareAdj(rda1) # adjusted Rsquared 
sqrt(vif.cca(rda1))
vif.cca(rda1) # variance inflation factor (<10 OK)
ggord(rda1, obslab = FALSE, ptslab = TRUE)
ordiplot(rda1)
ordiplot(rda1,type="t",scaling=2)
ordiplot(rda1,type="t",scaling=1)

####### plot #######
#  add the genera variables as arrows
arrowmatob = vegan::scores(rda1, display = "sp")
# Add labels, make a data.frame
arrowdfob <- data.frame(labels = rownames(arrowmatob), arrowmatob)
# Define the arrow aesthetic mapping
arrowmatob = aes(xend = 1.1 * RDA1, yend =  1.1 * RDA2, x = 0, y = 0, shape = NULL,
          color ="#369592", label = labels)
label_map = aes(x = 1.2 * RDA1, y = 1.12 * RDA2, shape = NULL,
          color = "#369592", label = labels)
# Make a new graphic
arrowheadob = arrow(length = unit(2, "mm"))
# plot
ggord(rda1, grp_in = Guam$site, axes = c("1", "2"), ellipse= FALSE, 
  arrow = 0.1, txt = 3, color = "#369592",
  repel = TRUE, size=2, xlim=c(-0.85,0.8), ylim=c(-0.45,0.75)) +
  theme(legend.position = "none") +
  theme(axis.title.x = element_text(size = 12, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size=12, color = "black", face = "bold")) +
  theme(plot.title = element_text(size=12,face = "bold", color = "black", hjust = 0.5)) +
  theme(panel.background = element_rect(fill=NA, colour = "black", size=1, linetype = "solid")) +
  labs(title = "RDA -- Structure ~ Traits") + 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
  panel.grid.minor = element_blank())  +
  guides(colour = guide_legend(override.aes = list(size=5, stroke=1))) 

tyring size things

## try first to identify the mode for each plot, so if the mode was 31-40, then that plot receives the dominant size value of 35.  This probably won't work if the size class distribution is not unimodal though...  

# calculate number of colonies per size category for each plot 
sizes_mode <- RawBenthicData %>% group_by(plot,size) %>%
  tally() %>%
  spread(size, n) %>%
  replace(is.na(.), 0)

sizes_mode$Max <- apply(sizes_mode[,-1], 1, function(x) paste0(names(sizes_mode[,-1])[x == max(x)],collapse = '_'))

## We may also be able to quantify the number of coral size classes per plot, and that gives us size class variability or even calculate diversity indices using size class data.
sizes_plot <- RawBenthicData %>% group_by(plot,size) %>%
  tally() %>%
  spread(size, n) %>%
  replace(is.na(.), 0)

size_diversity <- ddply(sizes_plot,~plot,function(x) {
  data.frame(size_shannon=diversity(x[-1], index="shannon")) })

## complexity vs size diversity
sizestructure <- size_diversity %>% 
  inner_join(ComplexityData, by="plot") 

ggplot(sizestructure, aes(x = size_shannon, y = surf_complex)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

environmental variables

######## substrate type #########
Guam$habitat_type <- as.factor(Guam$habitat_type)
# distribution per level
ggplot(data=Guam,aes(x = surf_complex, fill = habitat_type)) +
  geom_density(alpha=0.15)  +
  facet_wrap(habitat_type ~ .)+theme_bw()
# normality test - rock & boulder not normal
Guam %>%
  group_by(habitat_type) %>%
  shapiro_test(surf_complex)
# try log transformation
Guam <- Guam %>% 
  mutate(log_surf = log(surf_complex))
# check if normal
Guam %>%
   group_by(habitat_type) %>%
  shapiro_test(log_surf)
# ANOVA test
surfhabaov <- aov(log_surf ~ habitat_type, data = Guam)
summary(surfhabaov)
# get mean & sd
surfhab <- Guam %>%
  group_by(site, habitat_type) %>%
  dplyr::summarise(sd = sd(log_surf), len = mean(log_surf))
# surf complex ~ habitat
habitat <- ggplot(surfhab, mapping= aes(x = reorder(habitat_type, len), y = len), size=1) +
  geom_boxplot(fill ='#F5BABA') +
  theme_classic() +
  theme(legend.position = "none") +
  labs(tag="A") + 
   xlab("Habitat Type") + ylab("Log Surface Complexity") +
  theme(axis.text = element_text(size=10, color="black")) +
  theme(axis.title.x = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold")) +
  scale_x_discrete(labels=c("Pavement" = "Pavement", "Aggregate Reef" = "Aggregate", "Rock/Boulder" = "Rock/Boulder"))


######## wave exposure #########
Guam$wave_exposure <- as.factor(Guam$wave_exposure)
# distribution per level
ggplot(data=Guam,aes(x = surf_complex, fill = wave_exposure)) +
  geom_density(alpha=0.15)  +
  facet_wrap(wave_exposure ~ .)+theme_bw()
# normality test -- not normal with raw data so using the log-transformed data
Guam %>%
  group_by(wave_exposure) %>%
  shapiro_test(surf_complex)
Guam %>%
  group_by(wave_exposure) %>%
  shapiro_test(log_surf)
# ANOVA test
surfwaveaov <- aov(log_surf ~ wave_exposure, data = Guam)
summary(surfwaveaov)
# get mean & sd
surfwave <- Guam %>%
  group_by(site, wave_exposure) %>%
  dplyr::summarise(sd = sd(log_surf), len = mean(log_surf))
# wave exposure ~ habitat
wave <- ggplot(surfwave, mapping= aes(x = reorder(wave_exposure, len), y = len), size=1) +
   geom_boxplot(fill ='#F5BABA') +
  theme_classic() +
  theme(legend.position = "none") +
  labs(tag="B") + 
  xlab("Wave Exposure") + ylab("Log Surface Complexity") +
  scale_x_discrete(labels=c("3" = "High", "2" = "Moderate", "1" = "Low")) +
  theme(axis.text = element_text(size=10, color="black")) +
  theme(axis.title.x = element_text(size =10, color = "black", face = "bold")) +
  theme(axis.title.y = element_text(size =10, color = "black", face = "bold"))



#### combined ######
gridExtra::grid.arrange(habitat, wave, ncol=2)