Sediment cover analysis
Load data
# Load CPCe data (see tidy_cpce_data.Rmd)
load("data/processed/cpce.RData")
# Rename dataset to 'cpce'
cpce <- cpce_subQuantify sand cover from cpce data 2013-2016
# Count sand vs. not sand, and create numeric weeks column
cpce <- cpce %>%
mutate(sumpts = rowSums(select_if(., is.integer)), # sum counts of all categories
notsand = as.integer(sumpts - sand_sa)) %>% # sum counts of not sand
mutate_if(is.character, as.factor) %>%
mutate(weekn = (date - min(date)) / dweeks(1), # create numeric weeks column
dayn = as.numeric(date - min(date)))Visualize mean sand cover over time at each site
# Calculate mean proportion sand cover at each site on each date
# (averaged across all frames from all transects)
sand.summ <- cpce %>%
mutate(propsand = sand_sa / sumpts) %>%
group_by(reef, dir, channel, site, date) %>%
summarise(meanpropsand = mean(propsand)) %>%
ungroup()
# Plot by site
ggplot(sand.summ, aes(x = date, y = meanpropsand, color = channel, group = site)) +
facet_wrap(~ site) +
geom_point() +
geom_line() +
geom_vline(xintercept = c(dredge.start.date, dredge.end.date), linetype = 3, lwd = 0.25)
# Plot by reef and direction
ggplot(sand.summ, aes(x = date, y = meanpropsand, color = channel, group = site,
linetype = channel)) +
facet_grid(dir ~ reef) +
geom_point() +
geom_line() +
geom_vline(xintercept = c(dredge.start.date, dredge.end.date), linetype = 1, lwd = 0.25)
Filter out HBN1 and HBN2
Sites HBN1 and HBN2 have a high amount of sand cover to begin with, and they are classified in DCA reports as “Scattered Coral/Rock in Sand”. The rest of the sites are all categorized as “Colonized Pavement”, “Ridge Reef”, “Linear Reef”, or “Spur and Groove”, and generally have lower amounts of sand. Therefore, the HBN1 and HBN2 sites should be omitted.
cpce <- cpce %>%
filter(!site %in% c("HBN1", "HBN2"))
# Calculate mean proportion sand cover at each site on each date
# (averaged across all frames from all transects)
sand.summ <- cpce %>%
mutate(propsand = sand_sa / sumpts) %>%
group_by(reef, dir, channel, site, date) %>%
summarise(meanpropsand = mean(propsand)) %>%
ungroup()Model sand cover
Fit GAMMs for each reef and direction (site as random effect)
# Filter dataset to only fit over regions with less than 30-week gaps between data points
cpce.f <- cpce %>%
nest(-reef, -dir, -channel) %>% # For each reef
mutate(dates = map(data, ~ sort(unique(.$date))), # Which weeks have data
runs = map(dates, ~ split(.x, cumsum(c(TRUE, diff(.x) > 140)))), # Find runs w/o 20wk gaps
longest = map(runs, ~ .x[[which.max(sapply(.x, length))]]), # Get longest run
data.f = map2(data, longest, ~ filter(.x, .$date %in% .y))) %>% # Get data from those dates
unnest(data.f)
# Add up the total number of points counted as sand and not sand across all frames from each transect/date
cpce.f <- cpce.f %>%
group_by(reef, dir, channel, site, transect, dayn) %>%
summarise(sand = sum(sand_sa),
notsand = sum(notsand)) %>%
ungroup()
# Fit GAMM models for each reef and direction
set.seed(456)
gamms.reef.dir <- cpce.f %>%
#filter(reef == "R2", dir == "S") %>% # UNCOMMENT TO LOOK ONLY AT R2S (TESTING...)
nest(-reef, -dir) %>% # Nest by reef and direction
mutate(
# Fit gamm by channel (in parallel)
mod = future_map(data,
~ gamm4(cbind(sand, notsand) ~ s(dayn, by = channel) + channel, ##, k = 8, bs = "cs"
random = ~(1|site/transect), data = ., family = "binomial")$gam),
# Get newdata frame for predict() using the appropriate range of weeks
newdata = map(data, ~ expand.grid(dayn = min(.$dayn):max(.$dayn),
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)))
)
gamms.reef.dir.fits <- gamms.reef.dir %>%
unnest(newdata, fit, lci, uci) %>%
mutate(date = min(cpce$date) + dayn)
save(sand.summ, gamms.reef.dir.fits, file = "data/processed/sedcov.gamms.RData")Plot fitted models
ggplot(gamms.reef.dir.fits, aes(x = date, y = fit, color = channel)) +
facet_grid(dir ~ reef) +
geom_point(data = sand.summ, aes(y = meanpropsand), 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 = 3, lwd = 0.25) +
labs(title = "Sand cover", x = "Date", y = "Percent sand")
Assess ‘permanent’ effects
# Show percent sand at each site for TRUE baseline (must be before Nov. 20, 2013) and late 2016
ba <- cpce %>%
filter(date < as.Date("2013-11-20") | date > as.Date("2016-01-01")) %>%
group_by(reef, dir, channel, site, date, transect) %>%
summarize(propsand = mean(sand_sa / sumpts)) %>%
group_by(reef, dir, channel, site, year = factor(year(date))) %>%
summarize(meanpropsand = mean(propsand),
sd = sd(propsand))
ggplot(ba, aes(x = year, y = meanpropsand, fill = year)) +
facet_wrap(~ site) +
geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = meanpropsand - sd, ymax = meanpropsand + sd),
width = 0.5, lwd = 0.3) +
scale_fill_manual(values = c("cornflowerblue", "coral")) +
theme(legend.position = c(0.75, 0.075))
There are only three sites for which CPCE data exist both before dredging began, and then in 2016: R2N1, R2S1, and R2SC1.
ba2 <- ba %>%
filter(site %in% c("R2N1", "R2S1", "R2SC1"))
ba2 %>%
ggplot(aes(x = year, y = meanpropsand, fill = year)) +
facet_wrap(~ site) +
geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = meanpropsand - sd, ymax = meanpropsand + sd),
width = 0.5, lwd = 0.3) +
scale_fill_manual(values = c("cornflowerblue", "coral")) +
theme(legend.position = c(0.75, 0.075))
ba2## # A tibble: 6 x 7
## # Groups: reef, dir, channel, site [3]
## reef dir channel site year meanpropsand sd
## <fct> <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 R2 N channelside R2N1 2013 0.0126 0.00854
## 2 R2 N channelside R2N1 2016 0.194 0.00175
## 3 R2 S channelside R2S1 2013 0.00216 0.00206
## 4 R2 S channelside R2S1 2016 0.340 0.0284
## 5 R2 S control R2SC1 2013 0.218 0.0119
## 6 R2 S control R2SC1 2016 0.215 0.0822# Percent increase at R2N1
ba2 %>%
filter(site == "R2N1") %>%
summarise(fc = meanpropsand[year == 2016] / meanpropsand[year == 2013] - 1) %>%
select(fc)## # A tibble: 1 x 4
## # Groups: reef, dir, channel [1]
## reef dir channel fc
## <fct> <fct> <fct> <dbl>
## 1 R2 N channelside 14.4# Percent increase at R2S1
ba2 %>%
filter(site == "R2S1") %>%
summarise(fc = meanpropsand[year == 2016] / meanpropsand[year == 2013] - 1) %>%
select(fc)## # A tibble: 1 x 4
## # Groups: reef, dir, channel [1]
## reef dir channel fc
## <fct> <fct> <fct> <dbl>
## 1 R2 S channelside 156.Mean daily percent sediment cover by distance from channel
dates1 <- gamms.reef.dir.fits %>%
filter(date > dredge.start.date, date < dredge.end.date) %>%
group_by(reef, dir, channel) %>%
summarise(d = list(unique(as.Date(date))))
dates_all <- zoo::as.Date.numeric(Reduce(intersect, dates1$d))
sandcov <- gamms.reef.dir.fits %>%
filter(date %in% dates_all) %>%
add_dist_cat() %>%
group_by(dist_cat) %>%
summarise(mean_sandcov = mean(fit),
sd_sandcov = sd(fit))
sandcov## # A tibble: 3 x 3
## dist_cat mean_sandcov sd_sandcov
## <fct> <dbl> <dbl>
## 1 < 50 m 0.611 0.160
## 2 1.25 - 2.5 km 0.358 0.216
## 3 9.38 km 0.154 0.0846