Analysis of dredge plume at permanent monitoring sites

Import dredge plume data

plume <- readxl::read_xls(path = "data/plume/dredge_plume.xls", sheet = "plume presence absence") %>%
  gather(date, plume, -NAME, -LATITUDE, -LONGITUDE, -idx) %>%
  clean_names() %>%
  mutate(site = as_factor(name),
         date = as_date(date),
         plume = as.logical(plume),
         reef = str_sub(site, 1, 2),
         dir  = str_sub(site, 3, 3),
         channel = if_else(grepl("C", site), "control", "channelside")) %>%
  filter(date > dredge.start.date, date < dredge.end.date)



plume <- plume %>%
  filter(!grepl("50", name), !grepl("75", name), !(grepl("00", name)))# %>%


ggplot(plume, aes(x = date, y = plume)) +
  facet_wrap(~ name) +
  geom_point()

# Filter data to only includ unique pixels (no duplicates, e.g. if multiple sites were in the same pixel)
plume.f <- plume %>% distinct(idx, date, plume, reef, dir, channel)

# Fit GAMMs for each reef --- but only use multiple sites if diff pixels
set.seed(9483)
gamms1 <- plume.f %>%
  filter(reef %in% c("HB", "R2", "R3"), dir %in% c("N", "S")) %>%
  # For each reef and direction
  nest(-reef, -dir, -channel) %>%
  mutate(
    # Fit a binomial GAMM for sediment stress over time at channelside vs. control, site as random
    mod = future_map(data, ~ gam(plume ~ s(as.numeric(date)), family = "binomial", data = .)),
    # Get newdata frame for predict() using the appropriate range of weeks
    newdata = map(data, ~ expand.grid(date = as.numeric(min(.$date):max(.$date)),
                                      channel = c("channelside", "control"))),
    # 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 = map(sim, ~ logit2prob(apply(., 1, quantile, 0.025))),
    uci = map(sim, ~ logit2prob(apply(., 1, quantile, 0.975)))
  )


gamm.fits <- gamms1 %>% 
  unnest(newdata, fit, lci, uci) %>%
  mutate(date = zoo::as.Date.numeric(date))

### Plot fitted models
ggplot(gamm.fits, aes(x = date, y = fit, color = channel)) +
  facet_grid(dir + channel ~ reef) +
  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(x = "Date", y = "Probability dredge plume present") +
  geom_point(data = filter(plume, reef %in% c("HB", "R2", "R3"), dir %in% c("N", "S")),
             aes(y = as.numeric(plume)), pch = "|")

# Summarize plume presence
plume.summ <- plume %>%
  group_by(reef, dir, channel) %>%
  summarise(nobs = n(), nyes = sum(plume), nno = nobs - nyes, pct = nyes / nobs)

save(plume, plume.summ, gamm.fits, file = "data/processed/plume.RData")

Reported metrics

# Plume frequency in late 2013
gamm.fits %>%
  filter(date < "2014-01-01") %>%
  group_by(reef, dir, channel) %>%
  summarise(mean = mean(fit)) %>%
  arrange(-mean)

## # A tibble: 12 x 4
## # Groups:   reef, dir [6]
##    reef  dir   channel      mean
##    <chr> <chr> <chr>       <dbl>
##  1 HB    N     channelside 0.918
##  2 HB    S     channelside 0.876
##  3 R2    N     channelside 0.783
##  4 R2    S     channelside 0.768
##  5 R3    N     channelside 0.584
##  6 R3    S     channelside 0.582
##  7 R2    S     control     0.503
##  8 HB    N     control     0.497
##  9 R3    S     control     0.439
## 10 HB    S     control     0.326
## 11 R2    N     control     0.171
## 12 R3    N     control     0.139

# Plume frequency in 2014 - 2015
gamm.fits %>%
  filter(date >= "2014-01-01") %>%
  group_by(reef, dir, channel) %>%
  summarise(mean = mean(fit)) %>%
  arrange(-mean)

## # A tibble: 12 x 4
## # Groups:   reef, dir [6]
##    reef  dir   channel       mean
##    <chr> <chr> <chr>        <dbl>
##  1 HB    N     channelside 0.651 
##  2 R2    S     channelside 0.640 
##  3 R2    N     channelside 0.630 
##  4 HB    S     channelside 0.562 
##  5 R3    S     channelside 0.525 
##  6 R3    N     channelside 0.525 
##  7 R2    S     control     0.376 
##  8 HB    N     control     0.282 
##  9 R3    S     control     0.277 
## 10 HB    S     control     0.164 
## 11 R3    N     control     0.0418
## 12 R2    N     control     0.0312