---
title: "Statistics for *Characteristics of macroplastic debris in eight diverse river systems across four continents*"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, results='hide', fig.keep='all', message=FALSE, warning=FALSE)
```

**Packages**
```{r}
library(tidyverse)
library(corrplot)
library(lmerTest)
library(lme4)
library(vegan)
library(DescTools)
```

**Datasets**

fig2_data.csv

fig2_glmm.csv

fig3_data.csv

fig4_data.csv

fig5_data.csv


### Figure 2. Total Debris & Plastic Debris Collected Statistics


**GLMM Setup**

- Categorical technology rank (1-3)
- Categorical waste picker activity rank (1-3)
- Numerical nearest urban population size
- Numerical measure of river size (width and/or length - discharge data unavailable)

1. Plastic Debris
2. Total Debris
3. Plastic Debris:Total Debris

*Check for collinearity*
```{r}
#Using the monthly means Figure 2 data to make a matrix for a correlation plot of the predictor variables

#Load data
fig2_data <- read_csv("fig2_data.csv")

#Creating data frame to use filtering for plastic debris to remove duplicates and keeping only the predictor variables
fig2_cor <- fig2_data %>%
  filter(
    debris_type == "plastic_debris"
  ) %>% 
  select(urban_population, tech_rank, waste_picker_rank, river_width_m, river_length)

#Making river length km from m
fig2_cor$river_length <- fig2_cor$river_length/1000

#Creating a matrix                       
fig2_cor <- as.matrix(fig2_cor)

#Making the correlation plot
cor(fig2_cor)
corrplot(cor(fig2_cor))
```

River length and river width appear to be highly collinear (>88%). Therefore, we will run two GLMMs for each debris scenario: one that includes river width, and one that includes river length, but not both together.

*Run GLMM scenarios*
```{r}
#Using the monthly totals Figure 2 data with the predictor variables included

#Load data
fig2_glmm <- read_csv("fig2_glmm.csv")

#Checking the variable classes in the data frame. Technology ranking and waste picker ranking are categorical variables and therefore should be factors
sapply(fig2_glmm, class)
fig2_glmm$tech_rank <- as.factor(fig2_glmm$tech_rank)
fig2_glmm$waste_picker_rank <- as.factor(fig2_glmm$waste_picker_rank)

#Z-Score the continuous variables because if not, we receive a warning that we should consider rescaling
fig2_glmm$urban_population <- scale(fig2_glmm$urban_population)
fig2_glmm$river_width_m <- scale(fig2_glmm$river_width_m)
fig2_glmm$river_length_km <- scale(fig2_glmm$river_length_km)

#Plastic debris scenario (adding -1 at the end of the model removes the fitting to an intercept)
plastic_debris_model_width <- lmer(plastic_debris ~ (1|date) + (1|river_id) + urban_population + tech_rank + waste_picker_rank + river_width_m-1, data=fig2_glmm) #river width
plastic_debris_model_length <- lmer(plastic_debris ~ (1|date) + (1|river_id) + urban_population + tech_rank + waste_picker_rank + river_length_km-1, data=fig2_glmm) #river length
summary(plastic_debris_model_width)
summary(plastic_debris_model_length)

#All debris scenario
all_debris_model_width <- lmer(all_debris ~ (1|date) + (1|river_id) + urban_population + tech_rank + waste_picker_rank + river_width_m-1, data=fig2_glmm) #river width
all_debris_model_length <- lmer(all_debris ~ (1|date) + (1|river_id) + urban_population + tech_rank + waste_picker_rank + river_length_km-1, data=fig2_glmm) #river length
summary(all_debris_model_width)
summary(all_debris_model_length)

#Proportion of plastic:all debris scenario
prop_debris_model_width <- lmer(plastic_debris/all_debris ~ (1|date) + (1|river_id) + urban_population + tech_rank + waste_picker_rank + river_width_m-1, data=fig2_glmm) #river width
prop_debris_model_length <- lmer(plastic_debris/all_debris ~ (1|date) + (1|river_id) + urban_population + tech_rank + waste_picker_rank + river_length_km-1, data=fig2_glmm) #river length
summary(prop_debris_model_width)
summary(prop_debris_model_length)

#There is a warning when running these models if scale function not added to continuous variables before running: "some predictor variables are on very different scales: consider rescaling"

plastic_debris <- lmer(plastic_debris ~ (1|date) + river_id, data=fig2_glmm)
summary(plastic_debris)

anova <- aov(plastic_debris ~ river_id, data = fig2_glmm)
summary(anova)
```

**GLMM results**

Across all models, there are no significant effects from any of our predictor variables.

*Plastic Debris LMER (river width)*
```{r, results='markup'}
summary(plastic_debris_model_width)
```

*Plastic Debris LMER (river length)*
```{r, results='markup'}
summary(plastic_debris_model_length)
```

*All Debris LMER (river width)*
```{r, results='markup'}
summary(all_debris_model_width)
```

*All Debris LMER (river length)*
```{r, results='markup'}
summary(all_debris_model_length)
```

*Proportion Debris LMER (river width)*
```{r, results='markup'}
summary(prop_debris_model_width)
```

*Proportion Debris LMER (river length)*
```{r, results='markup'}
summary(prop_debris_model_length)
```

**ANOVA Setup**

- Want to explore whether plastic and total debris are significantly different across sites
- Option 1: ANOVA by location
- Option 2: LMER with date as random effect, pairwise test and compare across all sites where the reference level is set for each river

- Of the options, the ANOVA is the easiest to interpret and will be the most natural fit, and will therefore be using Option 1

1. Plastic Debris
2. Total Debris

*Run ANOVA scenarios*
```{r}
#ANOVA for plastic debris by river
anova_plastic <- aov(plastic_debris ~ river_id, data = fig2_glmm)
summary(anova_plastic)

#ANOVA for all debris by river
anova_debris <- aov(all_debris ~ river_id, data = fig2_glmm)
summary(anova_debris)
```

**ANOVA results**

*Plastic Debris ANOVA*
```{r, results='markup'}
summary(anova_plastic)
EtaSq(anova_plastic)

# Extract the p-value
p_value_plastic <- summary(anova_plastic)[[1]][["Pr(>F)"]][1]
print(p_value_plastic)
```

Plastic debris collected across sites differs significantly (p = 1.08712e-44).

*All Debris ANOVA*
```{r, results='markup'}
summary(anova_debris)
EtaSq(anova_debris)

# Extract the p-value
p_value_debris <- summary(anova_debris)[[1]][["Pr(>F)"]][1]
print(p_value_debris)
```

Total debris collected across sites differs significantly (p = 1.353249e-45).


### Figure 3. Polymer Composition of Plastic Debris Collected

### Figure 4. Proportion of Single-Use Plastic Item Categories in Plastic Debris Collected

### Figure 5. End-of-Life Fate of Plastic Debris Collected


**Chi Square Setup**

- Testing the null hypothesis that the proportions of these variables should be equal across sites
- A significant P value would indicate that these proportions are different across sites

1. Polymer Proportions across sites
2. Single-Use Item Proportions across subset of sites
3. End-of-Life Fate Proportions across sites

*Prepare Chi Square contingency tables*
```{r}
#Using proportion data for Figures 3, 4, and 5

#Load data
fig3_data <- read_csv("fig3_data.csv")
fig4_data <- read_csv("fig4_data.csv")
fig5_data <- read_csv("fig5_data.csv")

#Creating a contingency table of polymer proportions
fig3_data_contingency <- fig3_data %>%
  select(
    country, polymer, polymer_proportions
  ) %>% 
  pivot_wider(
    names_from = "polymer",
    values_from = "polymer_proportions"
  ) %>% 
  select(
    !country
  )

#Multiplying the proportions by 100 (%)
fig3_data_contingency <- round(mutate_if(fig3_data_contingency, is.numeric, ~ . * 100),1)

#Creating a contingency table of single-use item proportions
fig4_data_contingency <- fig4_data %>%
  select(
    country, item, item_proportions
  ) %>% 
  pivot_wider(
    names_from = "item",
    values_from = "item_proportions"
  ) %>% 
  select(
    !country
  )

#Multiplying the proportions by 100 (%)
fig4_data_contingency <- round(mutate_if(fig4_data_contingency, is.numeric, ~ . * 100),1)

#Creating a contingency table of end-of-life fate proportions
fig5_data_contingency <- fig5_data %>%
  select(
    country, fate, fate_proportions
  ) %>% 
  pivot_wider(
    names_from = "fate",
    values_from = "fate_proportions"
  )  %>% 
  select(
    !country
  )

#Multiplying the proportions by 100 (%)
fig5_data_contingency <- round(mutate_if(fig5_data_contingency, is.numeric, ~ . * 100),1)
```

*Run Chi Square tests*
```{r}
fig3_chi <- chisq.test(fig3_data_contingency)
fig4_chi <- chisq.test(fig4_data_contingency)
fig5_chi <- chisq.test(fig5_data_contingency)
```

**Chi Square results**

*Polymer Proportion Chi Square*
```{r, results='markup'}
fig3_chi
```

The polymer proportions across sites are significantly different (p = 9.388124e-131).

*Single-Use Item Proportion Chi Square*
```{r, results='markup'}
fig4_chi
```

The single-use item proportions across sites are significantly different (p = 7.86554e-21).

*End-of-Life Fate Proportion Chi Square*
```{r, results='markup'}
fig5_chi
```

The end-of-life fate proportions across sites are significantly different (p = 4.940656e-324)
