!!!!!!!!!!!!!
!!!READ ME!!!
!!!!!!!!!!!!!

This code is designed to be self-sufficient and requiring /minimal/ changes. 
Changes outside of the specified areas will lead to problems!
Before starting the code below:
1) verify all libraries are installed (in the console below) prior to running Setup
2) /carefully/ read the annotations 

#Setting-up and reading in data

```{r Setup, include = FALSE, echo = FALSE}

#libraries
library(devtools)
library(readxl)
library(tidyverse)
library(magrittr)
library(plotrix)
library(cowplot)
options(scipen = 999)

set.seed(666) #Do NOT change this :) 

#Import data
df <- read_xlsx(file.choose())
df <- data.frame(df)

#In df, there's >1 column for Uninvaded_measure, >1 column for Managed_measure; need to stack related variables
tidy.df <- df %>%
  select(-Entered, -Screened, -Title, -Invaded_SEM, -Invaded_CI, -Uninvaded_SEM, - Uninvaded_SEM2, -Managed_SEM, -Managed_SEM2, -Notes) %>%
  gather(., key = "Uninvaded_Type", value = "Uninvaded_measure", Uninvaded_measure, Uninvaded_measure2) %>%
  gather(., key = "UninvadedSD_Type", value = "Uninvaded_SD", Uninvaded_SD, Uninvaded_SD2) %>%
  gather(., key = "Managed_Type", value = "Managed_measure", Managed_measure, Managed_measure2) %>%
  gather(., key = "ManagedSD_Type", value = "Managed_SD", Managed_SD, Managed_SD2) %>%
  mutate(citation = paste(Author, Date)) %>% #Creating custom column for citation
  mutate(Binned_Lat = ifelse(Latitude < -35, -40,
                             ifelse(Latitude < -30, -35,
                                    ifelse(Latitude < -25, -30,
                                           ifelse(Latitude < -20, -25,
                                                  ifelse(Latitude < -15, -20,
                                                         ifelse(Latitude < -10, -15,
                                                                ifelse(Latitude < -5, -10,
                                                                       ifelse(Latitude < 0, -5,
                                                                              ifelse(Latitude < 5, 0,
                                                                                     ifelse(Latitude < 10, 5,
                                                                                            ifelse(Latitude < 15, 10,
                                                                                                   ifelse(Latitude < 20, 15,
                                                                                                          ifelse(Latitude < 25, 20,
                                                                                                                 ifelse(Latitude < 30, 25,
                                                                                                                        ifelse(Latitude < 35, 30,
                                                                                                                               ifelse(Latitude < 40, 35,
                                                                                                                                      ifelse(Latitude < 45, 40, 45)))))))))))))))))) %>%
  mutate(Binned_Lat = as.factor(Binned_Lat)) %>% #Randoms /cannot/ be numeric!!!!!!
  mutate_if(is.character, as.factor) #setting up as factors
  
```

#Simplify your life
Filter 1 = "abiotic" or "biotic"
Filter 2 = Metric_Type of "carbon", "nitrogen", "other nutrient", "soil characteristic", or "abundance"
Filter 3 = Finer_Metric_Type, Organism_Type, or Wetland_Classification

```{r Establish filters}

Filter1 <- "abiotic"
Filter2 <- "carbon"
colnames(tidy.df)[colnames(tidy.df) == "Wetland_Classification"] <- "Filter3" #CHANGE FINER_METRIC_TYPE BASED ON YOUR THIRD FILTER
Filter3 <- "Wetland_Classification"

```

#Starting analysis
PART A: Look at uninvaded-invaded ONLY
STEP 1: Apply filter (make edits *only* as indicated, read annotations and text)
STEP 2: Permute
STEP 3: Save permute results

PART B: Look at managed-invaded ONLY
STEP 1: Permute
STEP 2: Save permtue results

PART C: Graphing results (together)

#Starting Part A: uninvaded-invaded
Planning

```{r Planning for Part A}
#How many chunks will you need?
PlanA.Table <- tidy.df %>%
  filter(Biotic_or_Abiotic == Filter1) %>% 
  filter(Metric_Type == Filter2) %>% 
  filter(Uninvaded_measure != "NA") %>%
  filter(Uninvaded_SD != "NA") %>%
  filter(Uninvaded_SD > 0) %>%
  filter(Invaded_measure != "NA") %>%
  filter(Invaded_SD != "NA") %>%
  filter(Invaded_SD > 0) %>%
  group_by(Filter3) %>%
  summarise(n = length(Uninvaded_measure-Invaded_measure)) #This will show you the number of comparisons that can be made per group

PlanA.Table
#The number of categories listed will indicate the number of analysis chunks you will need below. Here in this example, there are five categories (c flux, c pool, mass, microbial, and organic) so analysis Part A will require 5 chunks.

#Do yourself a favour and define the levels as filters, formatted this way. If you have <5 levels, delete the unusued/unneeded lines; if you have >5 levels, continue with the same formatting. Recall that the text you enter for your levels must be formatted the same way (capitalisation, spacing).
Filter.PartA.1 <- "estuary"
Filter.PartA.2 <- "marsh"
Filter.PartA.3 <- "swamp"
#Filter.PartA.4 <- "microbial"
#Filter.PartA.5 <- "organic"
#Filter.PartA.6 <- "biomass ratio"
#Filter.PartA.7 <- "chlorophyll"
#Filter.PartA.8 <- "density"
#Filter.PartA.9 <- "invert biomass"
#Filter.PartA.10 <- "microbial"
```

#Part A

```{r Part A.1}

PartA.1.set <- tidy.df %>%
  filter(Biotic_or_Abiotic == Filter1) %>% 
  filter(Metric_Type == Filter2) %>% 
  filter(Filter3 == Filter.PartA.1) %>% 
  filter(Uninvaded_measure != "NA") %>%
  filter(Uninvaded_SD != "NA") %>%
  filter(Uninvaded_Sample.size != "NA") %>%
  filter(Invaded_measure != "NA") %>%
  filter(Invaded_SD != "NA") %>%
  filter(Invaded_Sample.size != "NA") %>%
  filter(Uninvaded_SD > 0) %>%
  filter(Invaded_SD > 0) %>%
  mutate(SDPool = sqrt((((Uninvaded_Sample.size-1)*Uninvaded_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Uninvaded_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Uninvaded_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Uninvaded_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))

##Saving this "real difference" as the object PartA.1.set.testobs, which we will use later
PartA.1.set.testobs <- mean(PartA.1.set$Hedges.G)

##Getting ready for randomisations
reps <- 10000 #Tells R we want 10,000 permutations
PartA.1.set.reps <- numeric(reps)

PartA.1.set %<>%
  select(-SDPool, -N, -Hedges.G) %>%
  gather(., key = "Measure_status", value = "Measure", Uninvaded_measure, Invaded_measure) %>%
  gather(., key = "SD_status", value = "SD", Uninvaded_SD, Invaded_SD) %>%
  gather(., key = "n_status", value = "n", Uninvaded_Sample.size, Invaded_Sample.size) %>%
  mutate(status = ifelse(Measure_status == "Uninvaded_measure", "Uninvaded", "Invaded")) %>%
  select(-Measure_status, -SD_status, -n_status)

for(i in 1:reps) {
  PartA.1.set$statusSHUFFLE <- sample(PartA.1.set$status) #shuffles the invasion status variable
  PartA.1.seta <- PartA.1.set %>%
    filter(statusSHUFFLE == "Uninvaded")
  PartA.1.setb <- PartA.1.set %>%
    filter(statusSHUFFLE == "Invaded")
  colnames(PartA.1.seta)[colnames(PartA.1.seta)=="Measure"] <- "Uninvaded_measure"
  colnames(PartA.1.seta)[colnames(PartA.1.seta)=="SD"] <- "Uninvaded_SD"
  colnames(PartA.1.seta)[colnames(PartA.1.seta)=="n"] <- "Uninvaded_Sample.size"
  PartA.1.seta$Invaded_measure <- PartA.1.setb$Measure
  PartA.1.seta$Invaded_SD <- PartA.1.setb$SD
  PartA.1.seta$Invaded_Sample.size <- PartA.1.setb$n
  PartA.1.seta %<>%
  mutate(SDPool = sqrt((((Uninvaded_Sample.size-1)*Uninvaded_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Uninvaded_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Uninvaded_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Uninvaded_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))
  PartA.1.set.reps[i] <- mean(PartA.1.seta$Hedges.G) #Take Hedges.G and save it as a permuted difference in teststat
}

##Generate a histogram of the permuted mean differences, add vertical line of our observed difference
ggplot() + geom_histogram(aes(x = PartA.1.set.reps), binwidth = 0.01) + geom_vline(aes(xintercept = PartA.1.set.testobs), colour = "red", linetype = "dashed", size = 1) + xlab("Permuted Hedges' G") + ylab("Frequency") + theme_cowplot()

##Generate the p-value (two-tailed)
options(digits = 3) #Adjust to see appropriate number of sig-digs
(sum(abs(PartA.1.set.reps) >= abs(PartA.1.set.testobs)) + 1)/(sum(reps) + 1) #p-value

#Step 3: Save permuted results
PartA.1.set.reps <- as.data.frame(PartA.1.set.reps)
fig.PartA.1.bindr <- PartA.1.set.reps %>%
  summarise(avg = mean(PartA.1.set.reps),
            lwr = avg - std.error(PartA.1.set.reps)*1.96,
            upr = avg + std.error(PartA.1.set.reps)*1.96)%>%
  add_column(obs.diff = PartA.1.set.testobs,
             var = Filter.PartA.1,                 
             p = (sum(abs(PartA.1.set.reps) >= abs(PartA.1.set.testobs)) + 1)/(sum(reps) + 1))

```

#Repeat Part A as needed
Copy code in previous chunk, paste into a text editor (*NOT* microsoft word, but notepad or an equivalent), search-replace PartA.1 and replace with PartA.2 for next step. Be sure to change the third filter to the next level of the third filter you are analysing!

```{r Part A.2}

PartA.2.set <- tidy.df %>%
  filter(Biotic_or_Abiotic == Filter1) %>% 
  filter(Metric_Type == Filter2) %>% 
  filter(Filter3 == Filter.PartA.2) %>% 
  filter(Uninvaded_measure != "NA") %>%
  filter(Uninvaded_SD != "NA") %>%
  filter(Uninvaded_Sample.size != "NA") %>%
  filter(Invaded_measure != "NA") %>%
  filter(Invaded_SD != "NA") %>%
  filter(Invaded_Sample.size != "NA") %>%
  filter(Uninvaded_SD > 0) %>%
  filter(Invaded_SD > 0) %>%
  mutate(SDPool = sqrt((((Uninvaded_Sample.size-1)*Uninvaded_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Uninvaded_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Uninvaded_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Uninvaded_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))

##Saving this "real difference" as the object PartA.2.set.testobs, which we will use later
PartA.2.set.testobs <- mean(PartA.2.set$Hedges.G)

##Getting ready for randomisations
reps <- 10000 #Tells R we want 10,000 permutations
PartA.2.set.reps <- numeric(reps)

PartA.2.set %<>%
  select(-SDPool, -N, -Hedges.G) %>%
  gather(., key = "Measure_status", value = "Measure", Uninvaded_measure, Invaded_measure) %>%
  gather(., key = "SD_status", value = "SD", Uninvaded_SD, Invaded_SD) %>%
  gather(., key = "n_status", value = "n", Uninvaded_Sample.size, Invaded_Sample.size) %>%
  mutate(status = ifelse(Measure_status == "Uninvaded_measure", "Uninvaded", "Invaded")) %>%
  select(-Measure_status, -SD_status, -n_status)

for(i in 1:reps) {
  PartA.2.set$statusSHUFFLE <- sample(PartA.2.set$status) #shuffles the invasion status variable
  PartA.2.seta <- PartA.2.set %>%
    filter(statusSHUFFLE == "Uninvaded")
  PartA.2.setb <- PartA.2.set %>%
    filter(statusSHUFFLE == "Invaded")
  colnames(PartA.2.seta)[colnames(PartA.2.seta)=="Measure"] <- "Uninvaded_measure"
  colnames(PartA.2.seta)[colnames(PartA.2.seta)=="SD"] <- "Uninvaded_SD"
  colnames(PartA.2.seta)[colnames(PartA.2.seta)=="n"] <- "Uninvaded_Sample.size"
  PartA.2.seta$Invaded_measure <- PartA.2.setb$Measure
  PartA.2.seta$Invaded_SD <- PartA.2.setb$SD
  PartA.2.seta$Invaded_Sample.size <- PartA.2.setb$n
  PartA.2.seta %<>%
  mutate(SDPool = sqrt((((Uninvaded_Sample.size-1)*Uninvaded_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Uninvaded_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Uninvaded_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Uninvaded_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))
  PartA.2.set.reps[i] <- mean(PartA.2.seta$Hedges.G) #Take Hedges.G and save it as a permuted difference in teststat
}

##Generate a histogram of the permuted mean differences, add vertical line of our observed difference
ggplot() + geom_histogram(aes(x = PartA.2.set.reps), binwidth = 0.01) + geom_vline(aes(xintercept = PartA.2.set.testobs), colour = "red", linetype = "dashed", size = 1) + xlab("Permuted Hedges' G") + ylab("Frequency") + theme_cowplot()

##Generate the p-value (two-tailed)
options(digits = 3) #Adjust to see appropriate number of sig-digs
(sum(abs(PartA.2.set.reps) >= abs(PartA.2.set.testobs)) + 1)/(sum(reps) + 1) #p-value

#Step 3: Save permuted results
PartA.2.set.reps <- as.data.frame(PartA.2.set.reps)
fig.PartA.2.bindr <- PartA.2.set.reps %>%
  summarise(avg = mean(PartA.2.set.reps),
            lwr = avg - std.error(PartA.2.set.reps)*1.96,
            upr = avg + std.error(PartA.2.set.reps)*1.96)%>%
  add_column(obs.diff = PartA.2.set.testobs,
             var = Filter.PartA.2,                 
             p = (sum(abs(PartA.2.set.reps) >= abs(PartA.2.set.testobs)) + 1)/(sum(reps) + 1))

```

#Repeat Part A as needed
Copy code in previous chunk, paste into a text editor (*NOT* microsoft word, but notepad or an equivalent), search-replace PartA.1 and replace with PartA.2 for next step. Be sure to change the third filter to the next level of the third filter you are analysing!

```{r Part A.3}

PartA.3.set <- tidy.df %>%
  filter(Biotic_or_Abiotic == Filter1) %>% 
  filter(Metric_Type == Filter2) %>% 
  filter(Filter3 == Filter.PartA.3) %>% 
  filter(Uninvaded_measure != "NA") %>%
  filter(Uninvaded_SD != "NA") %>%
  filter(Uninvaded_Sample.size != "NA") %>%
  filter(Invaded_measure != "NA") %>%
  filter(Invaded_SD != "NA") %>%
  filter(Invaded_Sample.size != "NA") %>%
  filter(Uninvaded_SD > 0) %>%
  filter(Invaded_SD > 0) %>%
  mutate(SDPool = sqrt((((Uninvaded_Sample.size-1)*Uninvaded_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Uninvaded_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Uninvaded_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Uninvaded_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))

##Saving this "real difference" as the object PartA.3.set.testobs, which we will use later
PartA.3.set.testobs <- mean(PartA.3.set$Hedges.G)

##Getting ready for randomisations
reps <- 10000 #Tells R we want 10,000 permutations
PartA.3.set.reps <- numeric(reps)

PartA.3.set %<>%
  select(-SDPool, -N, -Hedges.G) %>%
  gather(., key = "Measure_status", value = "Measure", Uninvaded_measure, Invaded_measure) %>%
  gather(., key = "SD_status", value = "SD", Uninvaded_SD, Invaded_SD) %>%
  gather(., key = "n_status", value = "n", Uninvaded_Sample.size, Invaded_Sample.size) %>%
  mutate(status = ifelse(Measure_status == "Uninvaded_measure", "Uninvaded", "Invaded")) %>%
  select(-Measure_status, -SD_status, -n_status)

for(i in 1:reps) {
  PartA.3.set$statusSHUFFLE <- sample(PartA.3.set$status) #shuffles the invasion status variable
  PartA.3.seta <- PartA.3.set %>%
    filter(statusSHUFFLE == "Uninvaded")
  PartA.3.setb <- PartA.3.set %>%
    filter(statusSHUFFLE == "Invaded")
  colnames(PartA.3.seta)[colnames(PartA.3.seta)=="Measure"] <- "Uninvaded_measure"
  colnames(PartA.3.seta)[colnames(PartA.3.seta)=="SD"] <- "Uninvaded_SD"
  colnames(PartA.3.seta)[colnames(PartA.3.seta)=="n"] <- "Uninvaded_Sample.size"
  PartA.3.seta$Invaded_measure <- PartA.3.setb$Measure
  PartA.3.seta$Invaded_SD <- PartA.3.setb$SD
  PartA.3.seta$Invaded_Sample.size <- PartA.3.setb$n
  PartA.3.seta %<>%
  mutate(SDPool = sqrt((((Uninvaded_Sample.size-1)*Uninvaded_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Uninvaded_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Uninvaded_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Uninvaded_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))
  PartA.3.set.reps[i] <- mean(PartA.3.seta$Hedges.G) #Take Hedges.G and save it as a permuted difference in teststat
}

##Generate a histogram of the permuted mean differences, add vertical line of our observed difference
ggplot() + geom_histogram(aes(x = PartA.3.set.reps), binwidth = 0.01) + geom_vline(aes(xintercept = PartA.3.set.testobs), colour = "red", linetype = "dashed", size = 1) + xlab("Permuted Hedges' G") + ylab("Frequency") + theme_cowplot()

##Generate the p-value (two-tailed)
options(digits = 3) #Adjust to see appropriate number of sig-digs
(sum(abs(PartA.3.set.reps) >= abs(PartA.3.set.testobs)) + 1)/(sum(reps) + 1) #p-value

#Step 3: Save permuted results
PartA.3.set.reps <- as.data.frame(PartA.3.set.reps)
fig.PartA.3.bindr <- PartA.3.set.reps %>%
  summarise(avg = mean(PartA.3.set.reps),
            lwr = avg - std.error(PartA.3.set.reps)*1.96,
            upr = avg + std.error(PartA.3.set.reps)*1.96)%>%
  add_column(obs.diff = PartA.3.set.testobs,
             var = Filter.PartA.3,                 
             p = (sum(abs(PartA.3.set.reps) >= abs(PartA.3.set.testobs)) + 1)/(sum(reps) + 1))

fig.PartA.3.bindr

```

#Starting Part B: managed-invaded

```{r Planning for Part B}

#How many chunks will you need for Part B?
PlanB.Table <- tidy.df %>%
  filter(Biotic_or_Abiotic == Filter1) %>% 
  filter(Metric_Type == Filter2) %>% 
  filter(Managed_measure != "NA") %>%
  filter(Managed_SD != "NA") %>%
  filter(Managed_SD > 0) %>%
  filter(Invaded_measure != "NA") %>%
  filter(Invaded_SD != "NA") %>%
  filter(Invaded_SD > 0) %>%
  group_by(Filter3) %>%
  summarise(n = length(Managed_measure-Invaded_measure)) #This will show you the number of comparisons that can be made per group

PlanB.Table
#The number of categories listed will indicate the number of analysis chunks you will need below. Here in this example, there is only one group where managed data exists (c flux), thus only one PartB code chunk is needed.

#You do not need to re-define objects

```

#Part B

```{r Part B.1}

PartB.1.set <- tidy.df %>%
  filter(Biotic_or_Abiotic == Filter1) %>% 
  filter(Metric_Type == Filter2) %>% 
  filter(Filter3 == Filter.PartA.1) %>% 
  filter(Managed_measure != "NA") %>%
  filter(Managed_SD != "NA") %>%
  filter(Managed_Sample.size != "NA") %>%
  filter(Invaded_measure != "NA") %>%
  filter(Invaded_SD != "NA") %>%
  filter(Invaded_Sample.size != "NA") %>%
  filter(Managed_SD > 0) %>%
  filter(Invaded_SD > 0) %>%
  mutate(SDPool = sqrt((((Managed_Sample.size-1)*Managed_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Managed_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Managed_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Managed_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))

##Saving this "real difference" as the object PartB.1.set.testobs, which we will use later
PartB.1.set.testobs <- mean(PartB.1.set$Hedges.G)

##Getting ready for randomisations
reps <- 10000 #Tells R we want 10,000 permutations
PartB.1.set.reps <- numeric(reps)

PartB.1.set %<>%
  select(-SDPool, -N, -Hedges.G) %>%
  gather(., key = "Measure_status", value = "Measure", Managed_measure, Invaded_measure) %>%
  gather(., key = "SD_status", value = "SD", Managed_SD, Invaded_SD) %>%
  gather(., key = "n_status", value = "n", Managed_Sample.size, Invaded_Sample.size) %>%
  mutate(status = ifelse(Measure_status == "Managed_measure", "Managed", "Invaded")) %>%
  select(-Measure_status, -SD_status, -n_status)

for(i in 1:reps) {
  PartB.1.set$statusSHUFFLE <- sample(PartB.1.set$status) #shuffles the invasion status variable
  PartB.1.seta <- PartB.1.set %>%
    filter(statusSHUFFLE == "Managed")
  PartB.1.setb <- PartB.1.set %>%
    filter(statusSHUFFLE == "Invaded")
  colnames(PartB.1.seta)[colnames(PartB.1.seta)=="Measure"] <- "Managed_measure"
  colnames(PartB.1.seta)[colnames(PartB.1.seta)=="SD"] <- "Managed_SD"
  colnames(PartB.1.seta)[colnames(PartB.1.seta)=="n"] <- "Managed_Sample.size"
  PartB.1.seta$Invaded_measure <- PartB.1.setb$Measure
  PartB.1.seta$Invaded_SD <- PartB.1.setb$SD
  PartB.1.seta$Invaded_Sample.size <- PartB.1.setb$n
  PartB.1.seta %<>%
  mutate(SDPool = sqrt((((Managed_Sample.size-1)*Managed_SD^2)+((Invaded_Sample.size-1)*Invaded_SD^2))/(Managed_Sample.size+Invaded_Sample.size-2))) %>%
  mutate(N = Managed_Sample.size + Invaded_Sample.size) %>%
  mutate(Hedges.G = ((Managed_measure-Invaded_measure)/SDPool) * ((N-3)/(N-2.25)) * sqrt((N-2)/N))
  PartB.1.set.reps[i] <- mean(PartB.1.seta$Hedges.G) #Take Hedges.G and save it as a permuted difference in teststat
}

##Generate a histogram of the permuted mean differences, add vertical line of our observed difference
ggplot() + geom_histogram(aes(x = PartB.1.set.reps), binwidth = 0.01) + geom_vline(aes(xintercept = PartB.1.set.testobs), colour = "red", linetype = "dashed", size = 1) + xlab("Permuted Hedges' G") + ylab("Frequency") + theme_cowplot()

##Generate the p-value (two-tailed)
options(digits = 3) #Adjust to see appropriate number of sig-digs
(sum(abs(PartB.1.set.reps) >= abs(PartB.1.set.testobs)) + 1)/(sum(reps) + 1) #p-value

#Step 3: Save permuted results
PartB.1.set.reps <- as.data.frame(PartB.1.set.reps)
fig.PartB.1.bindr <- PartB.1.set.reps %>%
  summarise(avg = mean(PartB.1.set.reps),
            lwr = avg - std.error(PartB.1.set.reps)*1.96,
            upr = avg + std.error(PartB.1.set.reps)*1.96)%>%
  add_column(obs.diff = PartB.1.set.testobs,
             var = Filter.PartA.1,                 
             p = (sum(abs(PartB.1.set.reps) >= abs(PartB.1.set.testobs)) + 1)/(sum(reps) + 1))

```

#Part C: Graphing

```{r Part C}

fig.bindr.A <- rbind(fig.PartA.1.bindr, fig.PartA.2.bindr, fig.PartA.3.bindr) #add or subtract based on how many you were able to permute
fig.bindr.A$comp <- "U-I"

fig.bindr.B <- rbind(fig.PartB.1.bindr) #add based on how many you were able to permute
fig.bindr.B$comp <- "M-I"

fig.bindr <- rbind(fig.bindr.A, fig.bindr.B)

fig.bindr %<>%
  mutate(sig = ifelse(p <= 0.001, "***",
                      ifelse(p <= 0.01, "**",
                             ifelse(p <= 0.05, "*",
                                    ifelse(p < 0.10, "near", "n.s.")))))

fig.bindr

Fig <- ggplot(fig.bindr, aes(x = var, y = avg, group = comp, colour = comp)) + 
  geom_hline(aes(yintercept = 0), colour = "red", linetype = "dashed", size = 1) + 
  geom_point(aes(x = var, y = avg, group = comp), size = 2, position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), size = 0.75, width = 0.2, position = position_dodge(width = 0.5)) +
  geom_point(aes(var, obs.diff, fill = comp), size = 3, colour = "black", shape = 24, position = position_dodge(width = 0.5)) +
  xlab(Filter3) + 
  ylab("Permutated Hedges' G") +
  coord_flip() +
  scale_fill_manual(values = c("steelblue1", "black"), breaks = c("U-I", "M-I"), labels=c("Uninvaded-Invaded", "Managed-Invaded")) +
  scale_color_manual(values = c("steelblue1", "black"), breaks = c("U-I", "M-I"), labels=c("Uninvaded-Invaded", "Managed-Invaded")) +
  scale_y_continuous(expand = c(0, 0), limits = c(-6, 6)) +
  theme(legend.position = c(0.90, 0.050),
        legend.title = element_blank(),
        axis.line.x = element_line(size = 1, colour = "black"), 
        axis.line.y = element_line(size = 1, colour = "black"),
        axis.ticks = element_line(colour = "black", size = 1),
        axis.ticks.length = unit(2.5, "mm"),
        legend.text = element_text(size = 8),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12, face = "bold"))

Fig

```