setwd("C:/Users/adikh/Desktop")
data=read.csv(file="Bleaching literature review subset.csv", header=TRUE, sep=",")

#subset by region
Atlantic=subset(data, subset = Geographic.region %in% "Atlantic")
Indian=subset(data, subset = Geographic.region %in% "Indian")
Pacific=subset(data, subset = Geographic.region %in% "Pacific")
Combo=rbind(Atlantic, Indian, Pacific)

#subset genera present in most or all regions
Genus.df <- subset(Combo, Genus %in% c("Acropora", "Astreopora", "Cyphastrea", "Dipsastraea", "Echinopora", "Favia", 
                                       "Favites", "Galaxea", "Goniastrea", "Goniopora", "Hydnophora","Lobophyllia", 
                                       "Montastraea", "Montipora", "Pavona", "Platygyra", "Pocillopora", 
                                       "Porites", "Stylophora", "Turbinaria"))

Genus.df$Bleaching.severity.score=as.numeric(Genus.df$Bleaching.severity.score)
Genus.df2=na.omit(Genus.df)

library(ggplot2)
library(EnvStats)
library(ggpubr)

#find observation sample sizes
give.n <- function(x){
  return(c(y = mean(x), label = length(x)))
}


library(RColorBrewer)
mycolors = c(brewer.pal(name="Set2", n = 3))
symnum.args <- list(cutpoints = c(0, 0.01, Inf), symbols = c("*", "ns"))

ggplot(Genus.df2, aes(x=Genus, y=Bleaching.severity.score, fill=Geographic.region)) + 
  theme_classic(base_size=14) + 
  #geom_boxplot(width=0.5, position = position_dodge(preserve = "single"))+
  geom_violin(trim=T, width=0.6, scale="width", position = position_dodge(width = 0.5))+
  stat_summary(aes(y = 4.1),geom = "text", fun.data = give.n, position = position_dodge(width = 0.75), hjust=0, size=3, angle=90)+ 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  ylab("Bleaching severity score")  + guides(fill=guide_legend(title="Region"))+
  scale_fill_manual(values = mycolors)+
  stat_compare_means(method="anova",label="p.signif",label.y = 4.5, hide.ns=T, symnum.args=symnum.args)


#subset morphologies present in most or all regions
Morph.df <- subset(Combo, Morphology %in% c("branching", "columnar", "encrusting", 
                                            "free-living", "massive", "massive / encrusting",
                                            "plating", "table / branching"))

Morph.df$Bleaching.severity.score=as.numeric(Morph.df$Bleaching.severity.score)
Morph.df2=na.omit(Morph.df)

ggplot(Morph.df2, aes(x=Morphology, y=Bleaching.severity.score, fill=Geographic.region)) + 
  theme_classic(base_size=14) + 
  #geom_boxplot(width=0.5, position = position_dodge(preserve = "single"))+
  geom_violin(trim=T, width=0.6, scale="width", position = position_dodge(width = 0.5))+
  stat_summary(aes(y = 4.1),geom = "text", fun.data = give.n, position = position_dodge(width = 0.75), hjust=0, size=3, angle=90)+ 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  ylab("Bleaching severity score")  + guides(fill=guide_legend(title="Region"))+
  scale_fill_manual(values = mycolors) +
  stat_compare_means(method="anova",label="p.signif",label.y = 4.5, hide.ns=T, symnum.args=symnum.args)

#two separate two-way ANOVAs because not all genera have all morphologies
#aligned rank transform ANOVAs because data are non-continuous
#install.packages("ARTool")

library(ARTool)
anova1=art(Bleaching.severity.score~factor(Geographic.region)*factor(Genus), Genus.df2)
anova(anova1, type=2)

anova2=art(Bleaching.severity.score~factor(Geographic.region)*factor(Morphology), Morph.df2)
anova(anova2, type=2)