Load sediment trap data

load("data/processed/sedtrap.RData")

# Convert g/d to mg/cm2/d (1" inner diameter PVC = 5.067 square centimeters)
sedtrap <- sedtrap %>%
  mutate(
    sedrate_mgcm2d = sedrate * 1000 / 5.067,
    finerate_mgcm2d = finerate * 1000 / 5.067
    )

Plot sedimentation rate over time for each permanent monitoring site

# Get average sedimentation rates for each permanent monitoring site
sedrate <- sedtrap %>%
  group_by(reef, dir, channel, site, middate) %>%
  summarise(meansed = mean(sedrate_mgcm2d),
            meanfine = mean(finerate_mgcm2d),
            sdsed = sd(sedrate_mgcm2d),
            sdfine = sd(finerate_mgcm2d),
            sesed = sdsed / sqrt(n()),
            sefine = sdfine / sqrt(n()))

# Plot average fine sedimentation rates for each site
ggplot(sedrate, aes(x = middate, y = meanfine)) +
  facet_wrap(~ site) +
  #geom_point(size = 1) +
  geom_errorbar(aes(ymin = meanfine - sdfine, ymax = meanfine + sdfine)) +
  geom_line() + 
  geom_vline(xintercept = c(dredge.start.date, dredge.end.date), linetype = 1, lwd = 0.25) +
  labs(title = "Fine sediments", x = "Date", y = "mg / cm2 / day")

Filter dataset

sedtrap <- sedtrap %>%
  filter(site != "HBN1")             # Omit HBN1 because was "buried" and not monitored for a long time
                                     # Both here and in tagged coral data

Plot sedimentation rate over time for each reef and direction

# Get average sedimentation rates for each reef/direction
sedrate <- sedtrap %>%
  group_by(reef, dir, channel, start, end, middate) %>%
  summarise(meansed = mean(sedrate_mgcm2d),
            meanfine = mean(finerate_mgcm2d),
            sdsed = sd(sedrate_mgcm2d),
            sdfine = sd(finerate_mgcm2d),
            sesed = sdsed / sqrt(n()),
            sefine = sdfine / sqrt(n())) %>%
  ungroup()

# Plot fine sediments
ggplot(sedrate, aes(x = middate, y = meanfine, color = channel)) +
  facet_grid(dir ~ reef) +
  geom_point(size = 1) +
  geom_errorbar(aes(ymin = meanfine - sdfine, ymax = meanfine + sdfine)) +
  geom_line() +
  labs(title = "Fine sediments", x = "Date", y = "mg / cm2 / day")

GAM models

Fine sediments

# No random effects -- each sediment trap within a monitoring area is treated as an independent replicate
set.seed(4324)
fine.gamms.reef.dir <- sedtrap %>%
  nest(-reef, -dir, -channel) %>%   # Nest by reef and direction
  mutate(
    # Fit gamm by channel (in parallel)
    mod = future_map(data,       
            ~ gam(total_dry * frac_fine/100 ~ s(as.numeric(middate)) +
                   offset(log(duration)),
                   data = ., family = "poisson",
                   verbosePQL = FALSE)),
    # Get newdata frame for predict() using the appropriate range of days
    newdata = map(data, ~ expand.grid(middate = as.numeric(min(.$middate):max(.$middate)),
                                      #channel = c("channelside", "control"),
                                      duration = 1)),
    # Get fitted values from GAMM using newdata frame
    fit = map2(mod, newdata, ~ predict(object = .x, newdata = .y, type = "response")),
    # Get confidence intervals by simulation
    Rbeta = map(mod, ~ mvrnorm(n = 1000, coef(.), vcov(.))),
    Xp = map2(mod, newdata, ~ predict(object = .x, newdata = .y, type = "lpmatrix")),
    sim = map2(Xp, Rbeta, function(x, y) x %*% t(y)),
    lci = map2(mod, sim, ~ .x$family$linkinv(apply(.y, 1, quantile, 0.025))),
    uci = map2(mod, sim, ~ .x$family$linkinv(apply(.y, 1, quantile, 0.975)))
  )

fine.gamms.reef.dir.fits <- fine.gamms.reef.dir %>% 
  unnest(newdata, fit, lci, uci) %>%
  mutate(middate = zoo::as.Date.numeric(middate)) %>%
  # Convert to mg per cm2 per day (from g per day)
  mutate_at(vars(fit, lci, uci), function(x) x * 1000 / 5.067)

save(fine.gamms.reef.dir.fits, sedtrap, file = "data/processed/sedtrap.fine.gamms.RData")

### Plot fitted models
ggplot(fine.gamms.reef.dir.fits, aes(x = middate, y = fit, color = channel)) +
  facet_grid(dir ~ reef) +
  geom_segment(data=sedtrap, aes(x = start, xend = end, 
                                 y = finerate_mgcm2d, yend = finerate_mgcm2d), alpha=0.4) +
  geom_ribbon(aes(ymin = lci, ymax = uci, linetype = NA), alpha = 0.2) + 
  geom_line(lwd = 1) +
  geom_vline(xintercept = c(dredge.start.date, dredge.end.date), linetype = 1, lwd = 0.25) +
  labs(title = "Fine sediments", x = "Date", y = "Sedimentation rate (mg/cm2/d)")

# Integrate to get total sediment deposited
## For reef 3 control sites, the earliest sediment trap sample had a midpoint of january 10-11, 2014. This cuts out almost 2 months of data after dredging began 11/20/13 from other sites for which we do have information on sedimentation rates. Therefore, to avoid excluding this data, we assume that the sedimentation rate at these sites from 11/20/13 until data are available is the same as the sedimentation rate at the beginning of data.
fine.gamms.reef.dir.fits.complete <- fine.gamms.reef.dir.fits %>%
  complete(nesting(reef, dir, channel), middate) %>%
  group_by(reef, dir, channel) %>%
  fill(fit, lci, uci, .direction = "up") %>%
  ungroup()

fine.depo <- fine.gamms.reef.dir.fits.complete %>%
  filter(middate >= "2013-11-20",              # Filter from start of dredging
         middate <= "2015-03-23") %>%          # Filter up to a week after date dredging ended
  group_by(reef, dir, channel) %>%
  summarise(int_min = sum(lci),
            int_fit = sum(fit),
            int_max = sum(uci))

ggplot(fine.depo, aes(x = channel, y = int_fit/1000, fill = channel)) +
  facet_grid(dir ~ reef) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = int_min/1000, ymax = int_max/1000), width = 0.25, lwd = 0.25) +
  labs(title = "Fine sediment deposition 2013-11-20 to 2015-03-23",
       y = "Grams / cm2", x = "")

Save sediment integrations as .RData

save(fine.depo, file = "data/processed/sed_deposited.RData")

Reported results

# Highest measurement at northern middle reef channelside
sedtrap %>% filter(finerate_mgcm2d > 600) %>% select(site, start, end, middate, duration, finerate_mgcm2d)
## # A tibble: 1 x 6
##   site  start      end        middate    duration finerate_mgcm2d
##   <chr> <date>     <date>     <date>        <dbl>           <dbl>
## 1 R2N2  2014-10-18 2014-12-08 2014-11-13       51            606.
# Measurement at northern middle reef control at the same time 
fine.gamms.reef.dir.fits %>% filter(reef == "R2", dir == "N", channel == "control",
                                    middate == "2014-11-13")
## # A tibble: 1 x 8
##   reef  dir   channel   fit   lci   uci middate    duration
##   <fct> <fct> <fct>   <dbl> <dbl> <dbl> <date>        <dbl>
## 1 R2    N     control  58.3  48.2  70.5 2014-11-13        1
# Final measurements
fine.gamms.reef.dir.fits %>%
  group_by(reef, dir, channel) %>%
  summarise(lastdate = last(middate),
            lastrate = fit[middate == lastdate])
## # A tibble: 12 x 5
## # Groups:   reef, dir [?]
##    reef  dir   channel     lastdate   lastrate
##    <fct> <fct> <fct>       <date>        <dbl>
##  1 HB    N     channelside 2015-07-06    39.3 
##  2 HB    N     control     2015-07-06    35.1 
##  3 HB    S     channelside 2015-07-03    30.8 
##  4 HB    S     control     2015-07-06    11.1 
##  5 R2    N     channelside 2015-07-01    32.0 
##  6 R2    N     control     2015-07-01    13.8 
##  7 R2    S     channelside 2015-07-01    35.6 
##  8 R2    S     control     2015-07-01    17.6 
##  9 R3    N     channelside 2015-07-04    13.6 
## 10 R3    N     control     2015-07-04     8.69
## 11 R3    S     channelside 2015-07-04    13.3 
## 12 R3    S     control     2015-07-04     7.47
# Mean throughout dredging
fine.gamms.reef.dir.fits %>%
  filter(middate > "2013-11-20", middate < "2015-03-16") %>%
  group_by(reef, dir, channel) %>%
  summarise(meanrate = mean(fit),
            medrate  = median(fit))
## # A tibble: 12 x 5
## # Groups:   reef, dir [?]
##    reef  dir   channel     meanrate medrate
##    <fct> <fct> <fct>          <dbl>   <dbl>
##  1 HB    N     channelside    117.   118.  
##  2 HB    N     control        105.    94.7 
##  3 HB    S     channelside    118.   116.  
##  4 HB    S     control         69.4   48.5 
##  5 R2    N     channelside    170.   155.  
##  6 R2    N     control         25.4   17.4 
##  7 R2    S     channelside    109.   105.  
##  8 R2    S     control         58.1   52.7 
##  9 R3    N     channelside     52.4   59.7 
## 10 R3    N     control         11.1    8.96
## 11 R3    S     channelside     36.9   23.5 
## 12 R3    S     control         24.2   11.3
load("data/processed/sed_deposited.RData")
load("data/processed/sedtrap.fine.gamms.RData")

add_dist_cat <- function(x) {
  x %>% 
    mutate(dist_cat = factor(case_when(
      channel == "channelside" ~ "< 50 m",
      reef == "HB" & channel == "control" ~ "1.25 - 2.5 km",
      reef %in% c("R2", "R3") & channel == "control" & dir == "S" ~ "1.25 - 2.5 km",
      reef %in% c("R2", "R3") & channel == "control" & dir == "N" ~ "9.38 km")))
}


fine.gamms.reef.dir.fits %>%
  filter(middate >= "2013-11-20", middate <= "2015-03-16") %>%
  nest(-reef, -dir, -channel) %>%
  mutate(       days = map(data, nrow),
              over25 = map(data, ~rle(.$fit > 25)),         # flag all days > 25 mg/cm2 threshold
          over25runs = map(over25, ~ .$lengths[.$values]),  # count consec days above threshold
         over25for30 = map(over25runs, ~ .[. > 30]),        # get periods over threshold for >30 days
         overthresh  = map(over25for30, ~ sum(.))) %>%      # Sum length of these periods 
  unnest(days, overthresh) %>%
  mutate(pctoverthresh = overthresh / days) %>%             # Percent of days in these periods over threshold
  arrange(pctoverthresh) %>%
  add_dist_cat() %>%
  group_by(dist_cat) %>%
  summarise(mean_pctoverthresh = mean(pctoverthresh))
## # A tibble: 3 x 2
##   dist_cat      mean_pctoverthresh
##   <fct>                      <dbl>
## 1 < 50 m                     0.840
## 2 1.25 - 2.5 km              0.667
## 3 9.38 km                    0.152