.Rhistory

geom_smooth(aes(y=mean_scars_m2_coral, x=mean_pct_coral_cover, color=region), method="lm")+
geom_jitter(aes(y=mean_scars_m2_coral, x=mean_pct_coral_cover, color=region), alpha=0.7)+
scale_colour_manual(values=region_palette)+
scale_y_continuous(trans = scales::log_trans(),
breaks= c(0.05, 0.5, 5, 50, 500),
labels= c(0.05, 0.5, 5, 50, 500),
limits=c(0.04, 700))+
scale_x_continuous(expand = c(0.05, 0))+ #values near 0 cutoff, expands x-axis to resolve issue
labs(y=expression(paste("Scar density (n ", m^-2, " coral)")), x="Coral cover (%)", color="Region")+
theme_bw()+
theme(legend.position = "top")+
rremove("x.title") +
rremove("x.ticks") +
rremove("x.text")+
rremove("y.title") +
rremove("y.ticks") +
rremove("y.text")
# coral area grazed ~ coral cover
site_level_coral_area_grazed_by_coral_cover <- fish_corals_site_level %>%
ggplot()+
geom_smooth(aes(y=mean_pct_coral_grazed, x=mean_pct_coral_cover, color=region), method="lm")+
geom_jitter(aes(y=mean_pct_coral_grazed, x=mean_pct_coral_cover, color=region), alpha=0.7)+
scale_colour_manual(values=region_palette)+
scale_y_continuous(trans = scales::log_trans(),
breaks= c(0.00025, 0.0025, 0.025, 0.25, 2.5),
labels= c("0.00025", 0.0025, 0.025, 0.25, 2.5),
limits = c(0.00007,3))+
scale_x_continuous(expand = c(0.05, 0))+ #values near 0 cutoff, expands x-axis to resolve issue
labs(y="Coral area preyed upon (%)", x="Coral cover (%)", color="Region")+
theme_bw()+
theme(legend.position = "top")+
rremove("y.title") +
rremove("y.ticks") +
rremove("y.text")
# making four panel plot in two parts to overcome weird spacing issues when making a four panel plot directly
site_level_corallivory_patterns_ac <- ggpubr::ggarrange(site_level_scar_density_by_corallivore_density +
theme(plot.margin =unit(c(t = 5, r = 5, b = 5, l = 25), "pt"))+
rremove("legend"),
site_level_coral_area_grazed_by_corallivore_density+
theme(plot.margin =unit(c(t = 5, r = 5, b = 5, l = 25), "pt"),
axis.title.x = element_text(margin = margin(t = 2.5, r = 0, b = 0, l = 0)))+
rremove("legend"), labels=c("a", "c"),
nrow=2, ncol=1, align="v", heights=c(1, 1.15))
site_level_corallivory_patterns_bd <- ggpubr::ggarrange(site_level_scar_density_by_coral_cover+
theme(plot.margin =unit(c(t = 0, r = 5, b = 5, l = 25), "pt"))+
rremove("legend"),
site_level_coral_area_grazed_by_coral_cover+
theme(plot.margin =unit(c(t = 5, r = 5, b = 0, l = 25), "pt"),
axis.title.x = element_text(
margin = margin(t = 7.5, r = 0, b = 0, l = 0)))+
rremove("legend"), labels=c("b", "d"),
nrow=2, ncol=1, align="v", heights=c(1, 1.15))
site_level_corallivory_patterns_temp <- ggpubr::ggarrange(site_level_corallivory_patterns_ac,
site_level_corallivory_patterns_bd+
theme(plot.margin =unit(c(5, 5, 5, 5), "pt")),
widths=c(1.15, 1))
common_legend <- ggpubr::get_legend(site_level_scar_density_by_coral_cover) %>%
as_ggplot()
site_level_corallivory_patterns <- ggpubr::ggarrange(common_legend,
site_level_corallivory_patterns_temp,
nrow=2, ncol=1, heights = c(1,10))
# image/icon placement will look off in the chunk output but correct in the saved image
site_level_corallivory_patterns_with_images <- cowplot::ggdraw() +
draw_plot(site_level_corallivory_patterns)+
draw_image(major_corallivore_icon,  x = -0.38, y = -0.47, scale=0.08) + #adds the major pf icon
draw_image(coral_icon, x = 0.175, y = -0.47, scale=0.055)+ #adds the coral
draw_image(corallivory_icon,  x = -0.448, y = 0.379, scale=0.055) + #adds corallivory icon
draw_image(corallivory_icon,  x = -0.448, y = -0.04, scale=0.055) #adds corallivory icon
site_level_corallivory_patterns_with_images
#saving as .jpg version (for presentations) and .pdf version (for publications)
ggsave(plot=site_level_corallivory_patterns_with_images,
filename=here("figures_tables/fig3_corallivory_by_coral_and_fish.jpg"),
width = 215, height = 160, dpi = 600, units = "mm")
ggsave(plot=site_level_corallivory_patterns_with_images,
filename=here("figures_tables/fig3_corallivory_by_coral_and_fish.pdf"),
width = 215, height = 160, dpi = 600, units = "mm")
# scar density ~ corallivore biomass
site_level_scar_density_by_corallivore_biomass <-
fish_corals_site_level %>%
ggplot()+
geom_smooth(aes(x=mean_pf_g_100m2, y=mean_scars_m2_coral, color=region),method="lm")+
geom_point(aes(x=mean_pf_g_100m2, y=mean_scars_m2_coral, color=region), alpha=0.7)+
scale_color_manual(values=region_palette)+
scale_y_continuous(trans = scales::log_trans(),
breaks= c(0.05, 0.5, 5, 50, 500),
labels= c(0.05, 0.5, 5, 50, 500))+
scale_x_continuous(expand = c(0.1, 0),
breaks= c(0, 500, 1000,1500, 2000, 2500),
labels= c(0, 500, 1000, 1500, 2000, 2500))+ #values near 0 cutoff, expands x-axis to resolve issued
labs(x=expression(paste("Major corallivore biomass (g 100 ", m^-2, ")")),
y=expression(paste("Scar density (n ", m^-2, " coral)")),
color="Region")  +
theme_bw()+
theme(legend.position = "top")+
rremove("x.title") +
rremove("x.ticks") +
rremove("x.text")+
rremove("y.title") +
rremove("y.ticks") +
rremove("y.text")
# coral area grazed ~ corallivore biomass
site_level_coral_area_grazed_by_corallivore_biomass <-
fish_corals_site_level %>%
ggplot()+
geom_smooth(aes(x=mean_pf_g_100m2, y=mean_pct_coral_grazed, color=region), method="lm")+
geom_point(aes(x=mean_pf_g_100m2, y=mean_pct_coral_grazed, color=region), alpha=0.7)+
scale_color_manual(values=region_palette)+
scale_y_continuous(trans = scales::log_trans(),
breaks= c(0.00025, 0.0025, 0.025, 0.25, 2.5),
labels= c("0.00025", 0.0025, 0.025, 0.25, 2.5))+
scale_x_continuous(expand = c(0.1, 0),
breaks= c(0, 500, 1000,1500, 2000, 2500),
labels= c(0, 500, 1000, 1500, 2000, 2500))+ #values near 0 cutoff, expands x-axis to resolve issue
labs(x=expression(paste("Major corallivore biomass (g 100 ", m^-2, ")")),
y=expression(paste("Coral area preyed upon (%)")),
color="Region") +
theme_bw()+
theme(legend.position = "top")+
rremove("y.title") +
rremove("y.ticks") +
rremove("y.text")
# making combo figure
fig_corallivory_by_pf_biomass_site <- ggarrange(site_level_scar_density_by_corallivore_biomass,
site_level_coral_area_grazed_by_corallivore_biomass,
nrow=2, ncol=1, legend="none", labels=c("b", "d"))
site_level_fish_density_biomass_temp <- ggpubr::ggarrange(site_level_corallivory_patterns_ac,
fig_corallivory_by_pf_biomass_site+
theme(plot.margin =unit(c(5, 5, 5, 5), "pt")),
widths=c(1.15, 1))
common_legend <- ggpubr::get_legend(site_level_scar_density_by_coral_cover) %>%
as_ggplot()
site_level_fish_density_biomass <- ggpubr::ggarrange(common_legend,
site_level_fish_density_biomass_temp,
nrow=2,
ncol=1,
heights = c(1,10))
# making four panel plot in two parts to overcome weird spacing issues when making a four panel plot directly
site_level_corallivory_patterns_ac <- ggpubr::ggarrange(site_level_scar_density_by_corallivore_density +
theme(plot.margin =unit(c(l=5, r=5, t=5, b=0), "pt"))+
rremove("legend"),
site_level_scar_density_by_corallivore_biomass+
theme(plot.margin =unit(c(l=5, r=5, t=5, b=5), "pt"))+
rremove("legend")+
rremove("ylab"),
site_level_coral_area_grazed_by_corallivore_density+
theme(plot.margin =unit(c(l=5, r=5, t=5, b=0), "pt"))+
rremove("legend"),
site_level_coral_area_grazed_by_corallivore_biomass+
theme(plot.margin =unit(c(l=5, r=5, t=5, b=5), "pt"))+
rremove("legend")+
rremove("ylab"),
nrow=2,
ncol=2,
labels=c("a","b","c","d"),
heights=c(1, 1.15, 1, 1.15),
common.legend = TRUE)
# making four panel plot in two parts to overcome weird spacing issues when making a four panel plot directly
site_level_corallivory_patterns_supp_bd <- ggpubr::ggarrange(site_level_scar_density_by_corallivore_biomass+
theme(plot.margin =unit(c(t = 0, r = 5, b = 5, l = 25), "pt"))+
rremove("legend"),
site_level_coral_area_grazed_by_corallivore_biomass+
theme(plot.margin =unit(c(t = 5, r = 5, b = 0, l = 25), "pt"),
axis.title.x = element_text(
margin = margin(t = 7.5, r = 0, b = 0, l = 0)))+
rremove("legend"),
nrow=2,
ncol=1,
align="v",
labels=c("b", "d"),
heights=c(1, 1.15))
site_level_corallivory_patterns_supp_temp <- ggpubr::ggarrange(site_level_corallivory_patterns_ac, #fig a-c made in previous chunk
site_level_corallivory_patterns_supp_bd+
theme(plot.margin =unit(c(5, 5, 5, 5), "pt")),
widths=c(1.15, 1))
site_level_corallivory_patterns_supp <- ggpubr::ggarrange(common_legend, #common legend object made in previous chunk
site_level_corallivory_patterns_temp,
nrow=2,
ncol=1,
heights = c(1,10))
ggsave(plot=site_level_corallivory_patterns_supp,
filename=here("figures_tables/supplements_figs/fig_corallivory_by_pf_biomass_density_site.jpg"),
width = 215, height = 160, dpi = 600, units = "mm")
set.seed(1945) #for reproducible results from random simulations of residuals
## random-slope and intercept model with interaction: warnings - model is not identifiable/won't converge (for both gaussian and t distributions)
m1_scar_pf_density_gauss <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_n_100m2*mean_pct_coral_cover + (mean_pf_n_100m2 + mean_pct_coral_cover | region), data = fish_corals_site_level, family=gaussian)
m1_scar_pf_density <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_n_100m2*mean_pct_coral_cover + (mean_pf_n_100m2 + mean_pct_coral_cover | region), data = fish_corals_site_level, family=t_family)
# random-intercept model model with interaction (no explanatory variable transformation)
# Gaussian model
m2_scar_pf_density_gauss <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_n_100m2*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level, family=gaussian)
m2_scar_pf_density_residuals_gauss <- DHARMa::simulateResiduals(m2_scar_pf_density_gauss, n = 1000, use.u = T) #evaluating model fit
plot(m2_scar_pf_density_residuals_gauss) #odd residual versus predicted pattern, trying transformation
# t-family model
m2_scar_pf_density <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_n_100m2*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level, family=t_family)
m2_scar_pf_density_residuals <- DHARMa::simulateResiduals(m2_scar_pf_density, n = 1000, use.u = T) #evaluating model fit
plot(m2_scar_pf_density_residuals) #odd residual versus predicted pattern, trying transformation. But, looks better than Gaussian.
# random-intercept model model with interaction
m3_scar_pf_density_gauss <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ sqrt(mean_pf_n_100m2)*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level, family=gaussian)
m3_scar_pf_density <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ sqrt(mean_pf_n_100m2)*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level, family=t_family)
m3_scar_pf_density_residuals <- DHARMa::simulateResiduals(m3_scar_pf_density, n = 1000, use.u = T) #evaluating model fit
plot(m3_scar_pf_density_residuals) #looks better
#looking at residuals by region - nothing seems very extreme
bind_cols(m3_scar_pf_density_residuals$fittedResiduals, fish_corals_site_level$region) %>%
rename(y = `...1`,
region = `...2`) %>%
ggplot(aes(x = 1:26, y = y, color = region)) +
geom_point()
m3_scar_pf_density %>% summary()
# random-intercept model model without interaction
m4_scar_pf_density <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ sqrt(mean_pf_n_100m2)+mean_pct_coral_cover + (1 | region), data = fish_corals_site_level, family=t_family)
m4_scar_pf_density_residuals <- DHARMa::simulateResiduals(m4_scar_pf_density, n = 1000, use.u = T) #evaluating model fit
plot(m4_scar_pf_density_residuals)
m4_scar_pf_density %>% summary() #model results similar to previous model
#tidying model output and calculating confidence intervals
m3_scar_pf_density_ci <- confint(m3_scar_pf_density) %>% as_tibble() %>% select(-Estimate)
scar_pf_density_model_output <- m3_scar_pf_density %>%
broom.mixed::tidy() %>%
filter(term != "sd__Observation") %>%
bind_cols(m3_scar_pf_density_ci) %>%
filter(effect!="ran_pars") %>%
janitor::clean_names() %>%
rename(ci_2.5 = x2_5_percent, ci_97.5=x97_5_percent) %>%
mutate(model="M1 - Scar density") %>%
select(model, term, statistic, p_value, estimate, ci_2.5, ci_97.5) %>%
mutate(across(where(is.double), ~round(.x, 3)))
scar_pf_density_model_output
set.seed(1945) #for reproducible results from random simulations of residuals
##Random-slope and intercept model with interaction: warnings - model is not identifiable/won't converge
#m1_pctgrazed_pfdensity <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ mean_pf_n_100m2*mean_pct_coral_cover + (mean_pf_n_100m2 + mean_pct_coral_cover | region), data = fish_corals_site_level, family=gaussian)
#Random-slope model with interaction
m2_pctgrazed_pfdensity <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ mean_pf_n_100m2*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level,family=gaussian)
m2_pctgrazed_pfdensity_residuals <- DHARMa::simulateResiduals(m2_pctgrazed_pfdensity, n = 1000, use.u = T) #evaluating model fit
plot(m2_pctgrazed_pfdensity_residuals) #KS test is significant
#looking at distributions of predictors
hist(fish_corals_site_level$mean_pf_g_100m2) #heavily skewed
hist(sqrt(fish_corals_site_level$mean_pf_g_100m2)) #more normal w/ slight skew
hist(fish_corals_site_level$mean_pct_coral_cover) #skewed, but not as much as parrotfish density
hist(asin(sqrt(fish_corals_site_level$mean_pct_coral_cover/100))) #more normal w/ some kurtosis
#Significant quantile deviations in previous model, first using sqrt transformation on pf density only to try to address this
m3_pctgrazed_pfdensity <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ sqrt(mean_pf_n_100m2)*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level,family=gaussian)
m3_pctgrazed_pfdensity_residuals <- DHARMa::simulateResiduals(m3_pctgrazed_pfdensity, n = 1000, use.u = T) #evaluating model fit
plot(m3_pctgrazed_pfdensity_residuals) #QQ plot looks better and no sig issues
#looking at residuals by region
#2 slightly more extreme points (1 in STX and 1 in PAN), but no consistent differences by region
bind_cols(m3_pctgrazed_pfdensity_residuals$fittedResiduals, fish_corals_site_level$region) %>%
rename(y = `...1`,
region = `...2`) %>%
ggplot(aes(x = 1:26, y = y, color = region)) +
geom_point()
m3_pctgrazed_pfdensity %>% summary() #viewing model results
#since interaction is NS, confirming that observed patterns are similar when interaction term is dropped
m4_pctgrazed_pfdensity <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ sqrt(mean_pf_n_100m2)+mean_pct_coral_cover + (1 | region), data = fish_corals_site_level,family=gaussian)
m4_pctgrazed_pfdensity %>% summary() #similar results when interaction is dropped
#tidying model output and calculating confidence intervals (of model 3, with the interaction)
pctgrazed_pfdensity_m3_ci <- confint(m3_pctgrazed_pfdensity) %>% as_tibble() %>% select(-Estimate)
proportion_grazed_pfdensity_model_output <- m3_pctgrazed_pfdensity %>%
broom.mixed::tidy() %>%
filter(term != "sd__Observation") %>%
bind_cols(pctgrazed_pfdensity_m3_ci) %>%
filter(effect!="ran_pars") %>%
janitor::clean_names() %>%
rename(ci_2.5 = x2_5_percent, ci_97.5=x97_5_percent) %>%
mutate(model="M3 - Coral area preyed upon (%)") %>%
select(model, term, statistic, p_value, estimate, ci_2.5, ci_97.5) %>%
mutate(across(where(is.double), ~round(.x, 3)))
set.seed(1945) #for reproducible results from random simulations of residuals
#random-slope and intercept model with interaction: warnings - model is not identifiable/won't converge (for both gaussian and t distribution)
m1_scar_density_biomass_gauss <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_g_100m2*mean_pct_coral_cover + (mean_pf_n_100m2 + mean_pct_coral_cover | region), data = fish_corals_site_level, family=gaussian)
m1_scar_density_biomass <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_g_100m2*mean_pct_coral_cover + (mean_pf_n_100m2 + mean_pct_coral_cover | region), data = fish_corals_site_level, family=t_family)
#random-intercept model
# Gaussian model
m2_scar_density_biomass_gauss <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_g_100m2*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level, family=t_family)
scar_density_biomass_m2_residuals_gauss <- DHARMa::simulateResiduals(m2_scar_density_biomass_gauss, n = 1000, use.u = T)
plot(scar_density_biomass_m2_residuals_gauss)
# t-family model
m2_scar_density_biomass <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~ mean_pf_g_100m2*mean_pct_coral_cover + (1 | region), data = fish_corals_site_level, family=t_family)
scar_density_biomass_m2_residuals <- DHARMa::simulateResiduals(m2_scar_density_biomass, n = 1000, use.u = T)
plot(scar_density_biomass_m2_residuals) #residuals look better than Gaussian (lines more horizontal)
m2_scar_density_biomass %>% summary() #viewing model results
hist(fish_corals_site_level$mean_pf_g_100m2)
hist(sqrt(fish_corals_site_level$mean_pf_g_100m2)) #not perfect, but less skewness with sqrt transformation
#rerunning model with square-root transformation of corallivore density to attempt to address KS test issues
#using transformation to make it consistent and more comparable to other models
m3_scar_density_biomass <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~sqrt(mean_pf_g_100m2)*(mean_pct_coral_cover)+ (1 | region), data = fish_corals_site_level, family=t_family)
#evaluating model fit
m3_scar_density_biomass_residuals <- DHARMa::simulateResiduals(m3_scar_density_biomass, n = 1000, use.u = T)
plot(m3_scar_density_biomass_residuals) #QQ plot looks good but quantile deviations detected, taking a closer look at residuals
#looking at residuals by region - nothing seems very extreme
bind_cols(m3_scar_density_biomass_residuals$fittedResiduals, fish_corals_site_level$region) %>%
rename(y = `...1`,
region = `...2`) %>%
ggplot(aes(x = 1:26, y = y, color = region)) +
geom_point()
m3_scar_density_biomass %>% summary() #viewing model results
#since interaction is NS, confirming that observed patterns are similar when interaction term is dropped
m4_scar_density_biomass <- glmmTMB::glmmTMB(log(mean_scars_m2_coral) ~sqrt(mean_pf_g_100m2)+(mean_pct_coral_cover)+ (1 | region), data = fish_corals_site_level, family=t_family)
m4_scar_density_biomass %>% summary() #similar results when interaction is dropped
#tidying model output and calculating confidence intervals
m3_scar_density_biomass_ci <- confint(m3_scar_density_biomass) %>% as_tibble() %>% select(-Estimate)
scar_density_biomass_model_output <- m3_scar_density_biomass %>%
broom.mixed::tidy() %>%
filter(term != "sd__Observation") %>%
bind_cols(m3_scar_density_biomass_ci) %>%
filter(effect!="ran_pars") %>%
janitor::clean_names() %>%
rename(ci_2.5 = x2_5_percent, ci_97.5=x97_5_percent) %>%
mutate(model="M2 - Scar density") %>%
select(model, term, statistic, p_value, estimate, ci_2.5, ci_97.5) %>%
mutate(across(where(is.double), ~round(.x, 3)))
scar_density_biomass_model_output
set.seed(1945) #for reproducible results from random simulations of residuals
##Random-slope and intercept model with interaction: warnings - model is not identifiable/won't converge
m1_pct_grazed_biomass <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ mean_pf_g_100m2*mean_pct_coral_cover + (mean_pf_g_100m2 + mean_pct_coral_cover | region), data = fish_corals_site_level, family=gaussian)
#Random-slope model with interaction
m2_pct_grazed_biomass <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ mean_pf_g_100m2*mean_pct_coral_cover +(1 | region), data = fish_corals_site_level, family=gaussian)
m2_pct_grazed_biomass_residuals <- DHARMa::simulateResiduals(m2_pct_grazed_biomass, n = 1000, use.u = T)
plot(m2_pct_grazed_biomass_residuals) #significant quantile deviations, QQ plot shows variation in center not tails
#Significant quantile deviations in previous model, using transformation of predictor to address this
m3_pct_grazed_biomass <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ sqrt(mean_pf_g_100m2)*mean_pct_coral_cover +(1 | region), data = fish_corals_site_level, family=gaussian)
m3_pct_grazed_biomass_residuals <- DHARMa::simulateResiduals(m3_pct_grazed_biomass, n = 1000, use.u = T)
plot(m3_pct_grazed_biomass_residuals) #still some quantile deviations but adjusted test NS, diagnostic plots look much better
#looking at residuals by region
#same as previous: 2 more extreme points (1 in STX and 1 in PAN), but no consistent differences by region
bind_cols(m3_pct_grazed_biomass_residuals$fittedResiduals, fish_corals_site_level$region) %>%
rename(y = `...1`,
region = `...2`) %>%
ggplot(aes(x = 1:26, y = y, color = region)) +
geom_point()
m4_pct_grazed_biomass <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~ sqrt(mean_pf_g_100m2)+mean_pct_coral_cover +(1 | region), data = fish_corals_site_level, family=gaussian)
m4_pct_grazed_biomass_residuals <- DHARMa::simulateResiduals(m4_pct_grazed_biomass, n = 1000, use.u = T)
plot(m4_pct_grazed_biomass_residuals) #still some quantile deviations but adjusted test NS, diagnostic plots look much better
m4_pct_grazed_biomass %>% summary() #same overall sig as findings to model w/ interaction, just different effect size
#viewing model results
m3_pct_grazed_biomass %>% summary()
#since interaction is NS, confirming that observed patterns are similar when interaction term is dropped
m4_pct_grazed_biomass <- glmmTMB::glmmTMB(log(mean_pct_coral_grazed) ~sqrt(mean_pf_g_100m2)+(mean_pct_coral_cover)+ (1 | region), data = fish_corals_site_level, family=gaussian)
m4_pct_grazed_biomass %>% summary() #similar results when interaction is dropped
#tidying model output and calculating confidence intervals
m3_pct_grazed_biomass_ci <- confint(m3_pct_grazed_biomass) %>% as_tibble() %>% select(-Estimate)
proportion_grazed_biomass_model_output <- m3_pct_grazed_biomass %>%
broom.mixed::tidy() %>%
filter(term != "sd__Observation") %>%
bind_cols(m3_pct_grazed_biomass_ci) %>%
filter(effect!="ran_pars") %>%
janitor::clean_names() %>%
rename(ci_2.5 = x2_5_percent, ci_97.5=x97_5_percent) %>%
mutate(model="M4 - Coral area preyed upon (%)") %>%
select(model, term, statistic, p_value, estimate, ci_2.5, ci_97.5) %>%
mutate(across(where(is.double), ~round(.x, 3)))
table_s4 <- scar_pf_density_model_output %>%
bind_rows(proportion_grazed_pfdensity_model_output) %>%
bind_rows(scar_density_biomass_model_output) %>%
bind_rows(proportion_grazed_biomass_model_output)
write_csv(table_s4, here("figures_tables/supplements_tables/table_s4_site_level_glmm_results.csv"))
set.seed(1945) #for reproducible results from random simulations of residuals
# total abundance of colonies per transect
total_coral_abundance <- coral_scar_df %>%
group_by(region, site, transect, depth_m) %>%
reframe(total_n_colonies=n())
# calculating transect-level coral diversity and coral cover
coral_taxa_abundance_cover <- coral_scar_df %>%
group_by(region, site, transect, depth_m, coral_species) %>%
reframe(n_colonies=n(), #transect-level abundance of coral colonies per species
pct_coral_cover_species=(sum(col_area_cm2)/300000)*100) %>% #coral cover by species
left_join(total_coral_abundance) #total coral abundance per transect
# calculating coral diversity and total coral cover per transect
coral_diversity_and_cover <- coral_taxa_abundance_cover %>%
mutate(prop_abundance=n_colonies/total_n_colonies, #proportional abundance of coral taxa
p_ln_prop_abundance=prop_abundance*log(prop_abundance)) %>% #for diversity calculations
group_by(region, site, transect, depth_m) %>%
reframe(coral_diversity=sum(p_ln_prop_abundance)*-1, #summarizing diversity per transect
pct_coral_cover=sum(pct_coral_cover_species)) #summarizing transect-level coral cover
# abundance of colonies per species per transect
n_colonies_species_transect <- coral_scar_df %>%
group_by(region, site, transect, depth_m, coral_species) %>%
dplyr::reframe(n_colonies=n()) %>% #transect-level abundance of coral colonies per species
left_join(total_coral_abundance, by = join_by(region, site, transect, depth_m))
# transect-level summary of the number of grazed colonies of four heavily predated taxa
transect_level_corallivory <- coral_scar_df %>%
group_by(region, site, transect, depth_m, coral_species) %>%
filter(n_new_scars>0) %>%
dplyr::reframe(n_colonies_grazed=n())
transect_corallivory_diversity_df <- coral_taxa_abundance_cover %>%
left_join(transect_level_corallivory, by = join_by(region, site, transect, depth_m, coral_species)) %>%
mutate(n_colonies_grazed=replace_na(n_colonies_grazed, 0)) %>%
select(region, site, transect, depth_m, coral_species, pct_coral_cover_species, n_colonies_grazed, n_colonies) %>%
filter(coral_species %in% c("Orbicella annularis", "Porites astreoides", "Branching Porites spp.", "Siderastrea siderea")) %>%
left_join(coral_diversity_and_cover, by = join_by(region, site, transect, depth_m)) %>%
mutate(proportional_cover=pct_coral_cover_species/pct_coral_cover)
set.seed(1945) #for reproducible results from random simulations of residuals
#looking at data distributions
transect_corallivory_diversity_df %>%
ggplot(aes(n_colonies_grazed))+
geom_histogram()
#coral cover is severely skewed
transect_corallivory_diversity_df %>%
ggplot(aes(pct_coral_cover_species))+
geom_histogram()
#checking that species percent cover does not have 0's
summary(transect_corallivory_diversity_df$pct_coral_cover_species)
#log-transformation normalized coral cover distribution
transect_corallivory_diversity_df %>%
ggplot(aes(log(pct_coral_cover_species)))+
geom_histogram()
#scar abundance per colony with poisson
n_grazed_diversity_poisson <- glmmTMB::glmmTMB(n_colonies_grazed ~ coral_diversity+ coral_species+coral_species:log(proportional_cover)+ depth_m+ region+
(1 | site/transect), #transects nested within sites within regions
family=poisson,
data = transect_corallivory_diversity_df)
#scar abundance per colony with negative binomial distribution (without zero-inflation)
n_grazed_diversity_nbinom <- glmmTMB::glmmTMB(n_colonies_grazed ~ coral_diversity+ coral_species+ coral_species:log(proportional_cover)+ depth_m+ region+
(1 | site/transect), #transects nested within sites within regions
zi=~0, #non zero-inflated
family=nbinom2(link = "log"),
data = transect_corallivory_diversity_df)
plot(DHARMa::simulateResiduals(n_grazed_diversity_nbinom, n = 1000, use.u = T))
#scar abundance per colony with negative binomial distribution (with zero-inflation)
n_grazed_diversity_nbinom_zi  <- glmmTMB::glmmTMB(n_colonies_grazed ~ coral_diversity+ coral_species+ coral_species:log(proportional_cover)+ depth_m+ region+
(1 | site/transect), #transects nested within sites within regions
zi=~1, # zero-inflated
family=nbinom2(link = "log"),
data = transect_corallivory_diversity_df)
plot(DHARMa::simulateResiduals(n_grazed_diversity_nbinom, n = 1000, use.u = T))
hist(log(transect_corallivory_diversity_df$proportional_cover))
hist(sqrt(transect_corallivory_diversity_df$proportional_cover))
#likelihood ratio tests suggest NS difference between poisson and negative binomial
anova(n_grazed_diversity_poisson, n_grazed_diversity_nbinom)
#likelihood ratio tests suggest NS difference between negative binomial w/ or w/o zero inflation... supports using non zi model
anova(n_grazed_diversity_nbinom_zi, n_grazed_diversity_nbinom)
#simulating residuals and evaluating model fit
n_grazed_diversity_poisson_residuals <- DHARMa::simulateResiduals(n_grazed_diversity_poisson, n = 1000, use.u = T)
plot(n_grazed_diversity_poisson_residuals) #diagnostic plot: KS test significant
n_grazed_diversity_nbinom_residuals <- DHARMa::simulateResiduals(n_grazed_diversity_nbinom, n = 1000, use.u = T)
plot(n_grazed_diversity_nbinom_residuals) #diagnostic plot indicates no sig problems: using nbinom distribution
DHARMa::testZeroInflation(n_grazed_diversity_nbinom_residuals) #secondary check that there's no evidence of zero inflation (in addition to likelihood ratio test)
n_grazed_diversity_nbinom %>% summary()
n_grazed_diversity_nbinom_ci <- confint(n_grazed_diversity_nbinom) %>%
as_tibble() %>%
select(-Estimate)
#tidying model output
transect_level_model_output <- n_grazed_diversity_nbinom %>%
broom.mixed::tidy() %>%
bind_cols(n_grazed_diversity_nbinom_ci) %>%
filter(effect!="ran_pars") %>%
janitor::clean_names() %>%
rename(ci_2.5 = x2_5_percent, ci_97.5=x97_5_percent) %>%
mutate(model="Colonies grazed") %>%
select(model, term, statistic, p_value, estimate, ci_2.5, ci_97.5) %>%
mutate(across(where(is.double), ~round(.x, 3))) %>%
mutate(term=str_replace(term, "coral_species", "")) %>%
mutate(term=str_replace(term, "region", ""))
#back transforming estimated effect sizes (see paper supplements for description)
transect_level_model_output %>%
#exponentiating all independent variables except proportional cover, which was log transformed
mutate(estimate_backtransformed = exp(estimate)) %>%
mutate(across(where(is.double), ~round(.x, 3)))
write_csv(transect_level_model_output, here("figures_tables/supplements_tables/table_s5_transect_level_grazed_coral_abundance_results.csv"))
# total number of coral colonies observed across all transects
nrow(coral_scar_df)
# total number of colonies grazed (for coral taxa with 3+ predated cols across observations)
n_grazed_cols <- coral_scar_df %>%
filter(n_new_scars>0) %>%
distinct(colony_id) %>%
nrow()
n_grazed_cols
# number of grazed colonies by taxa (for coral taxa with 3+ predated cols across observations)
coral_scar_df %>%
filter(n_new_scars>0) %>%
group_by(coral_species) %>%
tally() %>%
arrange(-n)
# total number of colonies grazed for top four predated species
n_top_4_grazed <- coral_scar_df %>%
filter(n_new_scars>0) %>%
filter(coral_species %in% c("Orbicella annularis", "Branching Porites spp.", "Porites astreoides", "Siderastrea siderea")) %>%
distinct(colony_id) %>%
nrow()
#percentage of predation scars observed on these four coral species compared to total number of recent scars
round(n_top_4_grazed/n_grazed_cols*100, 1)
corallivory_by_diversity_fig <- transect_corallivory_diversity_df %>%
mutate(coral_species=factor(coral_species, levels = c("Orbicella annularis", "Branching Porites spp.", "Porites astreoides", "Siderastrea siderea"))) %>%
mutate(region= factor(region, levels = c("Panamá", "Florida",  "St. Croix", "Bonaire"))) %>%
ggplot(aes(x=coral_diversity, y=n_colonies_grazed, color=region))+
geom_point(alpha=0.6)+
scale_y_continuous(breaks = c(0,2,4,6,8,10),
labels= c(0,2,4,6,8,10),
limits=c(0,10.2))+
theme_bw()+
scale_color_manual(values=region_palette)+
scale_fill_manual(values=region_palette)+
facet_wrap(~coral_species)+
labs(x = "Coral diversity (H')",
y = expression(paste("Colonies preyed upon per 30 m" ^2)),
color="Region",
fill="Region")+
theme_bw()+
theme(legend.position = "top")
#reading in coral genus icons
branching_porties_image <- magick::image_read(path=here("figures_tables/reef_images", "Branching_porties.png"), density = NULL, depth = NULL, strip = TRUE)
orbicella_image <- magick::image_read(path=here("figures_tables/reef_images", "Orbicella_annularis.png"), density = NULL, depth = NULL, strip = TRUE)
porites_astreoides_image <- magick::image_read(path=here("figures_tables/reef_images", "Porites_astreoides.png"), density = NULL, depth = NULL, strip = TRUE)
siderastrea_image <- magick::image_read(path=here("figures_tables/reef_images", "Siderastrea_siderea.png"), density = NULL, depth = NULL, strip = TRUE)
#adding coral images to figure (image placement is off in R, but good in ggsave generated png)
corallivory_by_diversity_fig_final <- cowplot::ggdraw() +
draw_plot(corallivory_by_diversity_fig)+ #adds the plot (order matters to overlay the layers correctly)
draw_image(orbicella_image, x = -0.385, y = 0.29, scale=0.075)+
draw_image(branching_porties_image, x = 0.081, y = 0.29, scale=0.08)+
draw_image(porites_astreoides_image, x = -0.385, y = -0.113, scale=0.062)+
draw_image(siderastrea_image, x = 0.08, y = -0.11, scale=0.063)
#saving figure
ggsave(plot=corallivory_by_diversity_fig_final, filename=here("figures_tables/supplements_figs/fig_s5_corallivory_by_coral_diversity.png"), width = 190, height = 140, dpi = 600, units = "mm")
corallivory_by_diversity_fig_final
??glmmTMB