Coral count data analysis
Import coral count data
# Import tidy count data
load("data/processed/all_counts.RData")
# Subset Scleractinian count data
scler <- filter(all_counts, category == "scleractinian") %>% # Get scleractinians only
filter(year %in% c(2010, 2013, 2016, 2017)) %>% # Get only before dredging (2010, 2013) and 2y postcon (2016, 2017)
dplyr::select(reef, dir, channel, site, site2, dist, year, date, # Reorder columns
transect, transect_area_m2, genus, species, diameter, count, perm)Coral size distributions: were small corals recorded?
ggplot(scler, aes(x = diameter)) +
geom_histogram(binwidth = 1) +
facet_wrap(~ year) +
scale_x_continuous(limits = c(0, 5))
Small corals (< 3cm) were recorded, but only in 2010, 2016, and 2017. Must filter out corals < 3 cm in order to compare data from all years.
Create subsetted dataset without small corals (< 3 cm)
## Many entries also do not include a size (diameter = NA)
## For coral density analysis, assume these are bigger than 3 cm and keep them
scler_3 <- scler %>%
filter(diameter >= 3 | is.na(diameter))Analyze coral density as a function of distance from channel
Decreases in coral density could have resulted from mortality due to disease, or mortality due to dredging sedimentation. Disease impacts should not show a relationship to distance from the dredging channel, but dredge impacts should. Therefore, analyzing changes in density as a function of distance from channel is one way to get at impacts due to dredging.
Change in density of larger corals (>3cm) – all pre-dredging vs. post-dredging count data
The full coral count dataset includes permanent sites and non-permanent sites, with data from 2010, 2013, 2015, and 2016/2017. Since corals <3 cm were not recorded in 2013 and 2015 surveys, we need to use the dataset excluding <3 cm corals.
# Sum number of large corals counted on each transect
large_counts <- scler_3 %>% # Use dataset that omits small corals
filter(dist < 1000) %>% # Use data within 1000m of the channel
group_by(reef, dir, channel, site, site2, dist, date, year, transect, transect_area_m2, perm) %>%
summarise(n = sum(count)) %>%
mutate(dens = n / transect_area_m2,
timepoint = cut(year, breaks = c(-Inf, 2014, Inf),
labels = c("pre-2014", "post-2014"))) %>%
ungroup()
# Plot density of corals for each transect
ggplot(large_counts, aes(x = dist, y = dens, color = timepoint)) +
facet_grid(dir ~ reef) +
geom_point() +
labs(y = "Corals / m2", x = "Distance from channel")
Reef 2 south post-dredging timepoint has an apparent outlier at 400m from channel. The transect was performed on 2017-03-29
# Subset data for only reef 2 south
r2s <- filter(large_counts, reef == "R2", dir == "S", timepoint == "post-2014")
# Fit a linear model with all points and calculate Cook's distance (cannot calculate Cook's for poisson glm)
mod2 <- lm(n/transect_area_m2 ~ dist, data = r2s)
nd2 <- crossing(dist = r2s$dist, transect_area_m2 = 1)
nd2$fit <- predict(mod2, newdata = nd2, re.form = NA, type = "response")
# Plot fitted lm
ggplot(r2s, aes(x = dist, y = dens, color = timepoint)) +
geom_point() +
geom_line(data = nd2, aes(x = dist, y = fit), inherit.aes = FALSE)
# Plot cook's distance ~ distance
plot(cooks.distance(mod2) ~ r2s$dist, xlab = "Distance from channel", ylab = "Cook's Distance")
text(r2s$dist - 30, cooks.distance(mod2),
labels = ifelse(cooks.distance(mod2) > 0.2, round(cooks.distance(mod2), 2), ""))
# Remove outlier point from R2S (it was the only data point from 2017-03-29)
large_counts <- filter(large_counts, date != "2017-03-29") Fit Poisson regression on count data
# Subset 1: locations with enough data at different distances from channel
large_counts_1 <- large_counts %>%
# only use where sufficient data exist to test dist
filter(reef == "R2" | (reef == "R3" & dir == "N") | (reef == "R3" & dir == "S" & timepoint == "pre-2014"))
# Fit model with distance from channel as predictor
mod_1 <- glmer(n ~ dist * reef:dir * timepoint + # count as response
offset(log(transect_area_m2)) + # offset accounts for diff transect sizes
(1|site/site2/transect) + (1|date), # random effects
data = large_counts_1, family = "poisson") # poisson model
# Get fitted responses
newdata_1 <- large_counts_1 %>%
tidyr::expand(nesting(reef, dir, timepoint), dist = 1:500,
transect_area_m2 = 1) # transect area of 1 m2 will predict count per m2
newdata_1$fit <- predict(mod_1, newdata_1, type = "response", re.form = NA)
# Get confidence intervals by simulation
set.seed(5234)
bootfit_1 <- bootMer(mod_1, FUN = function(x) predict(x, newdata_1, type = "response", re.form = NA),
nsim = 400, parallel = "multicore", ncpus = 8)
newdata_1$lci <- apply(bootfit_1$t, 2, quantile, 0.025)
newdata_1$uci <- apply(bootfit_1$t, 2, quantile, 0.975)
# Get fitted value of coral density adjacent to channel (i.e., dist = 20 m)
lsm20 <- lsmeans(mod_1, specs = c("reef", "dir", "timepoint"),
at = list(dist = 20, transect_area_m2 = 1), transform = "response")
# Get fitted value of coral density away from channel (i.e., dist = 300 m)
lsm300 <- lsmeans(mod_1, specs = c("reef", "dir", "timepoint"),
at = list(dist = 300, transect_area_m2 = 1), transform = "response")
# Gather fitted values for annotating plot
fit20 <- summary(lsm20) %>%
mutate(dist = 20) %>%
dplyr::select(reef, dir, timepoint, rate, dist) %>%
spread(timepoint, rate)
fit300 <- summary(lsm300) %>%
mutate(dist = 300) %>%
dplyr::select(reef, dir, timepoint, rate, dist) %>%
spread(timepoint, rate)
fits <- as_tibble(bind_rows(fit20, fit300))
# Test for significant differences by timepoint for fitted values
diff20 <- rbind(pairs(lsm20, by = c("reef", "dir")), adjust = "none")
diff300 <- rbind(pairs(lsm300, by = c("reef", "dir")), adjust = "none")
# Gather fitted values for annotating plot
sig20 <- summary(diff20) %>%
mutate(dist = 20) %>%
dplyr::select(reef, dir, estimate, p.value, dist)
sig300 <- summary(diff300) %>%
mutate(dist = 300) %>%
dplyr::select(reef, dir, estimate, p.value, dist)
sigs <- as_tibble(bind_rows(sig20, sig300))
# Generate data frame to use to annotate the plot
anno <- full_join(fits, sigs) %>%
mutate(rel = -(estimate) / `pre-2014` * 100,
rel2 = ifelse(p.value > 0.1, "", paste0(round(rel, 1), "%")),
sig = gtools::stars.pval(p.value),
sig2 = gsub(" ", "ns", sig))
# Subset 2: areas where data are only available close to channel
large_counts_2 <- large_counts %>%
filter(dist < 50) %>%
filter(reef == "HB" | (reef == "R3" & dir == "S")) %>% # & timepoint == "post-2014")) %>%
group_by(reef, dir) %>%
mutate(med = median(dist))
mod_2 <- glmer(n ~ reef * dir * timepoint + # count as response
offset(log(transect_area_m2)) + # offset accounts for diff transect sizes
(1|site/site2/transect) + (1|date), # random effects
data = large_counts_2, family = "poisson")
lsm <- lsmeans(mod_2, specs = c("reef", "dir", "timepoint"),
at = list(transect_area_m2 = 1), transform = "response")
sig <- rbind(pairs(lsm, by = c("reef", "dir")), adjust = "none")
sig <- summary(sig) %>%
mutate(dist = 25) %>%
dplyr::select(reef, dir, estimate, p.value, dist) %>% drop_na()
lsm <- summary(lsm, level = 0.84) %>%
mutate(dist = 25)
lsm2 <- lsm %>%
dplyr::select(reef, dir, timepoint, rate, dist) %>%
spread(timepoint, rate)
# Don't plot point for R3S pre, because already fitted with line
lsm <- lsm[-4,]
anno_2 <- full_join(lsm2, sig) %>%
mutate(rel = -(estimate) / `pre-2014` * 100,
rel2 = ifelse(p.value > 0.1, "", paste0(round(rel, 1), "%")),
sig = gtools::stars.pval(p.value),
sig2 = gsub(" ", "ns", sig))
anno <- bind_rows(anno, anno_2)
# Save as RData
save(newdata_1, lsm, anno, file = "data/processed/large.dens.mod.RData")
# PLOT ALL RESULTS FOR LARGE CORALS
ggplot(data = newdata_1, aes(x = dist, y = fit, color = timepoint)) +
facet_grid(dir ~ reef, scales = "free_x") +
geom_line() +
scale_x_continuous(limits = c(0, 500)) +
scale_y_continuous(limits = c(0.2, 1.6), expand = c(0,0)) +
geom_ribbon(aes(ymin = lci, ymax = uci, linetype = NA), alpha = 0.2) +
geom_segment(data = anno, inherit.aes = FALSE,
aes(x = dist, y = `pre-2014`, xend = dist, yend = `post-2014`),
arrow = arrow(length = unit(0.03, "npc"))) +
geom_text(data = anno, inherit.aes = FALSE, hjust = "left",
aes(x = dist, y = `post-2014` - 0.1, label = paste0(rel2, sig2))) +
labs(x = "Distance from channel", y = "# corals > 3cm per m2") +
geom_point(data = lsm, aes(x = dist, y = rate), alpha = 0.5) +
geom_errorbar(data = lsm, aes(x = dist, y = rate, ymin = asymp.LCL, ymax = asymp.UCL))
# Find distance at which difference is not significant
testlist <- map(1:500, function(x) data.frame(summary(rbind(pairs(lsmeans(mod_1, specs = c("reef", "dir", "timepoint"),
at = list(dist = x, transect_area_m2 = 1), transform = "response"), by = c("reef", "dir")), adjust = "none"))))
alltests <- bind_rows(testlist, .id = "dist")
alltests %>%
group_by(reef, dir) %>%
summarise(sig_dist = min(dist[!p.value < 0.05]))## # A tibble: 4 x 3
## # Groups: reef [?]
## reef dir sig_dist
## <chr> <chr> <chr>
## 1 R2 N 153
## 2 R2 S 1
## 3 R3 N 267
## 4 R3 S <NA>Change in density of small corals (1-2cm) (2010 vs. 2016-2017)
Transects at increasing distances from channel were recorded in 2010 baseline surveys and 2016/2017 cross-site surveys, when small corals were recorded.
# Subset 2010 baseline and 2016/2017 cross-site survey data <1000 m from channel
small_counts <- scler %>% # Select dataset that includes all corals, including small
filter(year %in% c(2010, 2016, 2017) & dist < 1000) %>%
filter(diameter < 3) %>% # ONLY 1-2cm CORALS
group_by(reef, dir, channel, site, site2, dist, date, year, transect, transect_area_m2, perm) %>%
summarise(n = sum(count)) %>%
mutate(dens = n / transect_area_m2,
timepoint = cut(year, breaks = c(-Inf, 2011, Inf),
labels = c("2010 baseline", "2016/2017 cross-site surveys"))) %>%
ungroup()
# Plot density of corals for each transect
ggplot(small_counts, aes(x = dist, y = dens, color = timepoint)) +
facet_grid(dir ~ reef) +
geom_point() +
labs(y = "Small corals / m2", x = "Distance from channel")
Fit poisson regression on count data
small_counts.f <- small_counts
# Fit model
mod_3 <- glmer(n ~ dist * reef * dir * timepoint + # count as response
offset(log(transect_area_m2)) + # offset accounts for diff transect sizes
(1|site/site2/transect) + (1|date), # random effects
data = small_counts.f, family = "poisson") # poisson model
# Get fitted responses
newdata_2 <- small_counts.f %>%
tidyr::expand(nesting(reef, dir, timepoint), dist = 1:500,
transect_area_m2 = 1) # transect area of 1 m2 will predict count per m2
newdata_2$fit <- predict(mod_3, newdata_2, type = "response", re.form = NA)
# Get confidence intervals by simulation
set.seed(9014)
bootfit_2 <- bootMer(mod_3, FUN = function(x) predict(x, newdata_2, type = "response", re.form = NA),
nsim = 400, parallel = "multicore", ncpus = 8)
newdata_2$lci <- apply(bootfit_2$t, 2, quantile, 0.025)
newdata_2$uci <- apply(bootfit_2$t, 2, quantile, 0.975)Plot fits and test for change in density near channel
# Subset areas where sufficient data exist to model based on distance from channel
newdata_2 <- filter(newdata_2, reef == "R2" | (reef == "R3" & dir == "N") | (reef == "R3" & dir == "S" & timepoint == "2010 baseline"))
# Plot raw data and fitted results
ggplot(newdata_2, aes(x = dist, y = fit, color = timepoint)) +
facet_grid(dir ~ reef) +
geom_ribbon(aes(ymin = lci, ymax = uci, linetype = NA), alpha = 0.2) +
geom_line() +
geom_point(data = small_counts, aes(y = dens)) +
labs(x = "Distance from channel", y = "# corals per m2")
# Get fitted value of coral density adjacent to channel (i.e., dist = 20 m)
lsm20 <- lsmeans(mod_3, specs = c("reef", "dir", "timepoint"),
at = list(dist = 20, transect_area_m2 = 1), transform = "response")
# Get fitted value of coral density away from channel (i.e., dist = 300 m)
lsm300 <- lsmeans(mod_3, specs = c("reef", "dir", "timepoint"),
at = list(dist = 300, transect_area_m2 = 1), transform = "response")
# Gather fitted values for annotating plot
fit20 <- summary(lsm20) %>%
mutate(dist = 20) %>%
select(reef, dir, timepoint, rate, dist) %>%
spread(timepoint, rate)
fit300 <- summary(lsm300) %>%
mutate(dist = 300) %>%
select(reef, dir, timepoint, rate, dist) %>%
spread(timepoint, rate)
fits <- as_tibble(bind_rows(fit20, fit300))
# Test for significant differences by timepoint for fitted values
diff20 <- rbind(pairs(lsm20, by = c("reef", "dir")), adjust = "none")
diff300 <- rbind(pairs(lsm300, by = c("reef", "dir")), adjust = "none")
# Gather fitted values for annotating plot
sig20 <- summary(diff20) %>%
mutate(dist = 20) %>%
select(reef, dir, estimate, p.value, dist)
sig300 <- summary(diff300) %>%
mutate(dist = 300) %>%
select(reef, dir, estimate, p.value, dist)
sigs <- as_tibble(bind_rows(sig20, sig300))
#anno <- full_join(fits, sigs)
anno <- full_join(fits, sigs) %>%
mutate(rel = -(estimate) / `2010 baseline` * 100,
rel2 = ifelse(p.value > 0.1, "", paste0(round(rel, 1), "%")),
sig = gtools::stars.pval(p.value),
sig2 = gsub(" ", "ns", sig))
anno <- anno %>%
mutate(`2010 baseline` = ifelse((reef == "HB" | (reef == "R3" & dir == "S")), NA, `2010 baseline`),
`2016/2017 cross-site surveys` = ifelse(reef == "HB", NA, `2016/2017 cross-site surveys`),
sig2 = ifelse((reef == "HB" | (reef == "R3" & dir == "S")), "", sig2))
# Subset 2: areas where data are only available close to channel
small_counts_2 <- small_counts %>%
filter(dist < 50) %>%
filter(reef == "HB" | (reef == "R3" & dir == "S")) %>%
group_by(reef, dir) %>%
mutate(med = median(dist))
ggplot(small_counts_2, aes(x = dist, y = dens, color = timepoint)) +
facet_grid(dir ~ reef) +
geom_point() +
labs(y = "Small corals / m2", x = "Distance from channel")
mod_2 <- glmer(n ~ reef * dir * timepoint + # count as response
offset(log(transect_area_m2)) + # offset accounts for diff transect sizes
(1|site/site2/transect) + (1|date), # random effects
data = small_counts_2, family = "poisson")
lsm <- lsmeans(mod_2, specs = c("reef", "dir", "timepoint"),
at = list(transect_area_m2 = 1), transform = "response")
sig <- rbind(pairs(lsm, by = c("reef", "dir")), adjust = "none")
sig <- summary(sig) %>%
mutate(dist = 25) %>%
dplyr::select(reef, dir, estimate, p.value, dist) %>% drop_na()
lsm <- summary(lsm, level = 0.84) %>%
mutate(dist = 25)
lsm2 <- lsm %>%
dplyr::select(reef, dir, timepoint, rate, dist) %>%
spread(timepoint, rate)
# Don't plot point for R3S pre, because already fitted with line
lsm <- lsm[-4,]
# Don't plot for HB before because no data
lsm <- lsm[-(c(1, 3)),]
anno_2 <- full_join(lsm2, sig) %>%
mutate(rel = -(estimate) / `2010 baseline` * 100,
rel2 = ifelse(p.value > 0.1, "", paste0(round(rel, 1), "%")),
sig = gtools::stars.pval(p.value),
sig2 = gsub(" ", "ns", sig))
anno <- bind_rows(drop_na(anno), anno_2)
anno <- drop_na(anno)
# PLOT ALL RESULTS FOR LARGE CORALS
ggplot(data = newdata_2, aes(x = dist, y = fit, color = timepoint)) +
facet_grid(dir ~ reef, scales = "free_x") +
geom_line() +
scale_x_continuous(limits = c(0, 500)) +
geom_ribbon(aes(ymin = lci, ymax = uci, linetype = NA), alpha = 0.2) +
geom_segment(data = anno, inherit.aes = FALSE,
aes(x = dist, y = `2010 baseline`, xend = dist, yend = `2016/2017 cross-site surveys`),
arrow = arrow(length = unit(0.03, "npc"))) +
geom_text(data = anno, inherit.aes = FALSE, hjust = "left",
aes(x = dist, y = `2016/2017 cross-site surveys` - 0.01, label = paste0(rel2, sig2))) +
labs(x = "Distance from channel", y = "# corals 1-2cm per m2") +
geom_point(data = lsm, aes(x = dist, y = rate), alpha = 0.5) +
geom_errorbar(data = lsm, aes(x = dist, y = rate, ymin = asymp.LCL, ymax = asymp.UCL)) +
coord_cartesian(ylim = c(0, 1.5))
save(newdata_2, lsm, anno, file = "data/processed/small.dens.mod.RData")
# Find distance at which difference is not significant
testlist <- map(1:500, function(x) data.frame(summary(rbind(pairs(lsmeans(mod_3, specs = c("reef", "dir", "timepoint"),
at = list(dist = x, transect_area_m2 = 1), transform = "response"), by = c("reef", "dir")), adjust = "none"))))
alltests <- bind_rows(testlist, .id = "dist")
alltests %>%
group_by(reef, dir) %>%
summarise(sig_dist = min(dist[!p.value < 0.05]))## # A tibble: 6 x 3
## # Groups: reef [?]
## reef dir sig_dist
## <chr> <chr> <chr>
## 1 HB N <NA>
## 2 HB S <NA>
## 3 R2 N 202
## 4 R2 S 1
## 5 R3 N 1
## 6 R3 S 1Analyze coral size frequency distributions
Look at tiny corals (< 3cm) abundance before and after dredging. Corals <3cm in diameter were only recorded in 2010 (within 500 m of channel) and in 2016/2017 at cross-sites and permanent monitoring sites. The surveys that are most comparable are those within 100m of channel (see DCA cross site report page 2). Unfortunately, there are no “control” locations at which we can compare size frequency distributions. We currently don’t have HB data from 2010, so we can only compare R2 and R3 locations.
### Look at map on page 2 of october 2017 cross_impact report
##### Reef 2 N:: 2010 transects below 100 m are all same general area as R2N1
##### Reef 2 S:: 2010 transects below 100 m are all same general area as R2S1
##### Reef 3 N:: 2010 transects below 100 m are all same general area as R3N1
##### Reef 3 S:: 2010 transects LR are all same general area as R3S2
# Filter: 2010 and 2016/2017, only dist < 100 m from channel
df <- filter(scler, year %in% c(2010, 2016, 2017),
dist <= 100,
site %in% c("R2N", "R2S", "R3N", "R3S", "R2N1", "R2S1", "R3N1", "R3S2")) %>%
mutate(timepoint = cut(year, breaks = c(-Inf, 2011, Inf),
labels = c("2010 baseline", "2016/2017 cross-site surveys")))
# Kruskal-Wallis non-parametric test for difference in mean
kwtest <- df %>%
nest(-reef, -dir) %>%
mutate(kw = map(data, ~ kruskal.test(diameter ~ timepoint, data = .)),
p = map_dbl(kw, ~ .$p.value))
kwtest <- kwtest %>% mutate(x = 15, y = 0.1)
save("df", "kwtest", file = "data/processed/sizefreq.RData")
# Plot size frequency distributions and annotate with KW test p values
ggplot(df, aes(x = diameter, fill = year <= 2010)) +
geom_histogram(aes(y = ..density..), # Scales each histogram so the total area = 1
position = "identity", binwidth = 1,
alpha = 0.3) +
facet_grid(dir ~ reef) +
labs(title = paste("Coral size frequency within", max(df$dist), "m of channel"),
x = "Coral diameter (cm)", y = "Proportion of individuals") +
scale_fill_manual(labels = c("after dredging", "before dredging"), values = c("#F8766D", "#00BFC4")) +
theme(legend.title = element_blank(),
legend.position = c(0.3, 0.9),
legend.background = element_blank()) +
xlim(0, 20) +
geom_text(data = kwtest, aes(fill = NULL, x = x, y = y,
label = paste("p =", formatC(p, format = "e", digits = 2))))
df %>%
group_by(reef, dir, timepoint) %>%
summarise(med = median(diameter),
mean = mean(diameter))## # A tibble: 8 x 5
## # Groups: reef, dir [?]
## reef dir timepoint med mean
## <chr> <chr> <fct> <dbl> <dbl>
## 1 R2 N 2010 baseline 3 5.09
## 2 R2 N 2016/2017 cross-site surveys 5 8.06
## 3 R2 S 2010 baseline 3 4.95
## 4 R2 S 2016/2017 cross-site surveys 5 7.68
## 5 R3 N 2010 baseline 8 10.5
## 6 R3 N 2016/2017 cross-site surveys 9 9.10
## 7 R3 S 2010 baseline 4 5.52
## 8 R3 S 2016/2017 cross-site surveys 5 6.51Estimate number of corals lost based on reef area data
# Import reef area data (from Maria's GIS analysis)
HBN <- readxl::read_xlsx(path = "data/reef_area/Distance from Channel by Reef Location.xlsx", sheet = "Inner North")
HBS <- readxl::read_xlsx(path = "data/reef_area/Distance from Channel by Reef Location.xlsx", sheet = "Inner South")
R2N <- readxl::read_xlsx(path = "data/reef_area/Distance from Channel by Reef Location.xlsx", sheet = "Middle North")
R2S <- readxl::read_xlsx(path = "data/reef_area/Distance from Channel by Reef Location.xlsx", sheet = "Middle South")
R3N <- readxl::read_xlsx(path = "data/reef_area/Distance from Channel by Reef Location.xlsx", sheet = "Outer North")
R3S <- readxl::read_xlsx(path = "data/reef_area/Distance from Channel by Reef Location.xlsx", sheet = "Outer South")
# Tidy reef area data
reef_areas <- bind_rows(HBN = HBN, HBS = HBS, R2N = R2N, R2S = R2S, R3N = R3N, R3S = R3S, .id = "reefdir") %>%
clean_names() %>%
rename(dist = distance_from_channel) %>%
group_by(reefdir) %>%
mutate(area = area_in_sq_m - lag(area_in_sq_m, default = 0),
reef = str_sub(reefdir, 1, 2),
dir = str_sub(reefdir, 3, 3)) %>%
ungroup() %>%
select(reef, dir, dist, area)
# Load change in large coral density data
load("data/processed/large.dens.mod.RData") # newdata_1 object contains fitted values for density before and after dredging
total_loss_large <- newdata_1 %>%
select(reef, dir, timepoint, dist, fit) %>%
spread(timepoint, fit) %>%
mutate(dens_change = `post-2014` - `pre-2014`) %>%
filter(dens_change < 0) %>%
left_join(reef_areas) %>%
mutate(abs_loss_per_m = dens_change * area) %>%
group_by(reef, dir) %>%
summarise(total_loss = sum(abs_loss_per_m),
out_to = max(dist))
total_loss_large## # A tibble: 3 x 4
## # Groups: reef [?]
## reef dir total_loss out_to
## <chr> <chr> <dbl> <dbl>
## 1 R2 N -72982. 349
## 2 R2 S -16912. 192
## 3 R3 N -58636. 500# Load change in small coral density data
load("data/processed/small.dens.mod.RData") # newdata_2 object contains fitted values for density before and after dredging
total_loss_small <- newdata_2 %>%
select(reef, dir, timepoint, dist, fit) %>%
spread(timepoint, fit) %>%
mutate(dens_change = `2016/2017 cross-site surveys` - `2010 baseline`) %>%
filter(dens_change < 0) %>%
left_join(reef_areas) %>%
mutate(abs_loss_per_m = dens_change * area) %>%
group_by(reef, dir) %>%
summarise(total_loss = sum(abs_loss_per_m),
out_to = max(dist))
total_loss_small## # A tibble: 3 x 4
## # Groups: reef [?]
## reef dir total_loss out_to
## <chr> <chr> <dbl> <dbl>
## 1 R2 N -203249. 500
## 2 R2 S -197917. 500
## 3 R3 N -11082. 500total_loss <- bind_rows(small = total_loss_small, large = total_loss_large, .id = "sizeclass")
total_loss## # A tibble: 6 x 5
## # Groups: reef [2]
## sizeclass reef dir total_loss out_to
## <chr> <chr> <chr> <dbl> <dbl>
## 1 small R2 N -203249. 500
## 2 small R2 S -197917. 500
## 3 small R3 N -11082. 500
## 4 large R2 N -72982. 349
## 5 large R2 S -16912. 192
## 6 large R3 N -58636. 500# Add info for total area within zone of loss
total_loss <- reef_areas %>%
mutate(total_area = cumsum(area)) %>% # calculate cumulative area at each distance
rename(out_to = dist) %>% # combine with total_loss df
right_join(total_loss) %>%
mutate(total_area = total_area / 1000000) # convert m2 to km2
total_loss_f <- total_loss %>%
ungroup() %>%
mutate(reef = recode(reef, R2 = "Inner", R3 = "Outer"),
dir = recode(dir, N = "Northern", S = "Southern"),
sizeclass = recode(sizeclass, small = "1-2 cm", large = "≥ 3 cm"),
total_loss = round(total_loss, 0)) %>%
select(dir, reef, sizeclass, total_loss, out_to, total_area) %>%
arrange(reef, dir)
write_tsv(total_loss_f, path = "figures/Table1.txt")
sum(total_loss_f$total_loss)## [1] -560778