###########################################################################################################
###   BIO252 Coursework -- Data Analysis Assignment -- R codes	  
###
###   Use this format to organise the R code to include in the Appendix -- obviously modify as needed for 
###   your submission
###
###   Title: Atlantic forest biodiversity - effects of altitude and bioregion                           
###
###													  
###   Author: Name Surname (ID)										  
###													  
###   Date: dd/mm/yyyy					# --> insert the date you hand in your report
###													  
###   R codes associated with ".....pdf" # --> (Name of the file of your report)				  
###													  
###   R version 3.4.2 (2017-09-28)		# --> change to the R version you used (type 'citation()' into R)
###########################################################################################################

###########################################################################################################
###   Content:											  	  
###													  
### 	1. Data used		# --> (provide links and brief explanation)								          
###												                                                                                               
###	2. Data management & preparation		# --> (use the codes I provide below - select only those for the dataset you will be using, of course!)						  
###													    												  
### 	3. Graphical analysis			  	 						  
###												                                                                                               
###	4. Statistical analysis										  
###													  
###	5. R codes for figures in the report		## --> (i.e. the codes used to generate the specific figures included in your report)						  
###													 
###########################################################################################################


###########################################################################################################

#################################################
# 1. Data used
#################################################

# Data obtained from: ATLANTIC: Data Papers from a biodiversity hotspot
# url: https://esajournals.onlinelibrary.wiley.com/doi/toc/10.1002/(ISSN)1939-9170.AtlanticPapers

# Students can choose among 2 datasets:

##########################
# Atlantic butterfly data
# dos Santos et al. (2018) Atlantic butterflies: a data set of fruit-feeding butterfly communities from the Atlantic forests. Ecology (in press)
# DOI: https://doi.org/10.1002/ecy.2507

# Data files: ATLANTIC_BUTTERFLIES_sites.csv & ATLANTIC_BUTTERFLIES_species.csv
# Metadata description: ecy2507-sup-0002-metadatas1.pdf
# Date last accessed: 26/11/2018
##########################

##########################
# Atlantic bird data
# Hasui et al. (2018) ATLANTIC BIRDS: a data set of bird species from the Brazilian Atlantic Forest. Ecology 99(2): 497-497
# DOI: https://doi.org/10.1002/ecy.2119

# Data file: ATLANTIC_BIRDS_quantitative.csv
# Metadata description: ecy2119-sup-0001-metadatas1.pdf
# Date last accessed: 26/11/2018
##########################




#################################################
# 2. Data management and preparation
#################################################

##########################################################################
# Butterfly data
# read in the data (or use the graphical interface in R Studio, as usual)

  species <- read.csv("ATLANTIC_BUTTERFLIES_species.csv", sep = ";")
  sites <- read.csv("ATLANTIC_BUTTERFLIES_sites.csv", sep = ";")
  
# calculate species richness for each site, then combine with altitude and bioregion information
# combine this information into a data frame, ready for the statistical analysis  
# using table() I count the number of records for each site
# using data.frame() I transform it into a dataframe with the columns required for the analyses   
  
  SpeciesRichness <- table(species$sites_ID)
  ButterflyDiversity <- data.frame(SiteID = factor(names(SpeciesRichness)), SpeciesRichness = as.numeric(SpeciesRichness))

# given that the 'siteID' variable is in the same order in the 'sites' file and the new dataframes, 
# we can easily add information from the former to the latter
# to verify this, check that the order is the same, using the all.equal() function
  
  all.equal(sites$sites_ID, ButterflyDiversity$SiteID)   # OK - same order

  ButterflyDiversity$Olsong200r <- sites$Olsong200r
  ButterflyDiversity$Altitude1km <- sites$Altitude1km
  ButterflyDiversity$Latitude <- sites$Latitude
  ButterflyDiversity$Longitude <- sites$Longitude

# transform column Olsong200r into a factor
  
  ButterflyDiversity$Olsong200r <- factor(ButterflyDiversity$Olsong200r)
  
  
# data frame now ready for the statistical analysis
# I suggest to use a log-transformation of the response variable   
# if the residual plots highlight some issues, discuss them, but do not use other transformations
##################################################################################################  
  
# to further investigate if the pattern differs between different taxonomic groups (Subfamilies), 
# we calculate species richness separately for each of the four subfamilies  
# I provide here the R code for you to do this, using the aggregate() function
# it is a useful function as it provides a dataframe as output, so no need for further changes
# however, see below why we have to add 'drop = FALSE'  
# read carefully until the end, as I do various operations below, then give a precise indication on
# which dataset to use  
  
  
  SubFamilyRichness <- aggregate(Species ~ Subfamily + sites_ID, data = species, FUN = length)
  # this does not consider sites where no species of a given subfamily was found ('zero records')
  # to change this, just add 'drop = FALSE'

  foo2 <- aggregate(Species ~ Subfamily + sites_ID, data = species, FUN = length, drop = FALSE)
  # this sets those cases to 'NA', see:
  
  summary(foo2$Species)
  # hence to transform to 'SpeciesRichness = 0', we need one more line of code
  
  foo2$Species[is.na(foo2$Species)] <- 0
  summary(foo2$Species)
  summary(foo2$Species)
  
  # this should include hence 19 values equal to zero - let's check:
  length(which(foo2$Species == 0))
  
# now, to combine with the spatial information from the 'sites' dataframe, we need a different approach
# the 'merge' function is very useful for this - it takes advantage of the fact that the 'sites_ID' column 
# is identical for both data frames and hence can be linked to uniquely combine information from the 
# 2 data frames, to create a new, merged data frame
  
  foo3 <- merge(foo2, sites, by = "sites_ID")

  # for the analyses for the report, however, I suggest you to focus on the non-zero data.
  # this is why I have called these latter dataframes with a name such as 'foo'
  # I name temporary objects, which I can remove, with names like 'foo', for simplicity
  # it avoids confusion with the 'ButterflyDiversity' dataframe, which you should use for the analyses    
  # let us hence remove the objects not needed, and use the below dataframe for analyses
  
  rm(foo2, foo3)  
  SubFamilyRichness2 <- merge(SubFamilyRichness, sites, by = "sites_ID")
  
  SubFamilyRichness2$Subfamily <- factor(SubFamilyRichness2$Subfamily)
  
  
# Thus, using the SubFamilyRichness2 dataframe, run 4 separate models/analyses, 
# for each subfamily, and investigate if the same results are found, and discuss the results accordingly
# IMPORTANT: you must read up some key information about these subfamilies and use it in your discussion 
#######################################################################################################################

  
  
##########################################################################  
##########################################################################
# Bird data
# read in the data (or use the graphical interface in R Studio, as usual)
  
  BirdSpecies <- read.csv("ATLANTIC_BIRDS_quantitative.csv")
  # this file has already both species and spatial information, hence we can generate the dataframe for 
  # the analyses with simply one operation of data aggregation 
  # (count of number of species detected per site sampled)
  
  BirdRichness <- aggregate(Species ~ Longitude_x + Latitude_y + OlsonG200r + Altitude, 
                                 data = BirdSpecies, FUN = length)
  
  BirdRichness$OlsonG200r <- factor(BirdRichness$OlsonG200r)
  
  summary(BirdRichness)
  # highlights 2 records without bioregion information
  # column OlsonG22r --> exclude those 2 records
  
  BirdRichness2 <- subset(BirdRichness, subset = OlsonG200r != "0")
  BirdRichness2 <- droplevels(BirdRichness2)
  summary(BirdRichness2)   # there appears to be a record with a very high species diversity value
  # inspect this further:
  hist(BirdRichness2$Species)  # suggests it may be best to exclude this outlying value above 400
  
  BirdRichness3 <- subset(BirdRichness2, subset = Species <= 400)
  # check again the distribution of diversity/richness values:
  hist(BirdRichness3$Species)   # this looks more reasonable now
  
# this dataset is now ready for the statistical analysis
# I suggest to use log tranformation for the response variable, as we deal with count data
# if the residual plots highlight some issues, discuss them, but do not use other transformations
#################################################################################################
  
  # to further investigate if the pattern differs between different taxonomic groups (Subfamilies), 
  # we calculate species richness separately for each of the four subfamilies  
  
  BirdRichnessMethod <- aggregate(Species ~ Longitude_x + Latitude_y + Methods + OlsonG200r + Altitude, 
                            data = BirdSpecies, FUN = length)
  
  # again, exclude the outlier
  BirdRichnessMethod2 <- subset(BirdRichnessMethod, subset = Species <= 400)

# convert categorical columns into a factor
  
  BirdRichnessMethod2$Methods <- factor(BirdRichnessMethod2$Methods)
  BirdRichnessMethod2$OlsonG200r <- factor(BirdRichnessMethod2$OlsonG200r)

# take out records without a value for the ecoregion (N = 2 records)
  
  BirdRichnessMethod3 <- subset(BirdRichnessMethod2, subset = OlsonG200r != "0")
  BirdRichnessMethod3 <- droplevels(BirdRichnessMethod3)
  summary(BirdRichnessMethod3)
  
# repeat now 3 separate analyses, fitting the same model but subsetting by method 
# (Point count vs. Mist net vs. Line transect)   
# compare and discuss the results to the analysis not separating by bird sampling method   
# IMPORTANT: you must read up a few key publications about these 3 bird monitoring methods  
#########################################################################################  

  
  
##################################################################################################  
##################################################################################################
# Here below you can continue to use the same structure, to insert your codes
# which you will use for all your analyses
# of course you can modify the structure/format overall
# importantly, it has to be consistent, well organised, easy to read and understand  
##################################################################################################
##################################################################################################  

  
    
#################################################
# 3. Graphical analysis to inspect the data
#################################################
  
# ...
  
  
#################################################
# 4. Statistical analysis
#################################################
  
# ...
  
#################################################
# 5. R codes for figures in the report
#################################################
  
# ...

  
#################################################
# END
#################################################