###################################################################
############### Module 6, Basics of data wrangling ################
###################################################################
### EXCERCISES (Lessons 12 and 13):

# First here is the code to load and cleanup the mammal exclosure vegetation data from the previous excercise file:
library(tidyverse)
mammal_exc_veg <- read_csv("data/mammal_exc_vegetation.csv")

veg_data <- select(mammal_exc_veg, year = YEAR, month=MONTH, site=SITE, block=BLOCK, treatment=TREATMENT,
                          plot=PLOT, subplot=COORDINATE, Abutilon:unknown_F) %>% 
  filter(treatment == "TOTAL" | treatment == "OPEN") %>% 
  select(year:subplot, everything())

### Quick reminder of what the columns mean:
# month = the month when the data were collected
# site = South, central, or north site
# block = 1, 2, or 3 (each site has 3 blocks)
# treatment = OPEN = open plots (no exclosure)
#             MEGA = megaherbivores excluded
#             MESO = meso- and megaherbivores excluded
#             TOTAL = all herbivores excluded;
#             OUT = near to, but outside experimental plots (rarely used).
# plot = unique plot identifier combining SITE, BLOCK, and TREATMENT
# sublot = marker for each of the 49 subplot gridpoint locations in each olot
# bare_ground = presence (B) or absence (0) of bare ground at that m^2 sublot location
# species = presence (B) or absence (0) of that species at that m^2 sublot location


# Now also load the mammal exclosure habitat data "mammal_exc_habitat.csv":
mammal_exc_hab <- read_csv("data/mammal_exc_habitat.csv")

# The columns in this dataset:
# LEVEL = same as site (south, central, north)
# BLOCK = same as block
# TREATMENT = same as treatment
# PLOT = same as plot
# REBAR = same as subplot
# TREE_CANOPY = Describes canopy cover above the subplot/rebar: 
#        Under = fully under canopy; by_under = under canopy, but near edge; 
#        by = at canopy edge; by_open = beside a tree canopy; 
#        open = not under or near a canopy.
# THICKET = When rebar is under or next to a tree, is it part of a thicket? Y/N (N if not near tree)
# HERBACEOUS = Describes density of herbaceous cover around the rebar (seven levels):
#        Bare = no herbaceous cover; few-bare, few, med-few, med, med-high, and high are increasing levels of cover.
# SLOPE = angle of the ground (flat, low, med, steep)
# ASPECT = Cardinal direction that the slope faces
# STONY, SANDY = substrate characters (Y/N)
# WATER_CHANNEL = Proximity of rebar to surface water (11 classes):
#       N = no water association
#       by = beside a surface water runoff channel
#       in = within a surface water runoff channel.
#       For rebar by or in a surface water runoff channel, qualifiers include:
#       Flow = water runs over the area without a distinct channel.
#       Shallow, med, and deep = levels of channel depth.
# species names = which tree species are part of the canopy cover:
#       Y = yes; N = No. Close = rebar in the open, but a canopy is nearby (1–2 m)
############################
# Exercises: 

# 1) Since the vegetation surveys were conducted through repeated surveys, filter the vegetation
# dataframe to only the final year of data:
veg_data_sub <- filter(veg_data, year == 2013)


# 2) Create a case_when function that takes in a species presence value (0, 'B', or NA),
# and converts those to (0, 1, or 0 (NAs are 0s)):

convert_presence <- function(spp_data){
  converted <- case_when(spp_data == 0 ~ 0,
                         spp_data == 'B' ~ 1,
                         is.na(spp_data) ~ 0)
  return(converted)
}
convert_presence(c(0, 'B', NA))

# 3) Use this function to convert all species presence columns:
# Hint: use the across function:
names(veg_data_sub)
veg_data_clean <- mutate(veg_data_sub, across(Abutilon:unknown_F, convert_presence))


# 4) Since there are 49 subplot observations for each plot, summarize those all by plot
# to create a frequency of occurrence of each species in each plot. Hint, you will need to
# use the group_by, summarize, and across functions. Another hint is that the mean of all
# 0 or 1 occurrence values is the same as its frequency:
table(veg_data_clean$plot)

veg_data_summ <- group_by(veg_data_clean, plot) %>% 
  summarize(across(Abutilon:unknown_F, mean))


# 5) Now go back and clean up the habitat data so that the columns in common with the vegetation data
# have the same names (case sensitive) and remove the tree species columns:
hab_data <- select(mammal_exc_hab, site=LEVEL, block=BLOCK, treatment=TREATMENT, plot=PLOT, 
                   subplot=REBAR, tree_canopy=TREE_CANOPY, thicket=THICKET, herbaceous=HERBACEOUS,
                   slope=SLOPE, aspect=ASPECT, water_channel=WATER_CHANNEL, stony=STONY,
                   sandy=SANDY)


# 6) Filter the habitat data to only those treatments that were either "OPEN" or "TOTAL" exclosure.
# Also remove the columns thicket through water_channel:
hab_data <- filter(hab_data, treatment == "OPEN" | treatment == "TOTAL") %>% 
  select(-c(thicket:water_channel))


# 7) Your next challenge is to summarize the tree_canopy, stony, and sandy columns by plot.
# Can you figure out a clever way to turn those variables into numeric summaries?
# Be sure to check out the solutions video on this exercise file to see an example
# solution if you are stumped on this one. Hint: keep the site, block, and treatment
# columns by also including them in the group_by function.

# always good to first look at tabular summaries to better understand these variables:
table(hab_data$tree_canopy) # looks like we could just simplify to open canopy versus 
# some canopy present (by, under, etc.)
# Stony and sandy are self explanatory yes' and nos

hab_plot_data <- group_by(hab_data, plot, site, block, treatment) %>% 
  summarize(prop_near_trees = mean(as.numeric(tree_canopy != "open")),
            stony_ground = mean(as.numeric(stony == "Y")),
            sandy_ground = mean(as.numeric(sandy == "Y"))) %>% 
  ungroup()

# Hmm, but why the NA? Let's subset the data to see whats up with
# the C2OPEN plot:
temp <- filter(hab_data, plot=="C2OPEN")

# ahh, it looks like there are some missing values there, so we can add in na.rm
# to our mean functions to ignore that one missing subplot

hab_plot_data <- group_by(hab_data, plot, site, block, treatment) %>% 
  summarize(prop_near_trees = mean(as.numeric(tree_canopy != "open"), na.rm=T),
            stony_ground = mean(as.numeric(stony == "Y"), na.rm=T),
            sandy_ground = mean(as.numeric(sandy == "Y"), na.rm=T)) %>% 
  ungroup()

# 8) Join the vegetation summary data from question 4 to the habitat summary data from the previous question:
all_data <- left_join(hab_plot_data, veg_data_summ)
all_data <- left_join(hab_plot_data, veg_data_summ, by="plot")

# 9) Create a species X plot matrix using the vegetation frequency data from question 4:
veg_data_summ
plot_names <- veg_data_summ$plot

spp_plot_table <- select(veg_data_summ, -plot)

spp_plot_mat <- as.matrix(spp_plot_table)
row.names(spp_plot_mat) <- plot_names

spp_plot_mat


# 10) Some species were not found in any plot for some reason, so remove all the species
# columns where the species had zero occurrences in the matrix. Hint: since it's a matrix
# you'll have to rely on some old-school indexing and the apply function.

# lets get the total sum of frequencies of each species:
spp_sums <- apply(spp_plot_mat, MARGIN = 2, FUN = sum)
# where the total sum is zero that means it does not occur at all, so remove that column, or keep only those
# that are greater than zero:
spp_plot_mat_fin <- spp_plot_mat[,spp_sums > 0]
ncol(spp_plot_mat_fin)
ncol(spp_plot_mat)

# 11) Which species has the highest frequency in any one plot?
sort(apply(spp_plot_mat_fin, 2, max))
# Cynodon_plechtostachyus

# 12) Which plot has the greatest number of species?
count_species <- function(x) sum(x > 0)
sort(apply(spp_plot_mat_fin, 1, count_species))
# S3OPEN

# 13) Create a scatterplot showing how the species richness in each plot is related to the 
# relative proportion of overstory trees in the habitat data, and a boxplot showing the 
# relationship between species richness and site.

# First calculate species richness at in each plot:
spp_rich <- apply(spp_plot_mat_fin, 1, FUN = function(x) sum(x > 0))

# create a tibble dataframe with these values and their name labels:
spp_rich_data <- tibble(spp_rich, plot = names(spp_rich))

# THen add this on to the habitat data:
data_fin <- left_join(hab_plot_data, spp_rich_data, by="plot")

# Then plot the figures:
plot(spp_rich ~ prop_near_trees, data=data_fin)
boxplot(spp_rich ~ site, data=data_fin)

# Make them a bit nicer:
plot(spp_rich ~ prop_near_trees, data=data_fin, pch=16, cex=1.5, ylab="Species Richness", xlab="Proportion of subplots\nwith overstory trees")
boxplot(spp_rich ~ site, data=data_fin, ylab="Species Richness", xlab = "Site", col=3)


# That's it!

