
##################################################
##################################################
##                                              ##
##       Descriptive statistics and plots       ##
##                                              ##
##################################################
##################################################

# R version used:
# > R.version
# "R version 4.4.2 (2024-10-31)"

# Clear previous data
rm(list = ls())

# Load packages
require(dplyr)
require(ggplot2)
require(openxlsx)
require(stringr)

# Loading data file
fishdat <- openxlsx::read.xlsx("Supplementary_file_S1.xlsx")


##################################################
# Plot of disturbance vs. prevalence
##################################################

# Select relevant variables
fishdat2 <- fishdat %>% dplyr::select(Disturbance_coded,
                                      Percent_injury_or_deformity) %>%
  na.omit()
fishdat2$Disturbance_coded <- as.factor(fishdat2$Disturbance_coded)

# Grouping data for adding sample size to plot
fishdat2_samp_size <- fishdat2 %>% 
  dplyr::group_by(Disturbance_coded) %>% 
  dplyr::tally()

# Produce plot
ggplot2::ggplot(fishdat2, aes(y = Percent_injury_or_deformity,
                              x = Disturbance_coded)) +
  geom_violin(width = 1) +
  geom_boxplot(width = 0.15, outlier.shape = NA) +
  geom_jitter(shape=21, width = 0.02, alpha = 1, fill = "lightgray", color = "darkgray") +
  xlab("Disturbance") +
  ylab("Percent of sampled indidividuals\nwith injury and/or deformity") +
  scale_x_discrete(labels = c("No", "Yes"),
                   expand = c(1, 0)) +
  scale_y_continuous(label = scales::percent_format(scale = 1)) +
  geom_text(data = fishdat2_samp_size,
            aes(Disturbance_coded, Inf, label = n),
            vjust = 1.2, hjust = 0.48) +
  theme_bw(base_size = 14)

# Save plot - using set.seed() for replication of random placement of 'geom_jitter' function
set.seed(42);ggplot2::ggsave("Injury_deformity.eps",
                             width = 7,
                             height = 5,
                             device = "eps")


##################################################
# Calculating means and SDs for disturbance
##################################################

fishdat2_only_disturbed <- fishdat2 %>%
  dplyr::filter(Disturbance_coded == 1)

fishdat2_only_undisturbed <- fishdat2 %>%
  dplyr::filter(Disturbance_coded == 0)

# Anthropogenic disturbance present
length(fishdat2_only_disturbed$Percent_injury_or_deformity) # sample size
round(mean(fishdat2_only_disturbed$Percent_injury_or_deformity), 1) # mean
round(sd(fishdat2_only_disturbed$Percent_injury_or_deformity), 1) # standard deviation

# Anthropogenic disturbance not present
length(fishdat2_only_undisturbed$Percent_injury_or_deformity) # sample size
round(mean(fishdat2_only_undisturbed$Percent_injury_or_deformity), 1) # mean
round(sd(fishdat2_only_undisturbed$Percent_injury_or_deformity), 1) # standard deviation


##################################################
# Calculating means and SDs for all data combined
##################################################

# Select relevant variable
fishdat_all_prevalence_estimates <- fishdat %>% dplyr::select(Percent_injury_or_deformity) %>%
  na.omit()

# All data
length(fishdat_all_prevalence_estimates$Percent_injury_or_deformity) # sample size
round(mean(fishdat_all_prevalence_estimates$Percent_injury_or_deformity), 1) # mean
round(sd(fishdat_all_prevalence_estimates$Percent_injury_or_deformity), 1) # standard deviation


##################################################
# Calculating means and SDs for injury/deformity
##################################################

fishdat3_only_injuries <- fishdat %>%
  dplyr::filter(Injury_or_deformity_coded == "1")

fishdat3_only_deformities <- fishdat %>%
  dplyr::filter(Injury_or_deformity_coded == "0")

# Injuries
length(fishdat3_only_injuries$Percent_injury_or_deformity) # sample size
round(mean(fishdat3_only_injuries$Percent_injury_or_deformity), 1) # mean
round(sd(fishdat3_only_injuries$Percent_injury_or_deformity), 1) # standard deviation

# Deformities
length(fishdat3_only_deformities$Percent_injury_or_deformity) # sample size
round(mean(fishdat3_only_deformities$Percent_injury_or_deformity), 1) # mean
round(sd(fishdat3_only_deformities$Percent_injury_or_deformity), 1) # standard deviation


##################################################
# Sample sizes for percent estimates
##################################################

# Extract sample sizes
sample_size_dat <- fishdat %>% dplyr::select(Sample_sizes_for_percent_estimates) %>%
  na.omit()
sample_size_dat2 <- stringr::str_split(sample_size_dat$Sample_sizes_for_percent_estimates, ",")
sample_size_dat3 <- as.numeric(unlist(sample_size_dat2))

round(mean(sample_size_dat3), 0) # mean
round(sd(sample_size_dat3), 0) # standard deviation
max(sample_size_dat3) # range upper
min(sample_size_dat3) # range lower
length(sample_size_dat3) # sample size of sample sizes
