Setup

# Load packages
library(tidyverse)
library(janitor)
library(lubridate)
library(lme4)
library(lsmeans)
library(ggrepel)
library(ggthemes)

Import data

# Import tidy tagged coral data
load("data/processed/tagged.RData")
tagged_bool <- tagged_bool %>%
  mutate(day = (date - min(date)) / ddays(1))

Prevalence of partial mortality due to sedimentation

The DCA November 2015 report (page 14) indicates that the “PM” condition code is used to indicate partial mortality associated with sedimentation. Analysis of the prevalence of “PM” codes recorded for all tagged corals will test the hypothesis that corals at channelside sites were more likely to experience partial mortality due to sedimentation than corals at control sites.

First we analyze the accumulation of the PM condition code across tagged corals at all sites. For this analysis, once a coral is recorded as PM, it is counted as PM at all future time points.

# Fill in missing observations (if a specific coral was missed on a day) 
tagged_bool_complete <- tagged_bool %>% 
  group_by(site) %>%
  complete(date, coral, fill = list(PM = FALSE)) %>%
  fill(reef, dir, channel) %>%
  ungroup()

# If PM ever observed for a coral, change all future observations of PM to TRUE
cumpm <- tagged_bool_complete %>%
  group_by(coral) %>%
  mutate(PM = ifelse(any(PM) & date > min(date[PM]), TRUE, PM)) %>%
  group_by(date, reef, dir, channel, site) %>%
  summarise(prop = sum(PM) / n())

# Plot cumulative proportion of corals with PM (ever) over time
#ggplot(cumpm, aes(x = date, y = prop)) + 
#  facet_wrap(~ site) +
#  geom_line(aes(color = channel, group = channel)) +
#  labs(title = "Cumulative PM", x = "Time (weeks)", y = "Proportion of corals")

ggplot(cumpm, aes(x = date, y = prop, color = channel, group = site)) + 
  geom_line() +
  labs(title = "Cumulative PM", x = "Date", y = "Proportion of corals") +
  geom_text_repel(data = cumpm %>% group_by(site) %>% filter(date == max(date)),
                  aes(label = site))

Filter dataset

Omit site HBN1, because this site was not monitored during much of the construction period because it had been ‘buried’ by a natural sand wave (DCA). Many of the corals were recorded as BUR before observations of this site ended.

tagged_bool <- tagged_bool %>%
  filter(site != "HBN1")                     # Remove site because corals all BUR => no opp for PM?

Calculate the total number of corals at each site that experienced partial mortality (up to end of March 2015, last month of dredging), and plot as a barplot.

# Quantify occurrence of PM during compliance and post-con:
pm <- tagged_bool %>%
  #filter(date < "2015-07-31") %>%     
  filter(date <= "2015-03-31") %>%     # Filter up through March 2015 (last month of dredging)
  group_by(reef, dir, channel, site, transect, species, coral) %>%
  summarise(PM = any(PM))

# Summarize proportion of tagged corals at each site that experienced PM
pm.summ <- pm %>%
  group_by(site, channel, reef, dir) %>%
  summarise(n = n(),
            pm = sum(PM),
            prop = pm / n)

ggplot(pm.summ, aes(x = site, y = prop, fill = interaction(reef, dir, channel))) +
  facet_wrap(dir ~ reef, scales = "free_x") +
  geom_bar(stat = "identity")

This plot of the proportion of tagged corals at each site show that channelside sites generally showed higher levels of partial mortality. We can test for a statistical effect of location using a generalized linear mixed model.

Was partial mortality due to sedimentation more common at channelside sites?

# Fit generalized linear mixed model
pm.mod1 <- glmer(PM ~ channel * reef * dir + (1|site/transect) + (1|species), 
                 family = "binomial", data = pm)

# Generate lsmeans for each reef/direction/channelside
pm.lsm1 <- lsmeans(pm.mod1, specs = c("channel", "reef", "dir"), type = "response")

# Test channelside vs. control mortality for each reef/direction
pm.lsm1.sig <- rbind(pairs(pm.lsm1, by = c("reef", "dir")), adjust = "none")
knitr::kable(summary(pm.lsm1.sig))
reef dir contrast odds.ratio SE df z.ratio p.value
HB N channelside - control 2.254076 1.805157 NA 1.014858 0.3101733
R2 N channelside - control 74.202913 51.854952 NA 6.162909 0.0000000
R3 N channelside - control 4.623587 3.735692 NA 1.895098 0.0580794
HB S channelside - control 7.918791 5.087288 NA 3.220944 0.0012777
R2 S channelside - control 8.333743 4.902544 NA 3.604280 0.0003130
R3 S channelside - control 3.731173 2.139299 NA 2.296509 0.0216468 
# Plot partial mortality prevalence
ggplot(summary(pm.lsm1), aes(x = channel, y = prob, fill = channel)) + 
  facet_grid(dir ~ reef) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2, lwd = 0.3) +
  geom_text(data = data.frame(summary(pm.lsm1.sig)), inherit.aes = FALSE,
            aes(x = 1.5, y = 0.8, label = ifelse(p.value < 0.05, "*", NA)), cex = 10)

save(pm.lsm1, pm.lsm1.sig, file = "data/processed/pm.lsm.RData")

We find a significantly higher probability of partial mortality due to sedimentation for corals at channelside sites relative to controls for the northern middle reef (75% vs. 4%, p < 0.0001), the southern middle reef (61% vs. 16%, p = 0.0003), the southern hardbottom (70% vs. 23%, p = 0.0013), and the southern outer reef (24% vs. 8%, p = 0.0216. There was an almost significant difference in probability of partial mortality at the northern outer reef channelside sites (68%) vs. control sites (31%; p = 0.0581).