###
# Figure 2 - Black Coral Density Comparisons
# Power analysis of A. Caribbeana density
###

rm(list=ls())

library('dplyr')

df1 <- read.csv("Data S2.csv")
covariates <- read.csv('Data S1.csv')

# setting variables as right type
df1$Site <- as.character(df1$Site)

# selecting just data with survey areas
df2 <- df1[which(is.na(df1$Transect.Area.m2)==FALSE), ]

# expanding for transects with 0 black corals recorded on during 2016
sites <- c("CO", "HE", "PJ", "PT", "PU", "SR", "VB", "TT")

df3 <- expand.grid(Year=2016, Site=sites, Transect=1:4, Area=NA, caribbeana=0, pennacea=0)

temp <- data.frame(Year=rep(1998, length(unique(df2$Site[which(df2$Year==1998)]))),
                   Site=unique(df2$Site[which(df2$Year==1998)]),
                   Transect=rep(1, length(unique(df2$Site[which(df2$Year==1998)]))),
                   Area=NA,
                   caribbeana=0,
                   pennacea=0)

####
## Removing San Francisco - as there's no black coral ID for these
####

temp <- temp[which(temp$Site!="San Francisco"), ]

df3 <- rbind(df3, temp)

# adding in the 2016 data
for(i in 1:nrow(df3))
{
  temp <- df2[which(df2$Year==df3$Year[i] &
                      df2$Site==df3$Site[i] &
                      df2$Transect==df3$Transect[i]), ]
  
  if (nrow(temp) > 0)
  {
    df3[i, "Area"] <- temp$Transect.Area.m2[1]
  } else {
    df3[i, "Area"] <- 120
  }
  df3[i, "caribbeana"] <- length(which(temp$Species=='caribbeana'))
  df3[i, "pennacea"] <- length(which(temp$Species=='pennacea'))
  
}


# Mean density per transect & changing mean density into per 100 m^2 so that numbers are not all super small


df3$caribbeana <- df3$caribbeana/df3$Area * 100
df3$pennacea <- df3$pennacea/df3$Area * 100

# Mean density per site and year
caribbeana <- aggregate(caribbeana ~ Year + Site, data=df3, FUN=mean)
pennacea <- aggregate(pennacea ~ Year + Site, data=df3, FUN=mean)

# Adding in contextual variables

caribbeana <- left_join(caribbeana, covariates, by='Site')
pennacea <- left_join(pennacea, covariates, by='Site')



###
# Differences in density between years
###

library('car')

# se function
se <- function(x) {sd(x)/ sqrt(length(x))}


# Caribbeana summary stats (for main text)
car.density <- aggregate(caribbeana ~ Year.x, data=caribbeana, FUN=mean)
car.sd <- aggregate(caribbeana ~ Year.x, data=caribbeana, FUN=sd)
car.se <- aggregate(caribbeana ~ Year.x, data=caribbeana, FUN=se)



# Pennacea summary stats (for main text)
pen.density <- aggregate(pennacea ~ Year.x, data=pennacea, FUN=mean)
pen.sd <- aggregate(pennacea ~ Year.x, data=pennacea, FUN=sd)
pen.se <- aggregate(pennacea ~ Year.x, data=pennacea, FUN=se)


##
# Power analysis of A. Caribbeana density
##

diff.in.mean <- car.density$caribbeana[1] - car.density$caribbeana[2]
pooled.sd <- sqrt((car.sd$caribbeana[1]^2 + car.sd$caribbeana[2]^2)/2)

pwr.t2n.test(d=(diff.in.mean/pooled.sd),
             n1=table(caribbeana$Year.x)[1], n2=table(caribbeana$Year.x)[2],
             sig.level=.05, alternative="two.sided")



###
# Plotting Figure 2
###

png("Fig 2.png", width = 4, height = 5, units = 'in', res = 300)


par(mfrow=c(1, 2))
par(oma=c(1, 1, 0, 0))

location <- c(0.7, 1.9)


# Caribbeana
barplot(car.density$caribbeana, axes=FALSE, ylim=c(0, 1))
axis(1, labels=c("", 1998, 2016, ""), at=c(-10, location , 100), las=2)
axis(2)

arrows(location, car.density$caribbeana, location, car.density$caribbeana+car.se$caribbeana, code=2, angle=90, length=0.1)
arrows(location, car.density$caribbeana, location, car.density$caribbeana-car.se$caribbeana, code=2, angle=90, length=0.1)

legend("topleft", "A", bty='n', adj=c(2, -0.5))

# Pennacea

barplot(pen.density$pennacea, axes=FALSE, ylim=c(0, 5))
axis(1, labels=c("", 1998, 2016, ""), at=c(-10, location , 100), las=2)
axis(2)

arrows(location, pen.density$pennacea, location, pen.density$pennacea+pen.se$pennacea, code=2, angle=90, length=0.1)
arrows(location, pen.density$pennacea, location, pen.density$pennacea-pen.se$pennacea, code=2, angle=90, length=0.1)

legend("topleft", "B", bty='n', adj=c(2, -0.5))

mtext("Year", side=1, outer=TRUE, line=-0.5)
mtext(expression("Colony Density (100m"^{2} ~")"), side=2, outer=TRUE, line=-0.5)

dev.off()




###
# Figure 3 - Black Coral height and width all depths
###

rm(list=ls())


library(sm)
library(KernSmooth)


df1 <- read.csv("Data S2.csv")
df1$Site <- as.character(df1$Site)



###
# Length comparisons - all data
###


###
# Caribbeana
###

df4 <- df1[which(df1$Species=="caribbeana"), ]

# selecting based on year
y1998 <- df4[which(df4$Year==1998), ]
y2016 <- df4[which(df4$Year==2016), ]

## Height

# selection of bandwidths
h1998.width <- dpik(y1998$Height.cm)
h2016.width <- dpik(y2016$Height.cm)

# generating the length distribution outlines
ch1998.outline <- bkde(y1998$Height.cm, bandwidth=h1998.width)
ch2016.outline <- bkde(y2016$Height.cm, bandwidth=h2016.width)

# concatinating all lengths
lengths <- c(y1998$Height.cm, y2016$Height.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

caribbeana.height <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                        nboot=9999, ngrid=500)

caribbeana.height.n <- length(lengths)

## Width

# selection of bandwidths
w1998.width <- dpik(y1998$Width.cm)
w2016.width <- dpik(y2016$Width.cm)

# generating the length distribution outlines
cw1998.outline <- bkde(y1998$Width.cm, bandwidth=w1998.width)
cw2016.outline <- bkde(y2016$Width.cm, bandwidth=w2016.width)

# concatinating all lengths
lengths <- c(y1998$Width.cm, y2016$Width.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

caribbeana.width <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                       nboot=9999, ngrid=500)

caribbeana.width.n <- length(lengths)

###
# Pennacea
###

df4 <- df1[which(df1$Species=="pennacea"), ]
df4 <- df4[(is.na(df4$Height.cm)==FALSE), ]

# selecting based on year
y1998 <- df4[which(df4$Year==1998), ]
y2016 <- df4[which(df4$Year==2016), ]

## Height

# selection of bandwidths
h1998.width <- dpik(y1998$Height.cm)
h2016.width <- dpik(y2016$Height.cm)

# generating the length distribution outlines
ph1998.outline <- bkde(y1998$Height.cm, bandwidth=h1998.width)
ph2016.outline <- bkde(y2016$Height.cm, bandwidth=h2016.width)

# concatinating all lengths
lengths <- c(y1998$Height.cm, y2016$Height.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

pennacea.height <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                      nboot=99, ngrid=500)

pennacea.height.n <- length(lengths)

## Width

# selection of bandwidths
w1998.width <- dpik(y1998$Width.cm)
w2016.width <- dpik(y2016$Width.cm)

# generating the length distribution outlines
pw1998.outline <- bkde(y1998$Width.cm, bandwidth=w1998.width)
pw2016.outline <- bkde(y2016$Width.cm, bandwidth=w2016.width)

# concatinating all lengths
lengths <- c(y1998$Width.cm, y2016$Width.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

pennacea.width <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                     nboot=9999, ngrid=500)

pennacea.width.n <- length(lengths)



###
# Plotting
###


png("Fig 3 - all depths size.png", width = 6, height = 6, units = 'in', res = 300)


par(mfrow=c(2, 2))
par(oma=c(3, 3, 0, 0))
par(mar=c(1, 1, 1, 1))

##
# Plotting Caribbeana height
##


plot(ch1998.outline, lwd=2, type='n', xlim=c(0, 250), ylim=c(0, 0.025),
     axes=FALSE)

axis(1, labels=FALSE)
axis(2, labels=TRUE)

polygon(c(caribbeana.height$eval.points, rev(caribbeana.height$eval.points)), c(caribbeana.height$upper, rev(caribbeana.height$lower)),
        col = "grey", border = NA)

lines(ch1998.outline, lwd=2, lty=1)
lines(ch2016.outline, lwd=2, lty=2)

text(210, 0.022, eval(paste("n = ", caribbeana.height.n)))
text(210, 0.020, eval(paste("p = ", round(caribbeana.height$p, 2))))
text(0, 0.024, substitute(paste("(A)  ", italic('A. caribbeana'), " height", sep=" ")), pos=4)


##
# Plotting pennacea Height
##


plot(ph1998.outline, lwd=2, type='n', xlim=c(0,250), ylim=c(0, 0.025),
     axes=FALSE)


axis(1, labels=FALSE)
axis(2, labels=FALSE)

polygon(c(pennacea.height$eval.points, rev(pennacea.height$eval.points)), c(pennacea.height$upper, rev(pennacea.height$lower)),
        col = "grey", border = NA)

lines(ph1998.outline, lwd=2, lty=1)
lines(ph2016.outline, lwd=2, lty=2)

text(210, 0.022, eval(paste("n = ", pennacea.height.n)))
text(210, 0.020, eval(paste("p < 0.001 ")))
text(0, 0.024, substitute(paste("(B)  ", italic('P. pennacea'), " height", sep=" ")), pos=4)

legend(150, 0.018, legend=c("1998", "2016"), lty=c(1, 2), lwd=2)



####
# Plotting Caribbeana width
####


plot(cw1998.outline, lwd=2, type='n', xlim=c(0, 250), ylim=c(0, 0.02),
     axes=FALSE)


axis(1, labels=TRUE)
axis(2, labels=TRUE)

polygon(c(caribbeana.width$eval.points, rev(caribbeana.width$eval.points)), c(caribbeana.width$upper, rev(caribbeana.width$lower)),
        col = "grey", border = NA)

lines(cw1998.outline, lwd=2, lty=1)
lines(cw2016.outline, lwd=2, lty=2)

text(210, 0.017, eval(paste("n = ", caribbeana.width.n)))
text(210, 0.0155, eval(paste("p = ", round(caribbeana.width$p, 2))))
text(0, 0.0185, substitute(paste("(C)  ", italic('A. caribbeana'), " width", sep=" ")), pos=4)




##
# Plotting pennacea width
##


plot(pw1998.outline, lwd=2, type='n', xlim=c(0, 250), ylim=c(0, 0.02),
     axes=FALSE)


axis(1, labels=TRUE)
axis(2, labels=FALSE)

polygon(c(pennacea.width$eval.points, rev(pennacea.width$eval.points)), c(pennacea.width$upper, rev(pennacea.width$lower)),
        col = "grey", border = NA)

lines(pw1998.outline, lwd=2, lty=1)
lines(pw2016.outline, lwd=2, lty=2)

text(210, 0.017, eval(paste("n = ", pennacea.width.n)))
text(210, 0.0155, eval(paste("p < 0.001 ")))
text(0, 0.0185, substitute(paste("(D)  ", italic('P. pennacea'), " width", sep=" ")), pos=4)



# Axis labels

mtext("Black Coral size (cm)", 1, line=1.5, outer=TRUE)
mtext("Probability Density", 2, line=1.5, outer=TRUE)


dev.off()




###
# Figure 4 - Black Coral Density Comparison by colony size
###

rm(list=ls())

df <- read.csv("Data S2.csv")
df$Site <- as.character(df$Site)

black.coral.density<- function(df1=data)
{
  # selecting just data with survey areas
  df2 <- df1[which(is.na(df1$Transect.Area.m2)==FALSE), ]
  
  # expanding for transects with 0 black corals recorded on during 2016
  sites <- c("CO", "HE", "PJ", "PT", "PU", "SR", "VB", "TT")
  
  df3 <- expand.grid(Year=2016, Site=sites, Transect=1:4, Area=NA, caribbeana=0, pennacea=0)
  
  temp <- data.frame(Year=rep(1998, length(unique(df2$Site[which(df2$Year==1998)]))),
                     Site=unique(df2$Site[which(df2$Year==1998)]),
                     Transect=rep(1, length(unique(df2$Site[which(df2$Year==1998)]))),
                     Area=NA,
                     caribbeana=0,
                     pennacea=0)
  
  ####
  ## Removing San Francisco - as there's no black coral ID for these
  ####
  
  temp <- temp[which(temp$Site!="San Francisco"), ]
  
  df3 <- rbind(df3, temp)
  
  # adding in the 2016 data
  for(i in 1:nrow(df3))
  {
    temp <- df2[which(df2$Year==df3$Year[i] &
                        df2$Site==df3$Site[i] &
                        df2$Transect==df3$Transect[i]), ]
    
    if (nrow(temp) > 0)
    {
      df3[i, "Area"] <- temp$Transect.Area.m2[1]
    } else {
      df3[i, "Area"] <- 120
    }
    df3[i, "caribbeana"] <- length(which(temp$Species=='caribbeana'))
    df3[i, "pennacea"] <- length(which(temp$Species=='pennacea'))
    
  }
  
  
  # Mean density per transect & changing mean density into per 100 m^2 so that numbers are not all super small
  
  
  df3$caribbeana <- df3$caribbeana/df3$Area * 100
  df3$pennacea <- df3$pennacea/df3$Area * 100
  
  # Mean density per site and year
  caribbeana <- aggregate(caribbeana ~ Year + Site, data=df3, FUN=mean)
  pennacea <- aggregate(pennacea ~ Year + Site, data=df3, FUN=mean)
  
  return(pennacea)
} # closes the black coral density function 


small <- df[which(df$Height.cm <75), ]
small <- black.coral.density(df1=small)

medium <- df[which(df$Height.cm >=75 & df$Height.cm < 150), ]
medium <- black.coral.density(medium)

large <- df[which(df$Height.cm <= 150), ]
large <- black.coral.density(large)

###
# Differences in density between years for different groups
###

# se function
se <- function(x) {sd(x)/ sqrt(length(x))}

mod1 <- wilcox.test(pennacea ~ Year, data=small)
mod1

small.den <- aggregate(pennacea ~ Year, data=small, FUN=mean)
small.se <- aggregate(pennacea ~ Year, data=small, FUN=se)


mod2 <- wilcox.test(pennacea ~ Year, data=medium)
mod2

medium.den <- aggregate(pennacea ~ Year, data=medium, FUN=mean)
medium.se <- aggregate(pennacea ~ Year, data=medium, FUN=se)

mod3 <- wilcox.test(pennacea ~ Year, data=large)
mod3

large.den <- aggregate(pennacea ~ Year, data=large, FUN=mean)
large.se <- aggregate(pennacea ~ Year, data=large, FUN=se)


# Aggregating data

density <- rbind(small.den, medium.den, large.den)
se.density <- rbind(small.se, medium.se, large.se)
density$se <- se.density$pennacea
density$size <- c("S", "S", "M", "M", "L", "L")

###
# Plotting density
###

png("Fig 4 - Density by colony size.png", width = 5, height = 4, units = 'in', res = 300)


par(mfrow=c(1, 1))

location <- c(1.3, 3.7, 6.1)
location2 <- seq(0.7, 7.2, 1.2)

# Pennacea

colours <- gray.colors(2, start = 0.5, end = 0.8)

barplot(density$pennacea, axes=FALSE, ylim=c(0, 5), col=colours)
axis(1, labels=c("", '<75', '75 - 150', '>150', ""), at=c(-10, location , 100), las=1)
axis(2)

arrows(location2, density$pennacea, location2, density$pennacea + density$se, code=2, angle=90, length=0.1)
arrows(location2, density$pennacea, location2, density$pennacea - density$se, code=2, angle=90, length=0.1)

legend("topleft", c("1998", "2016"), fill=colours, title='Year')

# adding significance stars
points(location, density$pennacea[which(density$Year=='1998')] + 0.8,
       pch="*", cex=2)

mtext("Colony Height (cm)", side=1, outer=TRUE, line=-1.5)
mtext(expression("Colony Density (100m"^{2} ~")"), side=2, outer=TRUE, line=-1.5)

dev.off()



###
# Table 1 - Black Coral Density Comparisons
###

rm(list=ls())

library('dplyr')
require(broom)

df1 <- read.csv("Data S2.csv")
covariates <- read.csv('Data S1.csv')

# setting variables as right type
df1$Site <- as.character(df1$Site)



# selecting just data with survey areas
df2 <- df1[which(is.na(df1$Transect.Area.m2)==FALSE), ]

# expanding for transects with 0 black corals recorded on during 2016
sites <- c("CO", "HE", "PJ", "PT", "PU", "SR", "VB", "TT")

df3 <- expand.grid(Year=2016, Site=sites, Transect=1:4, Area=NA, caribbeana=0, pennacea=0)

temp <- data.frame(Year=rep(1998, length(unique(df2$Site[which(df2$Year==1998)]))),
                   Site=unique(df2$Site[which(df2$Year==1998)]),
                   Transect=rep(1, length(unique(df2$Site[which(df2$Year==1998)]))),
                   Area=NA,
                   caribbeana=0,
                   pennacea=0)

####
## Removing San Francisco - as there's no black coral ID for these
####

temp <- temp[which(temp$Site!="San Francisco"), ]

df3 <- rbind(df3, temp)

# adding in the 2016 data
for(i in 1:nrow(df3))
{
  temp <- df2[which(df2$Year==df3$Year[i] &
                      df2$Site==df3$Site[i] &
                      df2$Transect==df3$Transect[i]), ]
  
  if (nrow(temp) > 0)
  {
    df3[i, "Area"] <- temp$Transect.Area.m2[1]
  } else {
    df3[i, "Area"] <- 120
  }
  df3[i, "caribbeana"] <- length(which(temp$Species=='caribbeana'))
  df3[i, "pennacea"] <- length(which(temp$Species=='pennacea'))
  
}


# Mean density per transect & changing mean density into per 100 m^2 so that numbers are not all super small


df3$caribbeana <- df3$caribbeana/df3$Area * 100
df3$pennacea <- df3$pennacea/df3$Area * 100

# Mean density per site and year
caribbeana <- aggregate(caribbeana ~ Year + Site, data=df3, FUN=mean)
pennacea <- aggregate(pennacea ~ Year + Site, data=df3, FUN=mean)

# Adding in contextual variables

caribbeana <- left_join(caribbeana, covariates, by='Site')
pennacea <- left_join(pennacea, covariates, by='Site')



###
# Table 1 - Differences in density between years accounting for covariates
###

library('car')

caribbeana$adjusted <- (caribbeana$caribbeana)^(1/2)

# Caribbeana density
mod1 <- lm(adjusted ~ Year.x + Marine.Threat + Current + Coastal.Dev + Lat +
             Year.x:Marine.Threat + Year.x:Current + Year.x:Coastal.Dev + Year.x:Lat,
           data=caribbeana)

step(mod1)

mod2 <- lm(formula = adjusted ~ Year.x + Current + Lat + Year.x:Current, 
           data = caribbeana)

summary(mod2)
tidy(mod2)

write.csv(tidy(aov(mod2)), 'Table 1A - Caribbeana density.csv', row.names=FALSE)


# Pennacea density
pennacea$adjusted <- (pennacea$pennacea)^(1/2)

mod3 <- lm(adjusted ~ Year.x + Marine.Threat + Current + Coastal.Dev + Lat +
             Year.x:Marine.Threat + Year.x:Current + Year.x:Coastal.Dev + Year.x:Lat,
           data=pennacea)

step(mod3)

mod4 <- lm(formula = adjusted ~ Year.x + Current + Coastal.Dev + Lat + 
             Year.x:Current + Year.x:Coastal.Dev + Year.x:Lat, data = pennacea)

summary(mod4)
tidy(mod4)

write.csv(tidy(aov(mod4)), 'Table 1B - Pennacea density.csv', row.names=FALSE)


###
# Figure S1 - Black Coral - colony size with depth comparison
###

rm(list=ls())


df1 <- read.csv("Data S2.csv")
df1$Site <- as.character(df1$Site)

# Selecting just 1998
df2 <- df1[which(df1$Year==1998), ]

# Just caribbeana

df3 <- df2[which(df2$Species=='caribbeana'), ]

plot(df3$Height.cm ~ df3$Depth.m)
plot(df3$Width.cm ~ df3$Depth.m)

# Testing for relationship between colony height and depth
mod1 <- lm(df3$Height.cm ~ df3$Depth.m)
summary(mod1)

mod2 <- lm(df3$Width.cm ~ df3$Depth.m)
summary(mod2)

# Just pennacea

df4 <- df2[which(df2$Species=='pennacea'), ]

plot(df4$Height.cm ~ df4$Depth.m)
plot(df4$Width.cm ~ df4$Depth.m)


# Testing for relationship between colony height and depth
mod3 <- lm(Height.cm ~ Depth.m, data=df4)
summary(mod3)

height.pred <- data.frame(Depth.m=seq(min(df4$Depth.m), max(df4$Depth.m), length.out=500),
                          Height.cm=seq(min(df4$Height.cm, na.rm=TRUE), max(df4$Height.cm, na.rm=TRUE), length.out=500))

predictions <- predict.lm(mod3, height.pred, se.fit=TRUE, type='response')

height.pred$p.height <- predictions$fit
height.pred$p.upp <- predictions$fit + (predictions$se.fit * 2)
height.pred$p.lwr <- predictions$fit - (predictions$se.fit * 2)


# Testing for relationship between colony width and depth
mod4 <- lm(Width.cm ~ Depth.m, data=df4)


width.pred <- data.frame(Depth.m=seq(min(df4$Depth.m), max(df4$Depth.m), length.out=500),
                         Width.cm=seq(min(df4$Width.cm, na.rm=TRUE), max(df4$Width.cm, na.rm=TRUE), length.out=500))

predictions <- predict.lm(mod4, width.pred, se.fit=TRUE, type='response')

width.pred$p.width <- predictions$fit
width.pred$p.upp <- predictions$fit + (predictions$se.fit * 2)
width.pred$p.lwr <- predictions$fit - (predictions$se.fit * 2)



# Pennacea clonies get smaller on deeper reefs


png("Fig S1 - Colony size with depth.png", width = 6, height = 3, units = 'in', res = 300)


par(mfrow=c(1, 2))
par(mar=c(2, 2, 1, 0.5))
par(oma=c(2, 2, 0, 0))

# Pennacea height
plot(df4$Height.cm ~ df4$Depth.m, axes=FALSE, ylim=c(0, 200), cex=0.8, xlim=c(15,80))
axis(1, at=c(-10, 20, 30, 40, 50, 60, 70, 80, 100),
     labels=c("", 20, "", 40, "", 60, "", 80, ""))
axis(2)

lines(p.height ~ Depth.m, data=height.pred,
      col='red')

lines(p.upp ~ Depth.m, data=height.pred,
      col='red', lty=2)

lines(p.lwr ~ Depth.m, data=height.pred,
      col='red', lty=2)

legend("topleft", "A", bty='n', adj=c(2, -0.5))


# Pennacea width
plot(df4$Width.cm ~ df4$Depth.m, axes=FALSE, ylim=c(0, 200), cex=0.8, xlim=c(15, 80))
axis(1, at=c(-10, 20, 30, 40, 50, 60, 70, 80, 100),
     labels=c("", 20, "", 40, "", 60, "", 80, ""))
axis(2, lab=NA)

lines(p.width ~ Depth.m, data=width.pred,
      col='red')

lines(p.upp ~ Depth.m, data=width.pred,
      col='red', lty=2)

lines(p.lwr ~ Depth.m, data=width.pred,
      col='red', lty=2)

legend("topleft", "B", bty='n', adj=c(2, -0.5))

mtext('Depth (m)', side=1, line=1, outer=TRUE)
mtext('Colony Size (cm)', side=2, line=1, outer=TRUE)

dev.off()


###
# Fig S2 - Black Coral - height and width 50-60 m depth range only
###

rm(list=ls())


library(sm)
library(KernSmooth)

df1 <- read.csv("Data S2.csv")
df1$Site <- as.character(df1$Site)



###
# Length comparisons - just corals between 50-60 m depth
###

# selecting data within the 50-60 m depth range
df1 <- df1[which(df1$Depth > 49.9), ]
df1 <- df1[which(df1$Depth < 60.1), ]


###
# Caribbeana
###

df4 <- df1[which(df1$Species=="caribbeana"), ]

# selecting based on year
y1998 <- df4[which(df4$Year==1998), ]
y2016 <- df4[which(df4$Year==2016), ]

## Height

# selection of bandwidths
h1998.width <- dpik(y1998$Height.cm)
h2016.width <- dpik(y2016$Height.cm)

# generating the length distribution outlines
ch1998.outline <- bkde(y1998$Height.cm, bandwidth=h1998.width)
ch2016.outline <- bkde(y2016$Height.cm, bandwidth=h2016.width)

# concatinating all lengths
lengths <- c(y1998$Height.cm, y2016$Height.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

caribbeana.height <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                        nboot=9999, ngrid=500)

caribbeana.height.n <- length(lengths)

## Width

# selection of bandwidths
w1998.width <- dpik(y1998$Width.cm)
w2016.width <- dpik(y2016$Width.cm)

# generating the length distribution outlines
cw1998.outline <- bkde(y1998$Width.cm, bandwidth=w1998.width)
cw2016.outline <- bkde(y2016$Width.cm, bandwidth=w2016.width)

# concatinating all lengths
lengths <- c(y1998$Width.cm, y2016$Width.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

caribbeana.width <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                       nboot=9999, ngrid=500)

caribbeana.width.n <- length(lengths)

###
# Pennacea
###

df4 <- df1[which(df1$Species=="pennacea"), ]
df4 <- df4[(is.na(df4$Height.cm)==FALSE), ]

# selecting based on year
y1998 <- df4[which(df4$Year==1998), ]
y2016 <- df4[which(df4$Year==2016), ]

## Height

# selection of bandwidths
h1998.width <- dpik(y1998$Height.cm)
h2016.width <- dpik(y2016$Height.cm)

# generating the length distribution outlines
ph1998.outline <- bkde(y1998$Height.cm, bandwidth=h1998.width)
ph2016.outline <- bkde(y2016$Height.cm, bandwidth=h2016.width)

# concatinating all lengths
lengths <- c(y1998$Height.cm, y2016$Height.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

pennacea.height <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                      nboot=9999, ngrid=500)

pennacea.height.n <- length(lengths)

## Width

# selection of bandwidths
w1998.width <- dpik(y1998$Width.cm)
w2016.width <- dpik(y2016$Width.cm)

# generating the length distribution outlines
pw1998.outline <- bkde(y1998$Width.cm, bandwidth=w1998.width)
pw2016.outline <- bkde(y2016$Width.cm, bandwidth=w2016.width)

# concatinating all lengths
lengths <- c(y1998$Width.cm, y2016$Width.cm)
groups <- c(rep(1, nrow(y1998)), rep(2, nrow(y2016)))

pennacea.width <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                     nboot=9999, ngrid=500)

pennacea.width.n <- length(lengths)



###
# Plotting
###


png("Fig S2 - 50-60m size.png", width = 6, height = 6, units = 'in', res = 300)


par(mfrow=c(2, 2))
par(oma=c(3, 3, 0, 0))
par(mar=c(1, 1, 1, 1))

##
# Plotting Caribbeana height
##


plot(ch1998.outline, lwd=2, type='n', xlim=c(0, 250), ylim=c(0, 0.03),
     axes=FALSE)

axis(1, labels=FALSE)
axis(2, labels=TRUE)

polygon(c(caribbeana.height$eval.points, rev(caribbeana.height$eval.points)), c(caribbeana.height$upper, rev(caribbeana.height$lower)),
        col = "grey", border = NA)

lines(ch1998.outline, lwd=2, lty=1)
lines(ch2016.outline, lwd=2, lty=2)

text(210, 0.022, eval(paste("n = ", caribbeana.height.n)))
text(210, 0.020, eval(paste("p = ", round(caribbeana.height$p, 2))))
text(0, 0.029, substitute(paste("(A)  ", italic('A. caribbeana'), " height", sep=" ")), pos=4)


##
# Plotting pennacea Height
##


plot(ph1998.outline, lwd=2, type='n', xlim=c(0, 250), ylim=c(0, 0.03),
     axes=FALSE)


axis(1, labels=FALSE)
axis(2, labels=FALSE)

polygon(c(pennacea.height$eval.points, rev(pennacea.height$eval.points)), c(pennacea.height$upper, rev(pennacea.height$lower)),
        col = "grey", border = NA)

lines(ph1998.outline, lwd=2, lty=1)
lines(ph2016.outline, lwd=2, lty=2)

text(210, 0.022, eval(paste("n = ", pennacea.height.n)))
text(210, 0.020, eval(paste("p < 0.001 ")))
text(0, 0.029, substitute(paste("(B)  ", italic('P. pennacea'), " height", sep=" ")), pos=4)

legend(150, 0.018, legend=c("1998", "2016"), lty=c(1, 2), lwd=2)



####
# Plotting Caribbeana width
####


plot(cw1998.outline, lwd=2, type='n', xlim=c(0, 250), ylim=c(0, 0.02),
     axes=FALSE)


axis(1, labels=TRUE)
axis(2, labels=TRUE)

polygon(c(caribbeana.width$eval.points, rev(caribbeana.width$eval.points)), c(caribbeana.width$upper, rev(caribbeana.width$lower)),
        col = "grey", border = NA)

lines(cw1998.outline, lwd=2, lty=1)
lines(cw2016.outline, lwd=2, lty=2)

text(210, 0.017, eval(paste("n = ", caribbeana.width.n)))
text(210, 0.0155, eval(paste("p = ", round(caribbeana.width$p, 2))))
text(0, 0.0185, substitute(paste("(C)  ", italic('A. caribbeana'), " width", sep=" ")), pos=4)




##
# Plotting pennacea width
##


plot(pw1998.outline, lwd=2, type='n', xlim=c(0, 250), ylim=c(0, 0.02),
     axes=FALSE)


axis(1, labels=TRUE)
axis(2, labels=FALSE)

polygon(c(pennacea.width$eval.points, rev(pennacea.width$eval.points)), c(pennacea.width$upper, rev(pennacea.width$lower)),
        col = "grey", border = NA)

lines(pw1998.outline, lwd=2, lty=1)
lines(pw2016.outline, lwd=2, lty=2)

text(210, 0.017, eval(paste("n = ", pennacea.width.n)))
text(210, 0.0155, eval(paste("p < 0.001 ")))
text(0, 0.0185, substitute(paste("(D)  ", italic('P. pennacea'), " width", sep=" ")), pos=4)



# Axis labels

mtext("Black Coral size (cm)", 1, line=1.5, outer=TRUE)
mtext("Probability Density", 2, line=1.5, outer=TRUE)


dev.off()


###
# Figure S3 - Black Coral - colony frequency with depth
###

rm(list=ls())


df1 <- read.csv("Data S2.csv")
df1$Site <- as.character(df1$Site)

# Selecting just 1998
df2 <- df1[which(df1$Year==1998), ]




png("Fig S3 - Colony frequency with depth.png", width = 6, height = 3, units = 'in', res = 300)

par(mfrow=c(1, 2))
par(mar=c(2, 2, 1, 0.5))
par(oma=c(2, 2, 0, 0))

# Just caribbeana

df3 <- df2[which(df2$Species=='caribbeana'), ]

hist(df3$Depth, breaks=c(0, 10, 20, 30, 40, 50, 60, 70, 80), main="", ylim=c(0, 20))

legend("topleft", "A", bty='n', adj=c(2, -0.5))



# Just pennacea

df4 <- df2[which(df2$Species=='pennacea'), ]

hist(df4$Depth, breaks=c(0, 10, 20, 30, 40, 50, 60, 70, 80), main='', ylim=c(0, 100))


legend("topleft", "B", bty='n', adj=c(2, -0.5))

mtext('Depth (m)', side=1, line=1, outer=TRUE)
mtext('Frequency', side=2, line=1, outer=TRUE)

dev.off()