# 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 MoRPH analysis
1. Setup
2. Data preparation and cleaning
# import data
raw_data <- read_csv("raw_data.csv")
# check data is present with no missing values
names(raw_data) [1] "...1" "CR1" "CR2" "CR3" "CR4" "CR5" "IR1" "IR2" "IR3" "IR4"
[11] "IR5" view(raw_data)
# restructure MoRPH data into a tidy format
morph_data <- raw_data %>%
rename(index = `...1`) %>%
pivot_longer(cols = CR1:IR5,
names_to = "module",
values_to = "value")
morph_data# A tibble: 140 × 3
index module value
<chr> <chr> <chr>
1 INDEX 1 - Number of flow types CR1 1
2 INDEX 1 - Number of flow types CR2 1
3 INDEX 1 - Number of flow types CR3 1
4 INDEX 1 - Number of flow types CR4 1
5 INDEX 1 - Number of flow types CR5 1
6 INDEX 1 - Number of flow types IR1 2
7 INDEX 1 - Number of flow types IR2 2
8 INDEX 1 - Number of flow types IR3 2
9 INDEX 1 - Number of flow types IR4 2
10 INDEX 1 - Number of flow types IR5 1
# ℹ 130 more rows# include the riparian condition and add numbers 1-5 for individual modules
morph_data <- morph_data %>%
mutate(condition = if_else(
str_starts(module, "CR"),
"Tree-lined",
"Non-tree-lined"),
module_position = as.numeric(str_extract(module, "\\d+")))
morph_data# A tibble: 140 × 5
index module value condition module_position
<chr> <chr> <chr> <chr> <dbl>
1 INDEX 1 - Number of flow types CR1 1 Tree-lined 1
2 INDEX 1 - Number of flow types CR2 1 Tree-lined 2
3 INDEX 1 - Number of flow types CR3 1 Tree-lined 3
4 INDEX 1 - Number of flow types CR4 1 Tree-lined 4
5 INDEX 1 - Number of flow types CR5 1 Tree-lined 5
6 INDEX 1 - Number of flow types IR1 2 Non-tree-lined 1
7 INDEX 1 - Number of flow types IR2 2 Non-tree-lined 2
8 INDEX 1 - Number of flow types IR3 2 Non-tree-lined 3
9 INDEX 1 - Number of flow types IR4 2 Non-tree-lined 4
10 INDEX 1 - Number of flow types IR5 1 Non-tree-lined 5
# ℹ 130 more rows# Isolating the indices to be analysed further
morph_index <- morph_data %>%
filter(index %in% c(
"INDEX 1 - Number of flow types",
"INDEX 3 - Number of bed material types",
"INDEX 7 - Extent of superficial bed siltation",
"INDEX 8 - Channel physical habitat complexity",
"INDEX 9 - Number of aquatic vegetation morphotypes",
"INDEX 10 - Riparian physical habitat complexity",
"INDEX 11 - Riparian vegetation structural complexity",
"INDEX 12 - Human pressure imposed by bank top land cover"))
view(morph_index)
unique(morph_index$index)[1] "INDEX 1 - Number of flow types"
[2] "INDEX 3 - Number of bed material types"
[3] "INDEX 7 - Extent of superficial bed siltation"
[4] "INDEX 8 - Channel physical habitat complexity"
[5] "INDEX 9 - Number of aquatic vegetation morphotypes"
[6] "INDEX 10 - Riparian physical habitat complexity"
[7] "INDEX 11 - Riparian vegetation structural complexity"
[8] "INDEX 12 - Human pressure imposed by bank top land cover"# checking values recorded for each index
morph_index %>%
group_by(index) %>%
summarise(values = paste(unique(value), collapse = ", "))# A tibble: 8 × 2
index values
<chr> <chr>
1 INDEX 1 - Number of flow types 1, 2
2 INDEX 10 - Riparian physical habitat complexity 3, 4.071428571, 3.75…
3 INDEX 11 - Riparian vegetation structural complexity 4, 5, 5.5, 7, 3.5, 4…
4 INDEX 12 - Human pressure imposed by bank top land cover 3.5, 5
5 INDEX 3 - Number of bed material types 2, 1, 3, 4
6 INDEX 7 - Extent of superficial bed siltation 7.65, 6.8, 6.7, 3.55…
7 INDEX 8 - Channel physical habitat complexity 2.083333333, 1.25, 2…
8 INDEX 9 - Number of aquatic vegetation morphotypes 1, 2, 3, 4 # shortening index names
morph_plot_data <- morph_index %>%
mutate(metric = case_when(
index == "INDEX 1 - Number of flow types" ~ "Number of flow types",
index == "INDEX 3 - Number of bed material types" ~ "Number of bed material types",
index == "INDEX 7 - Extent of superficial bed siltation" ~ "Bed siltation",
index == "INDEX 8 - Channel physical habitat complexity" ~ "Channel habitat complexity",
index == "INDEX 9 - Number of aquatic vegetation morphotypes" ~ "Aquatic vegetation morphotypes",
index == "INDEX 10 - Riparian physical habitat complexity" ~ "Riparian habitat complexity",
index == "INDEX 11 - Riparian vegetation structural complexity" ~ "Riparian vegetation complexity",
index == "INDEX 12 - Human pressure imposed by bank top land cover" ~ "Bank-top human pressure" ))
view(morph_plot_data)
# checking the data
table(morph_plot_data$metric, useNA = "ifany")
Aquatic vegetation morphotypes Bank-top human pressure
10 10
Bed siltation Channel habitat complexity
10 10
Number of bed material types Number of flow types
10 10
Riparian habitat complexity Riparian vegetation complexity
10 10 morph_plot_data %>% select(metric, condition, module_position, value)# A tibble: 80 × 4
metric condition module_position value
<chr> <chr> <dbl> <chr>
1 Number of flow types Tree-lined 1 1
2 Number of flow types Tree-lined 2 1
3 Number of flow types Tree-lined 3 1
4 Number of flow types Tree-lined 4 1
5 Number of flow types Tree-lined 5 1
6 Number of flow types Non-tree-lined 1 2
7 Number of flow types Non-tree-lined 2 2
8 Number of flow types Non-tree-lined 3 2
9 Number of flow types Non-tree-lined 4 2
10 Number of flow types Non-tree-lined 5 1
# ℹ 70 more rows# calculating integrated physical habitat complexity by combining MoRPH indices
# 8 and 10
morph_physical <- morph_plot_data %>%
filter(index %in% c("INDEX 8 - Channel physical habitat complexity",
"INDEX 10 - Riparian physical habitat complexity")) %>%
mutate(value = as.numeric(value)) %>%
group_by(condition, module_position) %>%
summarise(integrated_complexity = sum(value),
.groups = "drop")
morph_physical# A tibble: 10 × 3
condition module_position integrated_complexity
<chr> <dbl> <dbl>
1 Non-tree-lined 1 5.58
2 Non-tree-lined 2 6
3 Non-tree-lined 3 5.33
4 Non-tree-lined 4 5.58
5 Non-tree-lined 5 6.83
6 Tree-lined 1 5.08
7 Tree-lined 2 5.32
8 Tree-lined 3 5.83
9 Tree-lined 4 6.42
10 Tree-lined 5 6.5 # mean for each riparian condition
morph_physical_summary <- morph_physical %>%
group_by(condition) %>%
summarise(mean = mean(integrated_complexity),
sd = sd(integrated_complexity),
min = min(integrated_complexity),
max = max(integrated_complexity))
morph_physical_summary# A tibble: 2 × 5
condition mean sd min max
<chr> <dbl> <dbl> <dbl> <dbl>
1 Non-tree-lined 5.87 0.591 5.33 6.83
2 Tree-lined 5.83 0.634 5.08 6.5 # calculating descriptive statistics for each MoRPH index and riparian condition
morph_summary <- morph_plot_data %>%
group_by(metric, condition) %>%
summarise(Mean = mean(value, na.rm = TRUE),
SD = sd(value, na.rm = TRUE),
Minimum = min(value, na.rm = TRUE),
Maximum = max(value, na.rm = TRUE),
.groups = "drop")
morph_summary# A tibble: 16 × 6
metric condition Mean SD Minimum Maximum
<chr> <chr> <dbl> <dbl> <chr> <chr>
1 Aquatic vegetation morphotypes Non-tree-lined NA 0.837 2 4
2 Aquatic vegetation morphotypes Tree-lined NA 0.548 1 2
3 Bank-top human pressure Non-tree-lined NA 0.822 3.5 5
4 Bank-top human pressure Tree-lined NA 0 3.5 3.5
5 Bed siltation Non-tree-lined NA 1.22 0.1 3.35
6 Bed siltation Tree-lined NA 1.58 3.55 7.65
7 Channel habitat complexity Non-tree-lined NA 0.456 2.083333333 3.3333…
8 Channel habitat complexity Tree-lined NA 0.618 1.25 2.9166…
9 Number of bed material types Non-tree-lined NA 0.837 2 4
10 Number of bed material types Tree-lined NA 0.707 1 3
11 Number of flow types Non-tree-lined NA 0.447 1 2
12 Number of flow types Tree-lined NA 0 1 1
13 Riparian habitat complexity Non-tree-lined NA 0.310 2.833333333 3.5
14 Riparian habitat complexity Tree-lined NA 0.414 3 4.0714…
15 Riparian vegetation complexity Non-tree-lined NA 0.354 3.5 4.5
16 Riparian vegetation complexity Tree-lined NA 1.30 4 7 # export MoRPH descriptive statistics table for appendix
# format table
morph_summary_table <- morph_summary %>%
rename("MoRPH index" = metric,
"Riparian condition" = condition) %>%
mutate(across(where(is.numeric), ~ round(.x, 2)))
morph_summary_table# A tibble: 16 × 6
`MoRPH index` `Riparian condition` Mean SD Minimum Maximum
<chr> <chr> <dbl> <dbl> <chr> <chr>
1 Aquatic vegetation morphoty… Non-tree-lined NA 0.84 2 4
2 Aquatic vegetation morphoty… Tree-lined NA 0.55 1 2
3 Bank-top human pressure Non-tree-lined NA 0.82 3.5 5
4 Bank-top human pressure Tree-lined NA 0 3.5 3.5
5 Bed siltation Non-tree-lined NA 1.22 0.1 3.35
6 Bed siltation Tree-lined NA 1.58 3.55 7.65
7 Channel habitat complexity Non-tree-lined NA 0.46 2.0833… 3.3333…
8 Channel habitat complexity Tree-lined NA 0.62 1.25 2.9166…
9 Number of bed material types Non-tree-lined NA 0.84 2 4
10 Number of bed material types Tree-lined NA 0.71 1 3
11 Number of flow types Non-tree-lined NA 0.45 1 2
12 Number of flow types Tree-lined NA 0 1 1
13 Riparian habitat complexity Non-tree-lined NA 0.31 2.8333… 3.5
14 Riparian habitat complexity Tree-lined NA 0.41 3 4.0714…
15 Riparian vegetation complex… Non-tree-lined NA 0.35 3.5 4.5
16 Riparian vegetation complex… Tree-lined NA 1.3 4 7 # format as flextable
morph_summary_table <- morph_summary_table %>%
flextable() %>%
theme_booktabs() %>%
autofit() %>%
align(align = "center", part = "all")
morph_summary_tableMoRPH index | Riparian condition | Mean | SD | Minimum | Maximum |
|---|---|---|---|---|---|
Aquatic vegetation morphotypes | Non-tree-lined | 0.84 | 2 | 4 | |
Aquatic vegetation morphotypes | Tree-lined | 0.55 | 1 | 2 | |
Bank-top human pressure | Non-tree-lined | 0.82 | 3.5 | 5 | |
Bank-top human pressure | Tree-lined | 0.00 | 3.5 | 3.5 | |
Bed siltation | Non-tree-lined | 1.22 | 0.1 | 3.35 | |
Bed siltation | Tree-lined | 1.58 | 3.55 | 7.65 | |
Channel habitat complexity | Non-tree-lined | 0.46 | 2.083333333 | 3.333333333 | |
Channel habitat complexity | Tree-lined | 0.62 | 1.25 | 2.916666667 | |
Number of bed material types | Non-tree-lined | 0.84 | 2 | 4 | |
Number of bed material types | Tree-lined | 0.71 | 1 | 3 | |
Number of flow types | Non-tree-lined | 0.45 | 1 | 2 | |
Number of flow types | Tree-lined | 0.00 | 1 | 1 | |
Riparian habitat complexity | Non-tree-lined | 0.31 | 2.833333333 | 3.5 | |
Riparian habitat complexity | Tree-lined | 0.41 | 3 | 4.071428571 | |
Riparian vegetation complexity | Non-tree-lined | 0.35 | 3.5 | 4.5 | |
Riparian vegetation complexity | Tree-lined | 1.30 | 4 | 7 |
# save to appendix
save_as_docx(morph_summary_table, path = "outputs/tables/A_morph_descriptive.docx")3. Descriptive spatial plots
# creating a plot for descriptive exploration of indices
# change values to numbers and round axis labels to 1 decimal place for visual purposes
morph_plot_data <- morph_plot_data %>%
mutate(value = as.numeric(value))
morph_spatial_plot <- ggplot(morph_plot_data,
aes(x = module_position,
y = value,
group = condition,
colour = condition)) +
geom_line() +
geom_point(size = 2) +
facet_wrap(~ metric,
scales = "free_y",
ncol = 2) +
scale_colour_manual(values = c("Tree-lined" = "#F8766D",
"Non-tree-lined" = "#00BFC4")) +
scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) +
labs(x = "Module",
y = "Index value",
colour = NULL ) +
theme_light() +
theme(legend.position = "bottom",
legend.margin = margin(t = -10),
legend.title = element_blank(),
panel.spacing = unit(0.3, "lines"))
morph_spatial_plot
# export plot
ggsave("outputs/figures/1_morph_spatial.png",
plot = morph_spatial_plot,
width = 18,
height = 15,
units = "cm",
dpi = 300)4. Summarising index data
morph_summary <- morph_plot_data %>%
group_by(metric, condition) %>%
summarise(mean = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
min = min(value, na.rm = TRUE),
max = max(value, na.rm = TRUE),
.groups = "drop")
view(morph_summary)
# check five modules in each riparian condition for every selected MoRPH index
table(morph_plot_data$metric, morph_plot_data$condition)
Non-tree-lined Tree-lined
Aquatic vegetation morphotypes 5 5
Bank-top human pressure 5 5
Bed siltation 5 5
Channel habitat complexity 5 5
Number of bed material types 5 5
Number of flow types 5 5
Riparian habitat complexity 5 5
Riparian vegetation complexity 5 5morph_plot_data %>%
group_by(metric, condition) %>%
summarise(n = n(),
.groups = "drop")# A tibble: 16 × 3
metric condition n
<chr> <chr> <int>
1 Aquatic vegetation morphotypes Non-tree-lined 5
2 Aquatic vegetation morphotypes Tree-lined 5
3 Bank-top human pressure Non-tree-lined 5
4 Bank-top human pressure Tree-lined 5
5 Bed siltation Non-tree-lined 5
6 Bed siltation Tree-lined 5
7 Channel habitat complexity Non-tree-lined 5
8 Channel habitat complexity Tree-lined 5
9 Number of bed material types Non-tree-lined 5
10 Number of bed material types Tree-lined 5
11 Number of flow types Non-tree-lined 5
12 Number of flow types Tree-lined 5
13 Riparian habitat complexity Non-tree-lined 5
14 Riparian habitat complexity Tree-lined 5
15 Riparian vegetation complexity Non-tree-lined 5
16 Riparian vegetation complexity Tree-lined 5# compare mean index values between riparian conditions
morph_comparison <- morph_summary %>%
select(metric, condition, mean) %>%
pivot_wider(names_from = condition,
values_from = mean) %>%
mutate(difference = `Non-tree-lined` - `Tree-lined`)
morph_comparison# A tibble: 8 × 4
metric `Non-tree-lined` `Tree-lined` difference
<chr> <dbl> <dbl> <dbl>
1 Aquatic vegetation morphotypes 2.8 1.4 1.4
2 Bank-top human pressure 4.1 3.5 0.600
3 Bed siltation 1.26 6.28 -5.02
4 Channel habitat complexity 2.58 2.17 0.417
5 Number of bed material types 2.8 2 0.8
6 Number of flow types 1.8 1 0.8
7 Riparian habitat complexity 3.28 3.66 -0.381
8 Riparian vegetation complexity 4 5.7 -1.7 # rounding values
morph_comparison <- morph_comparison %>%
mutate(`Non-tree-lined` = round(`Non-tree-lined`, 2),
`Tree-lined` = round(`Tree-lined`, 2),
difference = round(difference, 2))
# rename headers
morph_comparison <- morph_comparison %>%
rename(`MoRPH index` = metric,
Difference = difference)
morph_comparison# A tibble: 8 × 4
`MoRPH index` `Non-tree-lined` `Tree-lined` Difference
<chr> <dbl> <dbl> <dbl>
1 Aquatic vegetation morphotypes 2.8 1.4 1.4
2 Bank-top human pressure 4.1 3.5 0.6
3 Bed siltation 1.26 6.28 -5.02
4 Channel habitat complexity 2.58 2.17 0.42
5 Number of bed material types 2.8 2 0.8
6 Number of flow types 1.8 1 0.8
7 Riparian habitat complexity 3.28 3.66 -0.38
8 Riparian vegetation complexity 4 5.7 -1.7 # export table for main text
morph_comparison_table <- morph_comparison %>%
flextable() %>%
theme_booktabs() %>%
autofit() %>%
align(align = "center", part = "all")
save_as_docx(morph_comparison_table, path = "outputs/tables/morph_comparison.docx")