#Install necessary packages.

install.packages("biomod2", dependencies=TRUE)

install.packages("raster", dependencies=TRUE)

install.packages("rgdal", dependencies=TRUE)



#Update packages.

#update.packages("biomod2", "raster", "rgdal", repos="http://R-Forge.R-project.org")



#Load packages.

library("biomod2", "raster", "rgdal")



#Create vector of species names.  Comment out those that you do not
#want to model, e.g. due to low sample sizes for evaluation data.

sp.names = c(
"fag.gra",
"ace.sac",
"til.ame",
"tsu.can",
"ulm.ame",
"fra.ame",
"mag.acu",
"cas.den",
"pin.str",
"fra.nig",
"que.vel",
"que.alb",
#"jug.cin",
#"lir.tul",
#"aln.inc",
#"ace.rub",
#"pla.occ",
#"abi.bal",
#"jug.nig",
"car.spp",
"bet.all",
"pru.ser")



#START of 'for' loop.  The 'for' loop repeats the same modeling
#procedures for each species in the above vector.

for(sp.n in sp.names){



#Get a single species.

myRespNameFinal <- sp.n



#Set working directory.

setwd("C:/Users/Steve/Desktop/biomod2/biomod2_results/")



#Create a suffix for your model outputs.  Helps organize outputted
#files.

MyDatasetName <- '_line-description_with_NAVs'



#Load the training data and evaluation data.

DataSpecies_lot_line_section <- read.table(
"C:/Users/Steve/Desktop/biomod2/biomod2_files_to_be_loaded/lot_line_section_with_predictors.csv",
sep=",",
header=TRUE)

lot_line_section_with_predictors <- read.table(
"C:/Users/Steve/Desktop/biomod2/biomod2_files_to_be_loaded/lot_line_section_with_predictors.csv",
sep=",",
header=TRUE)

township_bt_post_with_predictors <- read.table(
"C:/Users/Steve/Desktop/biomod2/biomod2_files_to_be_loaded/township_bt_post_with_predictors.csv",
sep=",",
header=TRUE)

township_line_section_with_predictors <- read.table(
"C:/Users/Steve/Desktop/biomod2/biomod2_files_to_be_loaded/township_line_section_with_predictors.csv",
sep=",",
header=TRUE)



#Create raster stacks, for both training models, and for projecting
#models in geographic space later.

myExpl_with <- stack(
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/aetoverpet07",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/cti",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/drainclass",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/percentclay",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/percentsand",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/ph",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/precipgs",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/slopedegrees",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/solradgs",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/temp01avg",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/tempavggs",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/historicvill",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/lwvill",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/trails")

myExpl_without <- stack(
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/aetoverpet07",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/cti",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/drainclass",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/percentclay",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/percentsand",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/ph",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/precipgs",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/slopedegrees",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/solradgs",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/temp01avg",
"C:/Users/Steve/Desktop/biomod2/biomod2_predictors/tempavggs")



#Select which raster stack of predictor variables you would like to
#use in models.  "myExpl_with" = with NAVs, "myExpl_without" =
#without NAVs.

myExplNameFinal <- myExpl_with



#Create a vector of presences/absences for the species you are
#modeling.

myResp_lot_line_section <- as.numeric(DataSpecies_lot_line_section[,myRespNameFinal])



#Create a vector, containing the geographic coordinates of the
#sample points.

myRespXY_lot_line_section <- DataSpecies_lot_line_section[,c("x_coord","y_coord")]



#Designate vector of presences/absences.

myRespFinal <- myResp_lot_line_section



#Designate vector of geographic coordinates.

myRespXYfinal <- myRespXY_lot_line_section



#Create data frame that contains presences/absences, and all
#predictor variables.  These data frames are used specifically to
#evaluate models.

lot_line_section_with_predictors_species_specific <-
lot_line_section_with_predictors[c(myRespNameFinal, "aetoverpet07", "cti",
"drainclass", "percentclay", "percentsand", "ph", "precipgs",
"slopedegrees", "solradgs", "temp01avg", "tempavggs", "historicvill",
"lwvill", "trails")]

township_bt_post_with_predictors_species_specific <-
township_bt_post_with_predictors[c(myRespNameFinal, "aetoverpet07", "cti",
"drainclass", "percentclay", "percentsand", "ph", "precipgs",
"slopedegrees", "solradgs", "temp01avg", "tempavggs", "historicvill",
"lwvill", "trails")]

township_line_section_with_predictors_species_specific <-
township_line_section_with_predictors[c(myRespNameFinal, "aetoverpet07", "cti",
"drainclass", "percentclay", "percentsand", "ph", "precipgs",
"slopedegrees", "solradgs", "temp01avg", "tempavggs", "historicvill",
"lwvill", "trails")]



#Set important information for models - response variables,
#predictor or 'explanatory' variables, etc.

myBiomodData <- BIOMOD_FormatingData(
resp.var = myRespFinal,
expl.var = myExplNameFinal,
resp.xy = myRespXYfinal,
resp.name = myRespNameFinal,
na.rm = TRUE)



#Set parameters for individual model algorithms.
#See 'biomod2' documentation for more details.

myBiomodOptions <- BIOMOD_ModelingOptions(
GAM = list(interaction.level = 1, algo = 'GAM_mgcv'),
GLM = list(interaction.level = 0, test = 'AIC'),
MARS = list(degree = 2, prune = TRUE),
)


#Run the models.
#See 'biomod2' documentation for more details.

myBiomodModelOut <- BIOMOD_Modeling(
myBiomodData,
models = c('GAM', 'GLM', 'MARS'),
models.options = myBiomodOptions,
NbRunEval=1,
DataSplit=100,
Yweights=NULL,
Prevalence=0.5,
VarImport=3,
models.eval.meth = c('ROC', 'TSS'),
SaveObj = TRUE,
rescal.all.models = TRUE,
do.full.models = TRUE,
modeling.id = myRespNameFinal)



#Project the models into geographic space.

myBiomodProj <- BIOMOD_Projection(
modeling.output = myBiomodModelOut,
new.env = myExplNameFinal,
proj.name = paste(myRespNameFinal, MyDatasetName, '_projection', sep=""),
selected.models = 'all',
binary.meth = NULL,
clamping.mask = FALSE,
output.format = '.img')



#Evaluate the models upon different datasets.

write.table(evaluate(myBiomodModelOut,
data=lot_line_section_with_predictors_species_specific,
stat=c('ROC', 'TSS')),
file = paste(myRespNameFinal, MyDatasetName, "_model_evaluations_upon_lot_line_section_data.csv", sep=""),
sep=",",
append=FALSE,
row.names=TRUE,
col.names=TRUE)

write.table(evaluate(myBiomodModelOut,
data=township_bt_post_with_predictors_species_specific,
stat=c('ROC', 'TSS')),
file = paste(myRespNameFinal, MyDatasetName, "_model_evaluations_upon_township_bt_post_data.csv", sep=""),
sep=",",
append=FALSE,
row.names=TRUE,
col.names=TRUE)

write.table(evaluate(myBiomodModelOut,
data=township_line_section_with_predictors_species_specific,
stat=c('ROC', 'TSS')),
file = paste(myRespNameFinal, MyDatasetName, "_model_evaluations_upon_township_line_section_data.csv", sep=""),
sep=",",
append=FALSE,
row.names=TRUE,
col.names=TRUE)



#Save current workspace before deleting all workspace objects and
#running the next model.  This procedure saves species-specific
#workspaces, so that they can be accessed later.

save.image(file=paste(myRespNameFinal, MyDatasetName, "biomod2.RData"))



#Remove all workspace objects.

rm(list=ls(all=TRUE))



#END of 'for' loop.

}