---
title: "Code to reproduce analyses"
author:
- Veronica Price-Jones
- Lien Reyserhove
date: "`r Sys.Date()`"
output:
  pdf_document:
    toc: yes
    toc_depth: '3'
  html_document:
    df_print: paged
    number_sections: yes
    toc: yes
    toc_depth: 3
    toc_float: yes
subtitle: ''
---

(short description of the goal of this document)

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

Load libraries:

```{r}
library(tidyverse)      # To do data science
library(magrittr)       # To use %<>% pipes
library(here)           # To find files
library(readxl)         # To read Excel files
library(janitor)        # To clean input data
library(ordinal)        # To perform clm
library(dplyr)          # To clean input data
library(splitstackshape)# To split rows for Chisq test
```

# Read source data

Create a data frame `input_data` from the source data:
 
```{r}
input_data <- read_csv(file = here::here("data", "processed_data.csv")) 
```

# Clean data

Remove empty rows (if present) and clean column headers:

```{r}
input_data %<>%
  remove_empty("rows") %>%       # Remove empty rows
  clean_names()                  # Have sensible (lowercase) column names
```

Show cleaned column names:

```{r}
colnames(input_data)
```

Preview data:

```{r}
input_data %>% head(n = 5)
```

# Fisher's Exact tests

## Audience/Taxa

Create data frame:

```{r}
df_aud_tax = data.frame(input_data$audience, input_data$taxa)
df_aud_tax
```

Split rows with multiple answers:

```{r}
df_aud_tax <- cSplit(df_aud_tax, "input_data.audience", sep=",", direction="long")
df_aud_tax <- cSplit(df_aud_tax, "input_data.taxa", sep=",", direction="long")
```

Fisher's exact test:

```{r}
fisher.test(df_aud_tax$input_data.audience, df_aud_tax$input_data.taxa,
            simulate.p.value=TRUE)
```

## Audience/Environment

Create data frame:

```{r}
df_aud_env = data.frame(input_data$audience, input_data$environment)
df_aud_env
```

Split rows with multiple answers:

```{r}
df_aud_env <- cSplit(df_aud_env, "input_data.audience", sep=",", direction="long")
df_aud_env <- cSplit(df_aud_env, "input_data.environment", sep=",", direction="long")
```

Fisher's exact test:

```{r}
fisher.test(df_aud_env$input_data.audience, df_aud_env$input_data.environment,
            simulate.p.value=TRUE)
```

## Region/Audience

Create data frame:

```{r}
df_reg_aud = data.frame(input_data$region, input_data$audience)
df_reg_aud
```

Remove rows without region:

```{r}
df_reg_aud <- df_reg_aud[-c(11, 29, 58),]
```

Split rows with multiple answers:

```{r}
df_reg_aud <- cSplit(df_reg_aud, "input_data.region", sep=",", direction="long")
df_reg_aud <- cSplit(df_reg_aud, "input_data.audience", sep=",", direction="long")
```

Fisher's exact test:

```{r}
fisher.test(df_reg_aud$input_data.region, df_reg_aud$input_data.audience,
            simulate.p.value=TRUE)
```

## Region/Taxa

Create data frame:

```{r}
df_reg_tax = data.frame(input_data$region, input_data$taxa)
df_reg_tax
```

Remove rows without region:

```{r}
df_reg_tax <- df_reg_tax[-c(11, 29, 58),]
```

Split rows with multiple answers:

```{r}
df_reg_tax <- cSplit(df_reg_tax, "input_data.region", sep=",", direction="long")
df_reg_tax <- cSplit(df_reg_tax, "input_data.taxa", sep=",", direction="long")
```

Fisher's exact test:

```{r}
fisher.test(df_reg_tax$input_data.region, df_reg_tax$input_data.taxa,
            simulate.p.value=TRUE)
```

## Region/Environment

Create data frame:

```{r}
df_reg_env = data.frame(input_data$region, input_data$environment)
df_reg_env
```

Remove rows without region:

```{r}
df_reg_env <- df_reg_env[-c(11, 29, 58),]
```

Split rows with multiple answers:

```{r}
df_reg_env <- cSplit(df_reg_env, "input_data.region", sep=",", direction="long")
df_reg_env <- cSplit(df_reg_env, "input_data.environment", sep=",", direction="long")
```

Fisher's exact test:

```{r}
fisher.test(df_reg_env$input_data.region, df_reg_env$input_data.environment,
            simulate.p.value=TRUE)
```

Chi square test:

```{r}
chisq.test(df_reg_env$input_data.region, df_reg_env$input_data.environment)
```

# Culmulative Links Models

## Pre-processing

### Number of invasive alien records

Rename variable `number_of_alien_or_invasive_alien_species_records` to `no_ias_records` to simplify the analysis:

```{r}
input_data <- input_data %>% 
  rename(no_ias_records = number_of_alien_or_invasive_alien_species_records) 
```

Inspect different values of `no_ias_records`:

```{r}
table(input_data$no_ias_records, useNA = "always")
```

Change `Info not available` and `Not applicable` to NA:

```{r}
input_data$no_ias_records <- na_if(input_data$no_ias_records, "Info not available") 
input_data$no_ias_records <- na_if(input_data$no_ias_records, "Not applicable") 
input_data$no_ias_records <- na_if(input_data$no_ias_records, "NA") 
```

Define `no_ias_records` as an ordered factor:

```{r}
input_data$no_ias_records <- factor(input_data$no_ias_records, ordered = TRUE)
```

Re-order `no_ias_records`:

```{r}
input_data$no_ias_records <- 
  ordered(input_data$no_ias_records, 
          levels = c("1 - 50", 
                     "51 - 100",
                     "101 - 500", 
                     "501 - 1000", 
                     "1001 - 2000", 
                     "2001 - 5000", 
                     "5001 - 10000",
                     "10001 - 50000", 
                     "50001 - 100000", 
                     "100001 - 500000",
                     "500001 - 1000000",
                     "> 1000001"))
```

Verify order:

```{r}
levels(input_data$no_ias_records)
```

### Number of participants

Inspect different values of `number_of_participants`:

```{r}
table(input_data$number_of_participants, useNA = "always")
```

Change `Info not available` and `Not applicable` to NA:

```{r}
input_data$number_of_participants <- na_if(input_data$number_of_participants, "Info not available") 
input_data$number_of_participants <- na_if(input_data$number_of_participants, "Not applicable") 
input_data$number_of_participants <- na_if(input_data$number_of_participants, "NA") 
```

Define `number_of_participants` as an ordered factor:

```{r}
input_data$number_of_participants <- factor(input_data $ number_of_participants, ordered = TRUE)
```

Re-order `number_of_participants`:

```{r}
input_data$number_of_participants <- 
  ordered (input_data$number_of_participants, 
           levels = c("1 - 50",
                      "51 - 100",
                      "101 - 500",
                      "501 - 1000",
                      "1001 - 2000",
                      "2001 - 5000",
                      "5001 - 10000",
                      "> 10000"))
```

Verify order:

```{r}
levels(input_data$number_of_participants)
```


### Number of publications

Inspect different values of `number_of_publications`:

```{r}
table(input_data$number_of_publications, useNA = "always")
```

Change `Unknown` and `Not applicable` to NA:

```{r}
input_data$number_of_publications <- na_if(input_data$number_of_publications, "Unknown") 
input_data$number_of_publications <- na_if(input_data$number_of_publications, "Not applicable") 
input_data$number_of_publications <- na_if(input_data$number_of_publications, "NA") 
```

Define `number_of_publications` as an ordered factor:

```{r}
input_data$number_of_publications <- factor(input_data $ number_of_publications, ordered = TRUE)
```

Re-order `number_of_participants`:

```{r}
input_data$number_of_publications <- 
  ordered (input_data$number_of_publications, 
           levels = c("0", 
                      "1 to 5", 
                      "6 to 10", 
                      "11 to 20", 
                      "21 - 50", 
                      "> 100"))
```

Verify order:

```{r}
str(input_data$number_of_publications)
```



## Analysis

### Inspect independent variables

The following independent variables were selected for the analyses:

- duration (`duration`)
- project design (`design`)
- use of social media (`sm_yn`)
- level of skill required (`skills`)
- expected contribution frequencey (`contribution_frequency`)
- provision of guidelines (`guidelines`)
- training (`training_yn`)
- sighting maps (`map`)
- active informing (`active_informing`)
- feedback (`feedback`)
- support (`support`) 

Inspect dependent variables on percentage of `NA`:

```{r}
input_data %>% 
  select(duration, design, sm_yn,  skills, contribution_frequency, guidelines, training_yn, map, active_informing, feedback, support) %>% 
  map(~mean(is.na(.)))
```

### Number of invasive alien records

The following variables are specified for analysis:

- `design`
- `duration`
- `sm_yn`
- `skills`
- `contribution_frequency`
- `guidelines`
- `training_yn`
- `map`
- `active_informing`
- `feedback`
- `support`

Cumulative Link Model:

```{r}
clm_ias_records <- 
  clm(
    no_ias_records ~ 
      duration + design + sm_yn + 
      skills + contribution_frequency + guidelines + 
      training_yn + map + active_informing + feedback + 
      support,
    data = input_data)
```

Model output:

```{r}
summary(clm_ias_records)
```

Output with p-values:

```{r}
anova(clm_ias_records)
```

Conclusion:
- significant positive effect of duration and provision of map
- significant negative effect of limited/no prior skills/knowledge and provision of feedback
- 64 deleted observations

### Number of participants

Cumulative Link Model

The following variables are specified for analysis:

- `design`
- `duration`
- `sm_yn`
- `skills`
- `contribution_frequency`
- `guidelines`
- `training_yn`
- `map`
- `active_informing`
- `feedback`
- `support`

Cumulative link model:

```{r}
clm_participants <- 
  clm(
    number_of_participants ~ duration + design + sm_yn + skills + contribution_frequency + guidelines + training_yn + map + active_informing + feedback + support,
    data = input_data)
```

Model output:

```{r}
summary(clm_participants)
```  

Output with p-values:

```{r}
anova(clm_participants)
```

Conclusion: 
- significant positive effect of duration and map
- 61 deleted observations

### Number of publications

Cumulative Link Model

The following variables are specified for analysis:

- `design`
- `duration`
- `sm_yn`
- `skills`
- `contribution_frequency`
- `guidelines`
- `training_yn`
- `map`
- `active_informing`
- `feedback`
- `support`

Cumulative link model:

```{r}
clm_publications <- 
  clm(
    number_of_publications ~ duration + design + sm_yn + skills + contribution_frequency + guidelines + training_yn + map + active_informing + feedback + support,
    data = input_data)
```

Model output:

```{r}
summary(clm_publications)
```  

Output with p-values:

```{r}
anova(clm_publications)
```

Conclusion:
- significant positive effect of duration, social media, guidelines, training, map and feedback
-significant negative effect of limited/no prior skills, support
-59 deleted observations
