# 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 Water Quality analysis
1. Setup
2. Data preparation and cleaning
# import data
raw_data <- read_csv("water_parameters_raw.csv")
# check data is present with no missing values
head(raw_data)# A tibble: 6 × 8
Module Riparian_condition Time Replicate Temperature_C EC_uS_cm TDS_ppm pH
<chr> <chr> <tim> <dbl> <dbl> <dbl> <dbl> <dbl>
1 CR1 CR 14:30 1 19.2 109. 54.5 7.31
2 CR1 CR 14:30 2 18.9 109. 54.4 7.41
3 CR1 CR 14:30 3 19.2 109. 54.3 7.43
4 CR2 CR 14:45 1 19.4 110. 54.6 7.36
5 CR2 CR 14:45 2 19.4 110. 54.8 7.38
6 CR2 CR 14:45 3 19.3 109 54.4 7.38tail(raw_data)# A tibble: 6 × 8
Module Riparian_condition Time Replicate Temperature_C EC_uS_cm TDS_ppm pH
<chr> <chr> <tim> <dbl> <dbl> <dbl> <dbl> <dbl>
1 IR4 IR 11:45 1 17 111. 55.4 7.51
2 IR4 IR 11:45 2 17 111. 55.5 7.51
3 IR4 IR 11:45 3 17 111. 55.5 7.54
4 IR5 IR 12:00 1 17 111. 55.5 7.46
5 IR5 IR 12:00 2 17 111. 55.4 7.44
6 IR5 IR 12:00 3 17 111. 55.5 7.46str(raw_data)spc_tbl_ [30 × 8] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ Module : chr [1:30] "CR1" "CR1" "CR1" "CR2" ...
$ Riparian_condition: chr [1:30] "CR" "CR" "CR" "CR" ...
$ Time : 'hms' num [1:30] 14:30:00 14:30:00 14:30:00 14:45:00 ...
..- attr(*, "units")= chr "secs"
$ Replicate : num [1:30] 1 2 3 1 2 3 1 2 3 1 ...
$ Temperature_C : num [1:30] 19.2 18.9 19.2 19.4 19.4 19.3 19.5 19.5 19.6 19.6 ...
$ EC_uS_cm : num [1:30] 109 109 109 110 110 ...
$ TDS_ppm : num [1:30] 54.5 54.4 54.3 54.6 54.8 54.4 54.7 54.7 54.5 54.4 ...
$ pH : num [1:30] 7.31 7.41 7.43 7.36 7.38 7.38 7.39 7.39 7.39 7.39 ...
- attr(*, "spec")=
.. cols(
.. Module = col_character(),
.. Riparian_condition = col_character(),
.. Time = col_time(format = ""),
.. Replicate = col_double(),
.. Temperature_C = col_double(),
.. EC_uS_cm = col_double(),
.. TDS_ppm = col_double(),
.. pH = col_double()
.. )
- attr(*, "problems")=<externalptr> nrow(raw_data)[1] 30colSums(is.na(raw_data)) Module Riparian_condition Time Replicate
0 0 0 0
Temperature_C EC_uS_cm TDS_ppm pH
0 0 0 0 # Data cleaning: removing empty rows
raw_data <- raw_data %>%
select(where(~ !all(is.na(.))))
# Check the cleaned structure
glimpse(raw_data)Rows: 30
Columns: 8
$ Module <chr> "CR1", "CR1", "CR1", "CR2", "CR2", "CR2", "CR3", "C…
$ Riparian_condition <chr> "CR", "CR", "CR", "CR", "CR", "CR", "CR", "CR", "CR…
$ Time <time> 14:30:00, 14:30:00, 14:30:00, 14:45:00, 14:45:00, …
$ Replicate <dbl> 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, …
$ Temperature_C <dbl> 19.2, 18.9, 19.2, 19.4, 19.4, 19.3, 19.5, 19.5, 19.…
$ EC_uS_cm <dbl> 109.3, 108.7, 108.7, 109.7, 109.5, 109.0, 109.5, 10…
$ TDS_ppm <dbl> 54.5, 54.4, 54.3, 54.6, 54.8, 54.4, 54.7, 54.7, 54.…
$ pH <dbl> 7.31, 7.41, 7.43, 7.36, 7.38, 7.38, 7.39, 7.39, 7.3…3. Replicate measurements
# Calculating module means and standard deviation of replicates per perimeter
module_summary <- raw_data %>%
group_by(Module, Riparian_condition) %>%
summarise(Temperature_mean = mean(Temperature_C),
Temperature_SD = sd(Temperature_C),
EC_mean = mean(EC_uS_cm),
EC_SD = sd(EC_uS_cm),
TDS_mean = mean(TDS_ppm),
TDS_SD = sd(TDS_ppm),
pH_mean = mean(pH),
pH_SD = sd(pH),
.groups = "drop")
module_summary# A tibble: 10 × 10
Module Riparian_condition Temperature_mean Temperature_SD EC_mean EC_SD
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 CR1 CR 19.1 0.173 109. 0.346
2 CR2 CR 19.4 0.0577 109. 0.361
3 CR3 CR 19.5 0.0577 109. 0.289
4 CR4 CR 19.5 0.100 109. 0.436
5 CR5 CR 19.8 0 109. 0.231
6 IR1 IR 16.9 0 110. 0.208
7 IR2 IR 16.9 0.0577 111. 0.0577
8 IR3 IR 17 0 111 0.608
9 IR4 IR 17 0 111. 0.115
10 IR5 IR 17 0 111. 0.252
# ℹ 4 more variables: TDS_mean <dbl>, TDS_SD <dbl>, pH_mean <dbl>, pH_SD <dbl># export descriptive stats table of module-level summary for appendix
ds_module_table <- module_summary %>%
mutate(Riparian_condition = recode(Riparian_condition,
"CR" = "Tree-lined", "IR" = "Non-tree-lined"),
Temperature_mean = round(Temperature_mean, 2),
Temperature_SD = round(Temperature_SD, 2),
EC_mean = round(EC_mean, 2),
EC_SD = round(EC_SD, 2),
TDS_mean = round(TDS_mean, 2),
TDS_SD = round(TDS_SD, 2),
pH_mean = round(pH_mean, 2),
pH_SD = round(pH_SD, 2)) %>%
rename("Module" = Module,
"Riparian condition" = Riparian_condition,
"Temperature\nmean (°C)" = Temperature_mean,
"Temperature\nSD" = Temperature_SD,
"EC\nmean (µS/cm)" = EC_mean,
"EC\nSD" = EC_SD,
"TDS\nmean (ppm)" = TDS_mean,
"TDS\nSD" = TDS_SD,
"pH\nmean" = pH_mean,
"pH\nSD" = pH_SD)
ds_module_table# A tibble: 10 × 10
Module `Riparian condition` `Temperature\nmean (°C)` `Temperature\nSD`
<chr> <chr> <dbl> <dbl>
1 CR1 Tree-lined 19.1 0.17
2 CR2 Tree-lined 19.4 0.06
3 CR3 Tree-lined 19.5 0.06
4 CR4 Tree-lined 19.5 0.1
5 CR5 Tree-lined 19.8 0
6 IR1 Non-tree-lined 16.9 0
7 IR2 Non-tree-lined 16.9 0.06
8 IR3 Non-tree-lined 17 0
9 IR4 Non-tree-lined 17 0
10 IR5 Non-tree-lined 17 0
# ℹ 6 more variables: `EC\nmean (µS/cm)` <dbl>, `EC\nSD` <dbl>,
# `TDS\nmean (ppm)` <dbl>, `TDS\nSD` <dbl>, `pH\nmean` <dbl>, `pH\nSD` <dbl># create table
ds_module_table <- ds_module_table %>%
flextable() %>%
theme_booktabs() %>%
align(align = "center", part = "header") %>%
align(j = 1:10, align = "center", part = "body") %>%
bold(part = "header") %>%
fontsize(size = 9, part = "header") %>%
fontsize(size = 9, part = "body") %>%
padding(padding = 2, part = "all") %>%
autofit()
ds_module_tableModule | Riparian condition | Temperature | Temperature | EC | EC | TDS | TDS | pH | pH |
|---|---|---|---|---|---|---|---|---|---|
CR1 | Tree-lined | 19.10 | 0.17 | 108.90 | 0.35 | 54.40 | 0.10 | 7.38 | 0.06 |
CR2 | Tree-lined | 19.37 | 0.06 | 109.40 | 0.36 | 54.60 | 0.20 | 7.37 | 0.01 |
CR3 | Tree-lined | 19.53 | 0.06 | 109.33 | 0.29 | 54.63 | 0.12 | 7.39 | 0.00 |
CR4 | Tree-lined | 19.50 | 0.10 | 109.20 | 0.44 | 54.50 | 0.10 | 7.38 | 0.01 |
CR5 | Tree-lined | 19.80 | 0.00 | 108.87 | 0.23 | 54.43 | 0.15 | 7.37 | 0.03 |
IR1 | Non-tree-lined | 16.90 | 0.00 | 110.37 | 0.21 | 55.43 | 0.06 | 7.32 | 0.06 |
IR2 | Non-tree-lined | 16.93 | 0.06 | 110.57 | 0.06 | 55.63 | 0.12 | 7.34 | 0.10 |
IR3 | Non-tree-lined | 17.00 | 0.00 | 111.00 | 0.61 | 55.73 | 0.06 | 7.53 | 0.04 |
IR4 | Non-tree-lined | 17.00 | 0.00 | 110.77 | 0.12 | 55.47 | 0.06 | 7.52 | 0.02 |
IR5 | Non-tree-lined | 17.00 | 0.00 | 111.13 | 0.25 | 55.47 | 0.06 | 7.45 | 0.01 |
save_as_docx(ds_module_table,path = "outputs/tables/A1_module_means.docx")4. Descriptive statistics
# renaming module variables so input and output columns have different names
module_values <- module_summary %>%
rename(Temperature = Temperature_mean,
EC = EC_mean,
TDS = TDS_mean,
pH = pH_mean)
module_values# A tibble: 10 × 10
Module Riparian_condition Temperature Temperature_SD EC EC_SD TDS
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 CR1 CR 19.1 0.173 109. 0.346 54.4
2 CR2 CR 19.4 0.0577 109. 0.361 54.6
3 CR3 CR 19.5 0.0577 109. 0.289 54.6
4 CR4 CR 19.5 0.100 109. 0.436 54.5
5 CR5 CR 19.8 0 109. 0.231 54.4
6 IR1 IR 16.9 0 110. 0.208 55.4
7 IR2 IR 16.9 0.0577 111. 0.0577 55.6
8 IR3 IR 17 0 111 0.608 55.7
9 IR4 IR 17 0 111. 0.115 55.5
10 IR5 IR 17 0 111. 0.252 55.5
# ℹ 3 more variables: TDS_SD <dbl>, pH <dbl>, pH_SD <dbl># calculating the mean and SD of perimeters for both riparian conditions
descriptive_stats <- module_values %>%
group_by(Riparian_condition) %>%
summarise(Temperature_mean = mean(Temperature),
Temperature_SD = sd(Temperature),
EC_mean = mean(EC),
EC_SD = sd(EC),
TDS_mean = mean(TDS),
TDS_SD = sd(TDS),
pH_mean = mean(pH),
pH_SD = sd(pH),
.groups = "drop")
descriptive_stats# A tibble: 2 × 9
Riparian_condition Temperature_mean Temperature_SD EC_mean EC_SD TDS_mean
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 CR 19.5 0.255 109. 0.245 54.5
2 IR 17.0 0.0471 111. 0.312 55.5
# ℹ 3 more variables: TDS_SD <dbl>, pH_mean <dbl>, pH_SD <dbl># export descriptive stats table of riparian conditions
ds_reach_table <- descriptive_stats %>%
mutate(Riparian_condition = recode(Riparian_condition,
"CR" = "Tree-lined", "IR" = "Non-tree-lined"))
# renaming and resizing columns
ds_reach_table <- ds_reach_table %>%
rename("Riparian Condition" = Riparian_condition,
"Temperature Mean" = Temperature_mean,
"Temperature SD" = Temperature_SD,
"EC Mean" = EC_mean,
"EC SD" = EC_SD,
"TDS Mean" = TDS_mean,
"TDS SD" = TDS_SD,
"pH Mean" = pH_mean,
"pH SD" = pH_SD)
# rounding values
ds_reach_table <- ds_reach_table %>%
mutate(across(where(is.numeric), ~ round(.x, 2)))
ds_reach_table# A tibble: 2 × 9
`Riparian Condition` `Temperature Mean` `Temperature SD` `EC Mean` `EC SD`
<chr> <dbl> <dbl> <dbl> <dbl>
1 Tree-lined 19.5 0.26 109. 0.25
2 Non-tree-lined 17.0 0.05 111. 0.31
# ℹ 4 more variables: `TDS Mean` <dbl>, `TDS SD` <dbl>, `pH Mean` <dbl>,
# `pH SD` <dbl># creating the table
ds_reach_table <- flextable(ds_reach_table) %>%
add_header_row(values = c("Riparian condition",
"Temperature (°C)",
"EC (µS/cm)",
"TDS (ppm)",
"pH"),
colwidths = c(1, 2, 2, 2, 2)) %>%
set_header_labels("Riparian condition" = "Riparian condition",
"Temperature Mean" = "Mean",
"Temperature SD" = "SD",
"EC Mean" = "Mean",
"EC SD" = "SD",
"TDS Mean" = "Mean",
"TDS SD" = "SD",
"pH Mean" = "Mean",
"pH SD" = "SD")
ds_reach_tableRiparian condition | Temperature (°C) | EC (µS/cm) | TDS (ppm) | pH | ||||
|---|---|---|---|---|---|---|---|---|
Riparian Condition | Mean | SD | Mean | SD | Mean | SD | Mean | SD |
Tree-lined | 19.46 | 0.26 | 109.14 | 0.25 | 54.51 | 0.10 | 7.38 | 0.01 |
Non-tree-lined | 16.97 | 0.05 | 110.77 | 0.31 | 55.55 | 0.13 | 7.43 | 0.10 |
save_as_docx(ds_reach_table,path = "outputs/tables/A2_descriptive.docx")
# visualising descritive statistics of each perimeter
# Temperature: jittered boxplot
ds_temperature <- ggplot(module_summary,
aes(x = Riparian_condition,
y = Temperature_mean)) +
geom_boxplot(outlier.shape = NA,
fill = "gray95") +
geom_jitter(aes(color = "Individual modules"),
width = 0.1, size = 1.5) +
labs(x = "Riparian condition",
y = "Water temperature (°C)",
color = NULL) +
scale_x_discrete(labels = c("CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
theme_light()
ds_temperature
# EC: jittered boxplot
ds_EC <- ggplot(module_summary,
aes(x = Riparian_condition,
y = EC_mean)) +
geom_boxplot(outlier.shape = NA,
fill = "gray95") +
geom_jitter(aes(color = "Individual modules"),
width = 0.1, size = 1.5) +
labs(x = "Riparian condition",
y = "Electrical conductivity (µS/cm)",
color = NULL) +
scale_x_discrete(labels = c("CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
theme_light()
ds_EC
# TDS: jittered boxplot
ds_TDS <- ggplot(module_summary,
aes(x = Riparian_condition,
y = TDS_mean)) +
geom_boxplot(outlier.shape = NA,
fill = "gray95") +
geom_jitter(aes(color = "Individual modules"),
width = 0.1, size = 1.5) +
labs(x = "Riparian condition",
y = "Total dissolved solids (ppm)",
color = NULL) +
scale_x_discrete(labels = c("CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
theme_light()
ds_TDS
# pH: jittered boxplot
ds_pH <- ggplot(module_summary,
aes(x = Riparian_condition,
y = pH_mean)) +
geom_boxplot(outlier.shape = NA,
fill = "gray95") +
geom_jitter(aes(color = "Individual modules"),
width = 0.1, size = 1.5) +
labs(x = "Riparian condition",
y = "pH",
color = NULL) +
scale_x_discrete(labels = c("CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
theme_light()
ds_pH
# combining the jittered boxplots to export for the main text
ds_modules_figure <- ds_temperature + ds_EC + ds_TDS + ds_pH +
plot_layout(ncol = 2, guides = "collect") &
theme(axis.title = element_text(size = 9),
axis.text.x = element_text(size = 8),
axis.text.y = element_text(size = 8),
legend.position = "bottom",
legend.text = element_text(size = 8),
legend.margin = margin(t = -5))
ds_modules_figure
# save figure
ggsave("outputs/figures/1_descriptive.png",
plot = ds_modules_figure,
width = 18,
height = 15,
units = "cm",
dpi = 300)5. Exploratory diagnostics
# Inspecting individual module values
module_values %>% select(Riparian_condition, Module, Temperature, EC, TDS, pH)# A tibble: 10 × 6
Riparian_condition Module Temperature EC TDS pH
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 CR CR1 19.1 109. 54.4 7.38
2 CR CR2 19.4 109. 54.6 7.37
3 CR CR3 19.5 109. 54.6 7.39
4 CR CR4 19.5 109. 54.5 7.38
5 CR CR5 19.8 109. 54.4 7.37
6 IR IR1 16.9 110. 55.4 7.32
7 IR IR2 16.9 111. 55.6 7.34
8 IR IR3 17 111 55.7 7.53
9 IR IR4 17 111. 55.5 7.52
10 IR IR5 17 111. 55.5 7.45# Looking at visual distribution of data using Normality diagnostic QQ plots
# Temperature QQ plot
qq_temperature <- ggplot(module_values, aes(sample = Temperature)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~ Riparian_condition,
labeller = as_labeller(c(
CR = "Tree-lined",
IR = "Non-tree-lined"))) +
labs(title = "Temperature",
x = "Theoretical quantiles",
y = "Sample quantiles") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
# EC QQ plot
qq_EC <- ggplot(module_values, aes(sample = EC)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~ Riparian_condition,
labeller = as_labeller(c(
CR = "Tree-lined",
IR = "Non-tree-lined"))) +
labs(title = "EC",
x = "Theoretical quantiles",
y = "Sample quantiles") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
# TDS QQ plot
qq_TDS <- ggplot(module_values, aes(sample = TDS)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~ Riparian_condition,
labeller = as_labeller(c(
CR = "Tree-lined",
IR = "Non-tree-lined"))) +
labs(title = "TDS",
x = "Theoretical quantiles",
y = "Sample quantiles") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
# pH QQ plot
qq_pH <- ggplot(module_values, aes(sample = pH)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~ Riparian_condition,
labeller = as_labeller(c(
CR = "Tree-lined",
IR = "Non-tree-lined"))) +
labs(title = "pH",
x = "Theoretical quantiles",
y = "Sample quantiles") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))
# Combine into one singular figure for export
qq_figure <- (qq_temperature | qq_EC) /
(qq_TDS | qq_pH)
qq_figure
# Export QQ plots for appendix
ggsave("outputs/figures/A1_qq.png",
plot = qq_figure,
width = 15,
height = 10,
units = "cm",
dpi = 300)6. Spatial analysis
# Creating a numeric spatial position from the module ID
spatial_data <- module_summary %>%
mutate(Module_position = as.numeric(str_extract(Module, "\\d+")))
spatial_data# A tibble: 10 × 11
Module Riparian_condition Temperature_mean Temperature_SD EC_mean EC_SD
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 CR1 CR 19.1 0.173 109. 0.346
2 CR2 CR 19.4 0.0577 109. 0.361
3 CR3 CR 19.5 0.0577 109. 0.289
4 CR4 CR 19.5 0.100 109. 0.436
5 CR5 CR 19.8 0 109. 0.231
6 IR1 IR 16.9 0 110. 0.208
7 IR2 IR 16.9 0.0577 111. 0.0577
8 IR3 IR 17 0 111 0.608
9 IR4 IR 17 0 111. 0.115
10 IR5 IR 17 0 111. 0.252
# ℹ 5 more variables: TDS_mean <dbl>, TDS_SD <dbl>, pH_mean <dbl>, pH_SD <dbl>,
# Module_position <dbl># Creating spatial module plots
# Temperature
spatial_temperature <- ggplot(spatial_data,
aes(x = Module_position,
y = Temperature_mean,
colour = Riparian_condition,
group = Riparian_condition)) +
geom_line() +
geom_point(size = 2) +
scale_x_continuous(breaks = 1:5,
labels = 1:5) +
scale_colour_discrete(labels = c(
"CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
labs(x = "Module",
y = "Water temperature (°C)",
colour = "Riparian condition") +
theme_light() +
theme(legend.position = "bottom",
legend.title = element_blank())
spatial_temperature
# EC
spatial_ec <- ggplot(spatial_data,
aes(x = Module_position,
y = EC_mean,
colour = Riparian_condition,
group = Riparian_condition)) +
geom_line() +
geom_point(size = 2) +
scale_x_continuous(breaks = 1:5,
labels = 1:5) +
scale_colour_discrete(labels = c(
"CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
labs(x = "Module",
y = "EC (µS/cm)",
colour = "Riparian condition") +
theme_light()
spatial_ec
# TDS
spatial_tds <- ggplot(spatial_data,
aes(x = Module_position,
y = TDS_mean,
colour = Riparian_condition,
group = Riparian_condition)) +
geom_line() +
geom_point(size = 2) +
scale_x_continuous(breaks = 1:5,
labels = 1:5) +
scale_colour_discrete(labels = c(
"CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
labs(x = "Module",
y = "TDS (ppm)",
colour = "Riparian condition") +
theme_light()
spatial_tds
# pH
spatial_ph <- ggplot(spatial_data,
aes(x = Module_position,
y = pH_mean,
colour = Riparian_condition,
group = Riparian_condition)) +
geom_line() +
geom_point(size = 2) +
scale_x_continuous(breaks = 1:5,
labels = 1:5) +
scale_colour_discrete(labels = c(
"CR" = "Tree-lined",
"IR" = "Non-tree-lined")) +
labs(x = "Module",
y = "pH",
colour = "Riparian condition") +
theme_light()
spatial_ph
# Combine the plots and export them
spatial_figure <- (spatial_temperature | spatial_ec) /
(spatial_tds | spatial_ph) +
plot_layout(guides = "collect") &
theme(legend.position = "bottom",
legend.title = element_blank())
spatial_figure
ggsave("outputs/figures/2_spatial_water_quality.png",
plot = spatial_figure,
width = 18,
height = 15,
units = "cm",
dpi = 300)