Setup

# Load packages
library(tidyverse)
library(lme4)
library(lsmeans)
library(janitor)

Load tidy sediment depth data

load(file = "data/processed/sediment_depth.RData")

Visualize all sediment depth measurements by reef area

ggplot(data = substrate, aes(x = dist_actual, y = sdepth_cm)) +
  facet_grid(dir ~ reef + reef2, scales = "free") + 
  geom_point()

Visualize data within 1000 meters of channel

All points

# Build ggplot
sdep <- ggplot(
  data = filter(substrate, dist < 1000), 
  mapping = aes(x = dist, y = sdepth_cm)) +
  facet_grid(dir ~ reef + reef2, scales = "free")

# Plot as points
sdep + geom_point()

Boxplots

sdep + geom_boxplot(aes(group = dist), width = 50)

RR showed lower sediment depth than LR, probably because it is a higher relief habitat type. R2 was sampled at LR and RR, but R3 was only sampled at LR. Therefore, comparisons between R2 and R3 are probably more accurate if only LR is included. R3 was also only sampled up to 500 m from the channel, while R2 was sampled out to ~800m. Comparisons should only include out to 500m.

# Subset LR habitat
LR <- substrate %>%
  filter(reef == "HB" | reef2 == "LR", dist < 1000)

# Build ggplot
sdep2 <- ggplot(data = LR, mapping = aes(x = dist, y = sdepth_cm)) +
  facet_grid(dir ~ reef) + 
  #geom_boxplot(aes(group = dist)) +
  geom_jitter(width = 20, size = 0.5, alpha = 0.5) +
  #geom_point(size = 0.5, alpha = 0.5) +
  coord_cartesian(xlim = c(0, 500))

sdep2

Visualize quantiles

# Get quantiles of sediment depth for each transect
df <- LR %>%
  filter(reef %in% c("R2", "R3")) %>% droplevels() %>%
  nest(-reef, -dir, -dist, -transect) %>%
  mutate(value = map(data, ~ quantile(.$sdepth_cm, c(0.99, 0.9, 0.9, 0.75, 0.5))),
         quantile = map(value, ~ factor(names(.)))) %>%
  unnest(quantile, value)

# Fit sediment depth quantiles as a function of distance from channel
mod <- lm(value ~ reef * dir * quantile * log(dist), data = df)
resdf <- tidyr::expand(df, crossing(reef, dir, quantile), dist = 15:500)
resdf$fit <- predict(mod, resdf)
## Convert depth quantile to "percent of reef covered"
resdf$pct_reef_covered <- factor(100 - as.numeric(str_sub(resdf$quantile, 1, 2)))

save(LR, resdf, file = "data/processed/sdepthmod.RData")

# Plot all sed depth measurements with fitted quantiles
ggplot(data = LR, mapping = aes(x = dist, y = sdepth_cm)) +
  facet_grid(dir ~ reef) + 
  geom_jitter(width = 10, size = 0.5, alpha = 0.5) +
  coord_cartesian(xlim = c(0, 500)) +
  geom_line(data = filter(resdf, !(reef == "R3" & dir == "S")),
                          aes(x = dist, y = fit, linetype = pct_reef_covered)) +
  labs(x = "Distance from channel (m)", y = "Sediment depth (cm)") +
  scale_linetype_manual(name = "% of area",
                        labels = c(1, 10, 25, 50), values = c(3,4,2,1)) +
  theme(legend.position = c (0.75, 0.3),
        legend.key.size = unit(0.25, "cm"),
        legend.key.width = unit(1,"cm"))

Quantile fitted at the 1 cm threshold

# Get quantiles of sediment depth for each transect
df <- LR %>%
  filter(reef %in% c("R2", "R3")) %>% droplevels() %>%
  nest(-reef, -dir, -dist, -transect) %>%
  mutate(value = map(data, ~ quantile(.$sdepth_cm, seq(0, 1, 0.005))),
         quantile = map(value, ~ factor(names(.)))) %>%
  unnest(quantile, value)

# Fit sediment depth quantiles as a function of distance from channel
mod <- lm(value ~ reef * dir * quantile * log(dist), data = df)
resdf <- tidyr::expand(df, crossing(reef, dir, quantile), dist = 15:500)
resdf$fit <- predict(mod, resdf)
## Convert depth quantile to "percent of reef covered"
resdf$pct_reef_covered <- factor(100 - as.numeric(str_sub(resdf$quantile, 1, 2)))

# Percent of reef covered in > 1 cm at 20m from channel
resdf %>%
  filter(dist == 20) %>%
  group_by(reef, dir) %>%
  mutate(thresh = abs(fit - 1)) %>%
  filter(thresh == min(thresh))
## # A tibble: 4 x 7
## # Groups:   reef, dir [4]
##   reef  dir   quantile  dist   fit pct_reef_covered   thresh
##   <fct> <fct> <fct>    <int> <dbl> <fct>               <dbl>
## 1 R2    N     53.5%       20 1.01  47               0.00767 
## 2 R2    S     94.0%       20 0.992 6                0.00827 
## 3 R3    N     44.5%       20 0.998 56               0.00167 
## 4 R3    S     86.5%       20 0.999 14               0.000654
# Percent of reef covered in > 1 cm at 300m from channel
resdf %>%
  filter(dist == 500) %>%
  group_by(reef, dir) %>%
  mutate(thresh = abs(fit - 1)) %>%
  filter(thresh == min(thresh))
## # A tibble: 4 x 7
## # Groups:   reef, dir [4]
##   reef  dir   quantile  dist   fit pct_reef_covered  thresh
##   <fct> <fct> <fct>    <int> <dbl> <fct>              <dbl>
## 1 R2    N     88.0%      500 1.01  12               0.00687
## 2 R2    S     87.0%      500 1.00  13               0.00430
## 3 R3    N     84.0%      500 0.969 16               0.0307 
## 4 R3    S     81.0%      500 1.00  19               0.00374