#######################################
##### Load the required libraries #####
#######################################

library(reshape2)
library(lme4)
library(vegan)
library(visreg)

##################################
##### step 1: Matrix reshape #####
##################################

# we need to reshape the original matrix downloaded from Lifewatch data repository (as csv file). We want to investigate the site vulnerability, distributed across different Eunis habitat, to different taxonomic groups. So we need to aggregate the data at site level by taxon name and Habita Eunis name

# first load raw data, in this example we'll used the marine dataset from Corriero et al. (2016)
#freshwater<-read.csv(file="Dataset_Biodiversity_AlienSpecies_Freshwaters_2015.csv",sep=";")
dataset<-read.csv(file="Dataset_Biodiversity_Marine_Lifewatch_2015.csv",sep=";")

# we are interested in investigating the Eunis habitat at level-1. We create a new variable.
dataset $EunisL1 <- as.factor(substr(dataset $eunishabitatstypecode, start = 1, stop = 2)) 

# ONLY FOR FRESHWATER Additional step.Eunis C3 and J5 were poorly represented. We will remove them from the dataset. Take, remove or change this step acording to the dataset
#dataset <-droplevels(subset(dataset,EunisL1=="c1"| EunisL1=="c2"))

# Subsetting
# we subset the dataset in alien and native
alien<-subset(dataset,alien=="1")
native<-subset(dataset,alien!="1")

# Reshaping and richness calculation 
## reshape alien and calculate the alien richness
alien.melt<-melt(alien)
alien_table<-dcast(alien.melt,locality+ EunisL1 + eunisspeciesgroups ~ scientificname)
alien_table[is.na(alien_table)]<-0
alien_richness<-specnumber(alien_table[sapply(alien_table, class)!="factor"])
alien_richness<-cbind(alien_table[sapply(alien_table, class)=="factor"],alien_richness)
## reshape native and calculate the native richness 
native.melt<-melt(native)
native_table<-dcast(native.melt, locality+ EunisL1 + eunisspeciesgroups ~ scientificname)
native_table[is.na(native_table)]<-0
native_richness<-specnumber(native_table[sapply(native_table, class)!="factor"])
native_richness<-cbind(native_table[sapply(native_table, class)=="factor"],native_richness)
##merge native+aliene
new_table<-merge(native_richness,alien_richness,all.x=T,all.y=T)
new_table[is.na(new_table)]<-0

# remove unnecessary files 
rm("alien_richness")
rm("native_richness")




#########################################################################################
#####  Step 2: Generalized Linear Mixed Model (GLMM) fitting usign the lme4 package #####
#########################################################################################


# now we are ready to fit our model. We will use a generalized linear mixed models in order to take into account the structure of our new dataset. Taxonomic group and locality are not the focus of our investigation but largely influence our sampling. We will include these two factor in the random effect. 

# First fit full model (a negative bionomial family is assumet for richness data)
gfit_Eu_Ri <- glmer.nb(alien_richness ~native_richness+ EunisL1 +(1| eunisspeciesgroups)+(1|locality), data= new_table)

# automatically calculate best model according to AIC
library(MuMIn)
options(na.action = "na.fail")
ms1<-dredge(gfit_Eu_Ri)
ms1; # the full model has the highest AICc support 

# fit the best model accordin to AIC
mod.fit<-glmer.nb(as.formula(getCall(ms1,1)), data = new_table)

# results
summary(mod.fit) #table.
visreg(mod.fit,trans=exp,nn=101,alpha=1,rug=F,partial=T) #graph



######################################
##### Step 3: RAREFACTION CURVES #####
######################################

library(vegan)

asfw<-droplevels(subset(freshwater, alien=="1"))
tabella_aspecie<-table(asfw $EunisL1, asfw $scientificname)
rarecurve(tabella_aspecie)