# load packages
library(tidyverse)
library(janitor)
library(flextable)
library(patchwork)
library(effectsize)
set.seed(42)
# create outputs folder
dir.create("outputs", showWarnings = TRUE)
dir.create("outputs/figures", showWarnings = TRUE)
dir.create("outputs/tables", showWarnings = TRUE)Welton Beck Macroinvertebrate Analysis
1. Setup
2. Data preparation and verification
# import data
raw_data <- read_csv("raw_riverfly.csv")
# check data is present with no missing values
head(raw_data)# A tibble: 6 × 12
site reach feature `phylum/group` class order family common_name
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 CR CR1 Tree-lined Crustacea Malacostraca Amphipo… Gamma… freshwater…
2 CR CR1 Tree-lined Crustacea Malacostraca Isopoda Asell… water slat…
3 CR CR1 Tree-lined Mollusca Gastropoda Hygroph… Plano… ramshorn s…
4 CR CR1 Tree-lined Mollusca Gastropoda Hygroph… Lymna… pond snail
5 CR CR1 Tree-lined Mollusca Gastropoda unident… unide… unidentifi…
6 CR CR1 Tree-lined Annelida Clitellata Hirudin… Gloss… leeches
# ℹ 4 more variables: BMWP_category <chr>, BMWP_score <dbl>, present <dbl>,
# BMWP_contribution <dbl>tail(raw_data)# A tibble: 6 × 12
site reach feature `phylum/group` class order family common_name
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 IR IR5 Non-tree-lined Arthropoda Insecta Hemiptera Corix… lesser wat…
2 IR IR5 Non-tree-lined Arthropoda Insecta Hemiptera Nepid… water scor…
3 IR IR5 Non-tree-lined Arthropoda Insecta Hemiptera Gerri… pond skater
4 IR IR5 Non-tree-lined Arthropoda Insecta Ephemero… Caeni… mayfly nym…
5 IR IR5 Non-tree-lined Arthropoda Insecta Ephemero… Baeti… mayfly nym…
6 IR IR5 Non-tree-lined Arthropoda Insecta Neuropte… Siali… alder flies
# ℹ 4 more variables: BMWP_category <chr>, BMWP_score <dbl>, present <dbl>,
# BMWP_contribution <dbl>str(raw_data)spc_tbl_ [240 × 12] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ site : chr [1:240] "CR" "CR" "CR" "CR" ...
$ reach : chr [1:240] "CR1" "CR1" "CR1" "CR1" ...
$ feature : chr [1:240] "Tree-lined" "Tree-lined" "Tree-lined" "Tree-lined" ...
$ phylum/group : chr [1:240] "Crustacea" "Crustacea" "Mollusca" "Mollusca" ...
$ class : chr [1:240] "Malacostraca" "Malacostraca" "Gastropoda" "Gastropoda" ...
$ order : chr [1:240] "Amphipoda" "Isopoda" "Hygrophila" "Hygrophila" ...
$ family : chr [1:240] "Gammaridae" "Asellidae" "Planorbidae" "Lymnaeidae" ...
$ common_name : chr [1:240] "freshwater shrimp" "water slater" "ramshorn snail" "pond snail" ...
$ BMWP_category : chr [1:240] "Gammaridae" "Asellidae" "Other Mollusca families" "Other Mollusca families" ...
$ BMWP_score : num [1:240] 6 3 3 3 3 3 3 4 5 5 ...
$ present : num [1:240] 1 1 1 1 0 1 0 1 0 1 ...
$ BMWP_contribution: num [1:240] 6 3 3 0 0 3 0 4 0 5 ...
- attr(*, "spec")=
.. cols(
.. site = col_character(),
.. reach = col_character(),
.. feature = col_character(),
.. `phylum/group` = col_character(),
.. class = col_character(),
.. order = col_character(),
.. family = col_character(),
.. common_name = col_character(),
.. BMWP_category = col_character(),
.. BMWP_score = col_double(),
.. present = col_double(),
.. BMWP_contribution = col_double()
.. )
- attr(*, "problems")=<externalptr> nrow(raw_data)[1] 240colSums(is.na(raw_data)) site reach feature phylum/group
0 0 0 0
class order family common_name
0 0 0 0
BMWP_category BMWP_score present BMWP_contribution
0 0 0 0 # verify experimental sampling structure
table(raw_data$site)
CR IR
120 120 table(raw_data$reach)
CR1 CR2 CR3 CR4 CR5 IR1 IR2 IR3 IR4 IR5
24 24 24 24 24 24 24 24 24 24 table(raw_data$feature)
Non-tree-lined Tree-lined
120 120 # check BMWP variables
summary(raw_data$BMWP_score) Min. 1st Qu. Median Mean 3rd Qu. Max.
2.000 3.000 5.000 4.792 5.000 10.000 summary(raw_data$BMWP_contribution) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 0.000 1.688 3.000 10.000 table(raw_data$present)
0 1
140 100 3. Module-level summary
# creating 10 module-level summaries
module_data <- raw_data %>%
group_by(site, reach, feature) %>%
summarise( BMWP = sum(BMWP_contribution),
NTAXA = sum(present),
.groups = "drop") %>%
mutate( ASPT = BMWP / NTAXA )
module_data# A tibble: 10 × 6
site reach feature BMWP NTAXA ASPT
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 CR CR1 Tree-lined 36 10 3.6
2 CR CR2 Tree-lined 17 5 3.4
3 CR CR3 Tree-lined 54 12 4.5
4 CR CR4 Tree-lined 54 12 4.5
5 CR CR5 Tree-lined 38 8 4.75
6 IR IR1 Non-tree-lined 28 7 4
7 IR IR2 Non-tree-lined 56 16 3.5
8 IR IR3 Non-tree-lined 52 12 4.33
9 IR IR4 Non-tree-lined 35 8 4.38
10 IR IR5 Non-tree-lined 35 10 3.5 # exporting module-level summary table for appendix
module_table <- module_data %>%
mutate(ASPT = round(ASPT, 2)) %>%
select(site, reach, feature, BMWP, NTAXA, ASPT) %>%
rename(Site = site,
"Module" = reach,
"Riparian condition" = feature,
BMWP = BMWP,
NTAXA = NTAXA,
ASPT = ASPT)
module_table# A tibble: 10 × 6
Site Module `Riparian condition` BMWP NTAXA ASPT
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 CR CR1 Tree-lined 36 10 3.6
2 CR CR2 Tree-lined 17 5 3.4
3 CR CR3 Tree-lined 54 12 4.5
4 CR CR4 Tree-lined 54 12 4.5
5 CR CR5 Tree-lined 38 8 4.75
6 IR IR1 Non-tree-lined 28 7 4
7 IR IR2 Non-tree-lined 56 16 3.5
8 IR IR3 Non-tree-lined 52 12 4.33
9 IR IR4 Non-tree-lined 35 8 4.38
10 IR IR5 Non-tree-lined 35 10 3.5 # format as flextable
module_table <- module_table %>%
flextable() %>%
theme_booktabs() %>%
autofit() %>%
align(align = "center", part = "all")
module_tableSite | Module | Riparian condition | BMWP | NTAXA | ASPT |
|---|---|---|---|---|---|
CR | CR1 | Tree-lined | 36 | 10 | 3.60 |
CR | CR2 | Tree-lined | 17 | 5 | 3.40 |
CR | CR3 | Tree-lined | 54 | 12 | 4.50 |
CR | CR4 | Tree-lined | 54 | 12 | 4.50 |
CR | CR5 | Tree-lined | 38 | 8 | 4.75 |
IR | IR1 | Non-tree-lined | 28 | 7 | 4.00 |
IR | IR2 | Non-tree-lined | 56 | 16 | 3.50 |
IR | IR3 | Non-tree-lined | 52 | 12 | 4.33 |
IR | IR4 | Non-tree-lined | 35 | 8 | 4.38 |
IR | IR5 | Non-tree-lined | 35 | 10 | 3.50 |
save_as_docx(module_table,path = "outputs/tables/A1_module_summary.docx")4. Descriptive Statistics
# summarising the 10 module-level BMWP, NTAXA, and ASPT data into the two
# riparian conditions
descriptive_stats <- module_data %>%
group_by(feature) %>%
summarise(n = n(),
BMWP_mean = mean(BMWP),
BMWP_sd = sd(BMWP),
BMWP_median = median(BMWP),
BMWP_min = min(BMWP),
BMWP_max = max(BMWP),
NTAXA_mean = mean(NTAXA),
NTAXA_sd = sd(NTAXA),
NTAXA_median = median(NTAXA),
NTAXA_min = min(NTAXA),
NTAXA_max = max(NTAXA),
ASPT_mean = mean(ASPT),
ASPT_sd = sd(ASPT),
ASPT_median = median(ASPT),
ASPT_min = min(ASPT),
ASPT_max = max(ASPT),
.groups = "drop")
view(descriptive_stats)
# creating a descriptive statistics table with mean and standard deviation of
# each riparian condition
ds_table <- descriptive_stats %>%
mutate(BMWP = paste0(round(BMWP_mean, 2), " ± ",
round(BMWP_sd,2)),
NTAXA = paste0(round(NTAXA_mean, 2), " ± ",
round(NTAXA_sd, 2)),
ASPT = paste0(round(ASPT_mean, 2), " ± ",
round(ASPT_sd,2))) %>%
select(feature, n, BMWP, NTAXA, ASPT)
ds_table# A tibble: 2 × 5
feature n BMWP NTAXA ASPT
<chr> <int> <chr> <chr> <chr>
1 Non-tree-lined 5 41.2 ± 12.11 10.6 ± 3.58 3.94 ± 0.43
2 Tree-lined 5 39.8 ± 15.34 9.4 ± 2.97 4.15 ± 0.61ds_table <- ds_table %>%
rename("Riparian condition" = feature,
"BMWP\n(mean ± SD)" = BMWP,
"NTAXA\n(mean ± SD)" = NTAXA,
"ASPT\n(mean ± SD)" = ASPT) %>%
flextable() %>%
theme_booktabs() %>%
autofit() %>%
align(align = "center", part = "all")
ds_tableRiparian condition | n | BMWP | NTAXA | ASPT |
|---|---|---|---|---|
Non-tree-lined | 5 | 41.2 ± 12.11 | 10.6 ± 3.58 | 3.94 ± 0.43 |
Tree-lined | 5 | 39.8 ± 15.34 | 9.4 ± 2.97 | 4.15 ± 0.61 |
save_as_docx(ds_table,path = "outputs/tables/1_Descriptive.docx")
# creating a jitter box plot presenting mean and standard deviation alongside
# individual module observations
# BMWP jittered boxplot
bmwp_ds <- ggplot( module_data,
aes( x = feature,
y = BMWP)) +
geom_boxplot(outlier.shape = NA,
fill = "gray95") +
geom_jitter( aes(color = "Individual 10m modules"),
width = 0.1, size = 1.5 ) +
labs( x = NULL,
y = "BMWP score",
color = NULL) +
theme_light()
bmwp_ds
# NTAXA jittered boxplot
ntaxa_ds <- ggplot( module_data,
aes( x = feature,
y = NTAXA)) +
geom_boxplot(outlier.shape = NA,
fill = "gray95") +
geom_jitter( aes(color = "Individual 10m modules"),
width = 0.1, size = 1.5 ) +
labs( x = "Riparian condition",
y = "Number of BMWP-scoring\ntaxa (NTAXA)",
color = NULL) +
theme_light()
ntaxa_ds
# ASPT jittered boxplot
aspt_ds <- ggplot( module_data,
aes( x = feature,
y = ASPT)) +
geom_boxplot(outlier.shape = NA,
fill = "gray95") +
geom_jitter( aes(color = "Individual 10m modules"),
width = 0.1, size = 1.5 ) +
labs( x = NULL,
y = "Average BMWP-score per\ntaxon (ASPT)",
color = NULL) +
theme_light()
aspt_ds
# merge 3 jitter boxplots into one line
plots_ds <- bmwp_ds + ntaxa_ds + aspt_ds +
plot_layout(ncol = 3, guides = "collect") &
theme(axis.title = element_text(size = 8),
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
legend.position = "bottom",
legend.text = element_text(size = 8),
legend.margin = margin(t = -10))
plots_ds
# save to folder
ggsave("outputs/figures/1_conditions.png",
plot = plots_ds,
width = 18,
height = 10,
units = "cm",
dpi = 300)5. Exploratory diagnostics
# Inspecting individual module values
module_data %>% select(feature, reach, BMWP, NTAXA, ASPT)# A tibble: 10 × 5
feature reach BMWP NTAXA ASPT
<chr> <chr> <dbl> <dbl> <dbl>
1 Tree-lined CR1 36 10 3.6
2 Tree-lined CR2 17 5 3.4
3 Tree-lined CR3 54 12 4.5
4 Tree-lined CR4 54 12 4.5
5 Tree-lined CR5 38 8 4.75
6 Non-tree-lined IR1 28 7 4
7 Non-tree-lined IR2 56 16 3.5
8 Non-tree-lined IR3 52 12 4.33
9 Non-tree-lined IR4 35 8 4.38
10 Non-tree-lined IR5 35 10 3.5 # Restructuring data for diagnostics
diagnostic_data <- module_data %>%
select(feature, BMWP, NTAXA, ASPT) %>%
pivot_longer(cols = c(BMWP, NTAXA, ASPT),
names_to = "metric",
values_to = "value")
# Looking at visual distribution of data using Normality diagnostic QQ plots
diagnostic_qqplot <- ggplot(diagnostic_data, aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_grid(feature ~ metric, scales = "free") +
labs(x = "Theoretical quantiles",
y = "Sample quantiles") +
theme_light()
diagnostic_qqplot
# Exporting as a supporting diagnostics in the Appendix
ggsave("outputs/figures/A1_qq.png",
plot = diagnostic_qqplot,
width = 8,
height = 6,
dpi = 300)6. Spatial analysis
# Creating a spatial module plot
# Setting the module order from downstream to upstream
module_data <- module_data %>%
mutate(reach = factor(reach,
levels = c("CR1", "CR2", "CR3", "CR4", "CR5",
"IR1", "IR2", "IR3", "IR4", "IR5")))
# Reordering data
spatial_data <- module_data %>%
mutate(module_position = c(1, 2, 3, 4, 5,
1, 2, 3, 4, 5)) %>%
select(feature, reach, module_position, BMWP, NTAXA, ASPT) %>%
pivot_longer(cols = c(BMWP, NTAXA, ASPT),
names_to = "metric",
values_to = "value")
# Plotting spatial variation in the three macroinvertebrate metrics
spatial_plot <- ggplot(spatial_data,
aes(x = reach,
y = value,
group = feature)) +
geom_line() +
geom_point(size = 2) +
facet_wrap(~ metric,
scales = "free_y",
ncol = 1) +
labs(x = "Sampling module",
y = NULL) +
theme_light() +
theme(axis.text.x = element_text(size = 11, margin = margin(t = 10)),
axis.text.y = element_text(size = 11),
strip.text = element_text(size = 11),
axis.title = element_text(size = 11))
spatial_plot
# Exporting spatial variation figure for appendix
ggsave("outputs/figures/A2_spatial_variation.png",
plot = spatial_plot,
width = 7,
height = 10,
dpi = 300)