#install.packages("raster")
#install.packages(c("raster", "ecospat", "Hmisc"))
#install.packages("rgdal")
#install.packages("MuMIn")
#install.packages("maptools")
#install.packages("PresenceAbsence")
#install.packages("Metrics")
library(raster)
library(ecospat)
library(Hmisc)
library(rgdal)
library(MuMIn)
library(maptools)
library(PresenceAbsence)
library(Metrics)

ex=c(66,97,7,38)##extent of the study

##Occurence data####
ag=read.csv("Occurence.27.7.17.csv")
head(ag)
plot(ag$Lon,ag$Lat,pch=0.5,cex=0.5,col="red")

###########################################
####WorldClim data-World Clim variables####
###########################################

####Reading World Clim 2 data#####

path.wc="F:/GIS DATA/World clim 2/wc2.0_30s_bio"
s.wc=list.files(path=paste(path.wc), pattern='tif', full.names=TRUE)
str(s.wc)
s.wc=stack(s.wc)
#plot(s.wc)

####extracting worldclim 2 values#

ag.clm.wc=extract(s.wc, cbind(ag$Lon,ag$Lat),cellnumbers=T)
head(ag.clm.wc)
ag.clm.wc=data.frame(ag.clm.wc)
ag.clm.wc = cbind(ag[,2:5],ag.clm.wc)
head(ag.clm.wc)


ag.clm.wc.sc= scale(ag.clm.wc)



### Cluster analysis:Variable selection for World Clim####

library(Hmisc)
plot( varclus(as.matrix(ag.clm.wc.sc[,6:24]), similarity="spearman",method="average") )
abline(a = 0.49, b = 0, col = "blue")##rho^2 = 0.49/rho=0.7


##Selected variables using correlation threshold of 0.7 (spearman's rho)####
#scale(wc_bio12)+scale(wc_bio14)+scale(wc_bio15)+scale(wc_bio3)+scale(wc_bio6)

##full model with selected variables####
m.f.wc=glm(Occurence~ scale(wc_bio12)+scale(wc_bio14)+scale(wc_bio15)+scale(wc_bio3)+scale(wc_bio6),
           family="binomial",data=ag.clm.wc)
summary(m.f.wc)

####multimodel inference####
library(MuMIn)

dr.wc=dredge(m.f.wc,options(na.action = "na.fail"))
summary(dr.wc)

top.dr.wc=get.models(dr.wc, subset=delta<2)
summary(top.dr.wc)
##
model.sel.table.wc=model.sel(top.dr.wc)

model.avr.wc=model.avg(top.dr.wc)

##

p.avr.wc=predict(model.avr.wc, type="response")#predicted values 

p.avr.wc.r=predict(s.wc,model.avr.wc,ext=ex,type="response")##prediction using the rasters layers for indian subcontinent

plot(p.avr.wc.r, cex.axis=1.5, cex.lab=1.5, main="WorldClim 2")

####

###########################################
####Chelsa data-Chelsa variables##########
###########################################

##Reading Chelsa 2 data#####
path.ch="F:/GIS DATA/Chesla1.2"
s.ch=list.files(path=paste(path.ch), pattern='tif', full.names=TRUE)
str(s.ch)
s.ch=stack(s.ch)
#plot(s.ch)


#extracting Chelsa 2 values#####

ag.clm.ch=extract(s.ch, cbind(ag$Lon,ag$Lat),cellnumbers=T)
head(ag.clm.ch)
ag.clm.ch=data.frame(ag.clm.ch)
ag.clm.ch= cbind(ag[,2:5],ag.clm.ch)
head(ag.clm.ch)


ag.clm.ch.sc=scale(ag.clm.ch)



### Cluster analysis:Variable selection for World Clim####

library(Hmisc)
plot( varclus(as.matrix(ag.clm.ch.sc[,6:24]), similarity="spearman",method="average") )
abline(a = 0.49, b = 0, col = "blue")##rho^2 = 0.49/rho=0.7

###Full model based on selected variables####
m.f.ch=glm(Occurence~scale(CHELSA_bio3)+scale(CHELSA_bio15)+scale(CHELSA_bio12)+scale(CHELSA_bio7)+scale(CHELSA_bio4)
           +scale(CHELSA_bio14)+scale(CHELSA_bio6),
           family="binomial",data=ag.clm.ch)

summary(m.f.ch)


####Multimodel inference:Chelsa####
library(MuMIn)
dr.ch=dredge(m.f.ch,options(na.action = "na.fail"))
summary(dr.ch)

top.dr.ch=get.models(dr.ch, subset=delta<2)

##
model.sel.table=model.sel(top.dr.ch)

d.ch.ch=data.frame(model.sel.table)
#write.csv(d.ch.ch, file="ch.ch.csv")

model.avr.ch=model.avg(top.dr.ch)



p.avr.ch=predict(model.avr.ch, type="response")

p.avr.ch.r=predict(s.ch,model.avr.ch,ext=ex,type="response")
plot(p.avr.ch.r)

###

###################################
###Chelsa data-worldclim variables
####################################

###Full model based on selected variables for worldclim####


m.f.ch.wc=glm(Occurence~scale(CHELSA_bio3)+scale(CHELSA_bio15)+scale(CHELSA_bio12)
              +scale(CHELSA_bio14)+scale(CHELSA_bio6),
              family="binomial",data=ag.clm.ch)

summary(m.f.ch.wc)


####Multimodel inference:Chelsa####
library(MuMIn)
dr.ch.wc=dredge(m.f.ch.wc,options(na.action = "na.fail"))
summary(dr.ch.wc)

top.dr.ch.wc=get.models(dr.ch.wc, subset=delta<2)

##
model.sel.table.ch.wc=model.sel(top.dr.ch.wc)
d.ch.wc=data.frame(model.sel.table.ch.wc)
#write.csv(d.ch.wc, file="ch.wc.csv")

(model.avr.ch.wc=model.avg(top.dr.ch.wc))



p.avr.ch.wc=predict(model.avr.ch.wc, type="response")

p.avr.ch.wc.r=predict(s.ch,model.avr.ch.wc,ext=ex,type="response")
plot(p.avr.ch.wc.r)

##########Model Evaluation####
Evaluate <- function(id, obs, pred){ 
  require(PresenceAbsence)
  require(Metrics)
  comp.prop <- data.frame(id, obs, pred)
  paa <- roc.plot.calculate(comp.prop) 
  paa$tss <- paa$sensitivity + paa$specificity - 1
  paa.grp <- as.vector(paa[which(paa$tss==max(paa$tss)), ])
  pred.spec <- as.numeric(pred >= paa.grp$threshold)
  mse=mean((obs - pred.spec)^2)
  metrics <- cbind(paa.grp,mse)
  metrics
}


##prediction data frame for all the models
pred.df=data.frame(cbind(id=(1:length(ag$Occurence)),obs=ag$Occurence, p.avr.wc,p.avr.ch,p.avr.ch.wc))
head(pred.df)
ev.wc=Evaluate(pred.df$id, pred.df$obs,pred.df$p.avr.wc)
ev.ch=Evaluate(pred.df$id, pred.df$obs,pred.df$p.avr.ch)
ev.ch.wc=Evaluate(pred.df$id, pred.df$obs,pred.df$p.avr.ch.wc)


ev.df= rbind(wc.wc=ev.wc,ch.ch=ev.ch, ch.wc=ev.ch.wc)


####################
##Boyce Index#######
####################

##########################################
###Boyce Index for external validation####
##########################################

#####calculating boyce index for central and eastern himalaya using presence only data###
library(ecospat)
library(rgdal)
library(raster)
##Importing convex hull for presence only data to crop the prediction raster#convex hull was generated in QGIS
shape=readOGR("D:/R/Sep 27_2017_3rd Paper/3rd Paper/Himalaya_Presence_polygon/Himalaya_Presence.shp")
plot(shape)


##WorldClim-WorldClim-Boyce##
pr=p.avr.wc.r#predicted raster from the model to be used for evaluation
masked= mask(pr,shape)
p.r=crop(masked,shape)
plot(p.r)
prs= read.csv("Presence only_Himalaya.csv")##importing presence only data from himalaya
obs1=cbind(prs$lon, prs$lat)
Boyce.p=ecospat.boyce (fit = p.r , obs=obs1, nclass=0,
                       window.w="default", res=100, PEplot = TRUE)

(wc.wc.b=Boyce.p$Spearman.cor)


##Chelsa-Chelsa-Boyce##

rm(Boyce.p)

pr=p.avr.ch.r#predicted raster from the model to be used for evaluation
masked= mask(pr,shape)
p.r=crop(masked,shape)
plot(p.r)
prs= read.csv("Presence only_Himalaya.csv")##importing presence only data from himalaya
obs1=cbind(prs$lon, prs$lat)
Boyce.p=ecospat.boyce(fit = p.r , obs=obs1, nclass=0,
                      window.w="default", res=100, PEplot = TRUE)

(ch.ch.b= Boyce.p$Spearman.cor)

##Chelsa-WorldClim-Boyce##
rm(Boyce.p)

pr=p.avr.ch.wc.r#predicted raster from the model to be used for evaluation
masked= mask(pr,shape)
p.r=crop(masked,shape)
plot(p.r)
prs= read.csv("Presence only_Himalaya.csv")##importing presence only data from himalaya
obs1=cbind(prs$lon, prs$lat)
Boyce.p=ecospat.boyce (fit = p.r , obs=obs1, nclass=0,
                       window.w="default", res=100, PEplot = TRUE)

(ch.wc.b=Boyce.p$Spearman.cor)


rm(Boyce.p)

(result.boyce.ext=rbind(wc.wc.b, ch.ch.b, ch.wc.b))

#########################################
##Boyce index for internal validation####
#########################################

##Importing convex hull for presence only data to crop the prediction raster#convex hull was generated in QGIS
shape.wh=readOGR("D:/R/Sep 27_2017_3rd Paper/3rd Paper/Western Himalaya polygon/ch.survey.wh.shp")##function from raster package
plot(shape.wh)

##WorldClim-WorldClim-Boyce##
rm(Boyce.p)

pr=p.avr.wc.r#predicted raster from the model to be used for evaluation
masked= mask(pr,shape.wh)
p.r=crop(masked,shape.wh)
plot(p.r)
prs.wh=subset(ag,ag$Occurence==1)##subsetting presence only data
obs.wh=cbind(prs.wh$Lon, prs.wh$Lat)
Boyce.p=ecospat.boyce (fit = p.r , obs=obs.wh, nclass=0,
                       window.w="default", res=100, PEplot = TRUE)
(wc.wc.b.wh=Boyce.p$Spearman.cor)

##Chelsa-Chelsa-Boyce##

rm(Boyce.p)

pr.wh=p.avr.ch.r#predicted raster from the model to be used for evaluation
masked= mask(pr.wh,shape.wh)
p.r.wh=crop(masked,shape.wh)
plot(p.r.wh)

prs.wh=subset(ag,ag$Occurence==1)##subsetting presence only data
obs.wh=cbind(prs.wh$Lon, prs.wh$Lat)
Boyce.p=ecospat.boyce(fit = p.r.wh, obs=obs.wh, nclass=0,window.w="default", res=100, PEplot = TRUE)
(ch.ch.b.wh=Boyce.p$Spearman.cor)


##Chelsa-WorldClim##
rm(Boyce.p)

pr=p.avr.ch.wc.r#predicted raster from the model to be used for evaluation
masked= mask(pr,shape.wh)
p.r=crop(masked,shape.wh)
plot(p.r)
prs.wh=subset(ag,ag$Occurence==1)##subsetting presence only data
obs.wh=cbind(prs.wh$Lon, prs.wh$Lat)
Boyce.p=ecospat.boyce (fit = p.r , obs=obs.wh, nclass=0,
                       window.w="default", res=100, PEplot = TRUE)

(ch.wc.b.wh=Boyce.p$Spearman.cor)


result.boyce.int= rbind(wc.wc.b.wh,
                        ch.ch.b.wh,
                        ch.wc.b.wh)

result.boyce.ext


###################
##plotting#######


par(mfrow=c(3,2))

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.wc.r, main="WorldClim 2 data")

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.r, main="Chelsa 1.2 data")

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.wc.r, main="Chelsa 1.2 data - WorldClim 2 variable selection")

##Binary maps

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.wc.r>=0.69, main="WorldClim 2 data")

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.r>=0.46, main="Chelsa 1.2 data")

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.wc.r>=0.54, main="Chelsa 1.2 data - WorldClim 2 variable selection")



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



par(mfrow=c(3,2))

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.wc.r, main="WorldClim 2 data")


par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.wc.r>=0.69, main="WorldClim 2 data")##Binary maps

##

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.r, main="Chelsa 1.2 data")

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.r>=0.46, main="Chelsa 1.2 data")##Binary maps

##

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.wc.r, main="Chelsa 1.2 data - WorldClim 2 variable selection")

par(mar = c(1.9, 1.9, 1.9, 1.9))
plot(p.avr.ch.wc.r>=0.54, main="Chelsa 1.2 data - WorldClim 2 variable selection")##Binary maps

##