##### GLMMs of Predictors of LHB Presence in Compartments ######
#### Correlation Matrix ####
# A. Osborne  January 2022-August 2023

# Script for analysis of Astley Habitat Survey dataset 
#with Astley Flight Point Data.
#167 Survey compartments
#The Astley release area Habitat Survey is based on a 10m x 10m grid 
  #of 144 quadrats. The release area is divided into 144 10m compartments
  #centered on each quadrat.
  #additional survey points to cover the larger 2021 (14) and 2022 (7) quadrats.

#Predicting C.tulia Flight Points 
  #with Fixed Effect Habitat Variables 
#and Flight Season (2020, 2021, 2022) as the Random Effect variable.

##### Install and load packages ######
#install.packages("lme4")
#install.packages("MuMIn")
#install.packages("corrplot")
#install.packages("Hmisc")

library(Hmisc)
library(corrplot)
library(lme4)
library(MuMIn)
library(performance)
library(glmm.hp)
#citation("lme4")
#citation("MuMIn")
options(na.action = "na.exclude") 

# Set working directory
setwd("~/Desktop/PC_Folders/LHB/Manuscripts/BC_JIC")
list.files(pattern=".csv")

options(scipen=999)

##########################################
#### GLMMs ####
##########################################

# Load in the Habitat Dispersal data 
H_D <-read.csv("Astley_Habitat_Dispersal_Data.csv")
# check variable structure
str(H_D)
# View variable names for reference
names(H_D)
nrow(H_D)
sum(H_D$C.t_FP)
############

#Remove compartments that arn't part of 2020, 2021 or 2022 polygon
H_D20 <-subset(H_D,H_D$Flt_Season==2020)
H_D20.P <-subset(H_D20, H_D20$Polygon==1)
sum(H_D20.P$C.t_FP);nrow(H_D20.P)# summary stats for data used in GLMM
H_D21 <-subset(H_D,H_D$Flt_Season==2021)
H_D21.P <-subset(H_D21, H_D21$Polygon==1)
sum(H_D21.P$C.t_FP);nrow(H_D21.P)# summary stats for data used in GLMM
H_D22 <-subset(H_D,H_D$Flt_Season==2022)
H_D22.P <-subset(H_D22, H_D22$Polygon==1)
sum(H_D22.P$C.t_FP);nrow(H_D22.P)# summary stats for data used in GLMM

H_D <-rbind(H_D20.P,H_D21.P,H_D22.P)
nrow(H_D)

#Quadrats within polygons by season
with(H_D, tapply(C.t_PA, Flt_Season,length))
# Quadrats with Presence/absence of C.tullia
nrow(H_D20.P)
nrow(subset(H_D20.P, H_D20.P$C.t_PA==1))
nrow(subset(H_D20.P, H_D20.P$C.t_PA==0))

nrow(H_D21.P)
nrow(subset(H_D21.P, H_D21.P$C.t_PA==1))
nrow(subset(H_D21.P, H_D21.P$C.t_PA==0))

nrow(H_D22.P)
nrow(subset(H_D22.P, H_D22.P$C.t_PA==1))
nrow(subset(H_D22.P, H_D22.P$C.t_PA==0))

####Pooled occupancy####
nrow(H_D)
H_D.pres <-subset(H_D, H_D$C.t_PA==1);nrow(H_D.pres)
H_D.abs <-subset(H_D, H_D$C.t_PA==0);nrow(H_D.abs)
#Remove duplicates 
H_D.pres.n <-nrow(H_D.pres[!duplicated(H_D.pres[c("Number")]), ]);H_D.pres.n#Present in at least 1 year
H_D.n <-nrow(H_D[!duplicated(H_D[c("Number")]), ]);H_D.n#number of compartments within the 3year polygon
H_D.n - H_D.pres.n #number absent across all 3 years
(H_D.pres.n/H_D.n)*100 #Percent of compartments present.

# Scale all variables 
# Debate over whether variables should be scaled
# https://statisticsbyjim.com/regression/standardize-variables-regression/ 
#Set class for Variables not to be scaled
H_D$Flt_Season <-as.factor(H_D$Flt_Season)
H_D$C.t_PA <-as.numeric(H_D$C.t_PA)
#H_D$C.t_FP <-as.numeric(H_D$C.t_FP)
#Scale (Not centered, in order to avoid negative values)
H_D$Sphagnum_total <-scale(H_D$Sphagnum_total, center=FALSE)
H_D$Bryophytes_mixed <-scale(H_D$Bryophytes_mixed, center=FALSE)
H_D$E.tetralix <-scale(H_D$E.tetralix, center=FALSE)
H_D$E.vaginatum <-scale(H_D$E.vaginatum, center=FALSE)
H_D$Ev_tussock_count <-scale(H_D$Ev_tussock_count, center=FALSE)
H_D$E.angustifolium <-scale(H_D$E.angustifolium, center=FALSE)
H_D$C.vulgaris <-scale(H_D$C.vulgaris, center=FALSE)
H_D$M.caerulea <-scale(H_D$M.caerulea, center=FALSE)
H_D$ORP <-scale(H_D$ORP, center=FALSE)
H_D$EC <-scale(H_D$EC, center=FALSE)
H_D$CGB_Percent <-scale(H_D$CGB_Percent, center=FALSE)

##### Overview Models of Cottongrass bed habitat use ####
library(lme4)
modCGB_1_0.FP<-glmer(C.t_FP ~ 
                       
                       CGB_1_0+
                       
                       (1|Flt_Season),data = H_D, family=poisson,
                     glmerControl(optimizer = "bobyqa", 
                                  optCtrl = list(maxfun = 100000)))
summary(modCGB_1_0.FP)

modCGB_P.FP<-glmer(C.t_FP ~ 
                     
                     CGB_Percent+
                     
                     (1|Flt_Season),data = H_D, family=poisson,
                   glmerControl(optimizer = "bobyqa", 
                                optCtrl = list(maxfun = 100000)))
summary(modCGB_P.FP)

modCGB_1_0.PA<-glmer(C.t_PA ~ 
                       
                       CGB_1_0+
                       
                       (1|Flt_Season),data = H_D, family=binomial,
                     glmerControl(optimizer = "bobyqa", 
                                  optCtrl = list(maxfun = 100000)))
summary(modCGB_1_0.PA)

modCGB_P.PA<-glmer(C.t_PA ~ 
                     CGB_Percent+
                     (1|Flt_Season),data = H_D, family=binomial,
                   glmerControl(optimizer = "bobyqa", 
                                optCtrl = list(maxfun = 100000)))
summary(modCGB_P.PA) #lowest AICc


####Binomial Models of C.t Presence/Absence####
set.seed(137)

Mod.A<-glmer(C.t_PA ~ 
               E.tetralix:E.vaginatum+
               E.angustifolium+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.A)#AIC 439.5
check_collinearity(Mod.A)

Mod.B<-glmer(C.t_PA ~ 
              E.tetralix+
               E.vaginatum:Ev_tussock_count+
               Bryophytes_mixed+
             (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.B)#AIC 431.8
check_collinearity(Mod.B)

Mod.C<-glmer(C.t_PA ~ 
              E.vaginatum+
              #E.tetralix+
               Sphagnum_total+
               #M.caerulea+
               #C.vulgaris+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.C)#AIC 428.5
check_collinearity(Mod.C)
?check_collinearity
glmm.hp(Mod.C)

Mod.D<-glmer(C.t_PA ~ 
               E.vaginatum:Ev_tussock_count+
               E.tetralix+
               Sphagnum_total+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.D)#AIC 424.9
check_collinearity(Mod.D)

Mod.E<-glmer(C.t_PA ~ 
              E.tetralix:Ev_tussock_count+
              E.vaginatum+
               C.vulgaris+
               #M.caerulea+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.E)#AIC 432.1
check_collinearity(Mod.E)

Mod.F<-glmer(C.t_PA ~
              Ev_tussock_count+
               M.caerulea+
              Sphagnum_total+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.F)#AIC 419.2 #Table 3.
check_collinearity(Mod.F)
r.squaredGLMM(Mod.F)
help(glmm.hp)
Mod.F.hp <-glmm.hp(Mod.F)
Mod.F.hp
plot.glmmhp(Mod.F.hp)

Mod.G<-glmer(C.t_PA ~ 
              Ev_tussock_count+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.G)#AIC 427.5
#check_collinearity(Mod.G)

Mod.H<-glmer(C.t_PA ~ 
              E.tetralix:Ev_tussock_count+
              E.vaginatum+
              E.angustifolium+
              M.caerulea+
              #C.vulgaris+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.H)#AIC 434.5
check_collinearity(Mod.H)

Mod.I<-glmer(C.t_PA ~ 
              E.vaginatum+
              E.angustifolium+
              E.tetralix:Ev_tussock_count+
              C.vulgaris+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.I)#AIC  433.9
check_collinearity(Mod.I)

Mod.J<-glmer(C.t_PA ~ 
              E.tetralix:Ev_tussock_count+
              E.vaginatum+
              Bryophytes_mixed+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.J)#AIC 432.9
check_collinearity(Mod.J)

Mod.K<-glmer(C.t_PA ~ 
              E.tetralix+
              Ev_tussock_count+
              E.vaginatum+
              Sphagnum_total+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.K)#AIC 423.1
check_collinearity(Mod.K)
glmm.hp(Mod.K)

Mod.L<-glmer(C.t_PA ~ 
              E.tetralix+
              Ev_tussock_count+
              Sphagnum_total+
              M.caerulea+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.L)#AIC 421.1 #Table 3.
check_collinearity(Mod.L)
r.squaredGLMM(Mod.L)
Mod.L.hp <-glmm.hp(Mod.L)
Mod.L.hp
plot.glmmhp(Mod.L.hp)

Mod.M<-glmer(C.t_PA ~ 
              E.tetralix+
              E.angustifolium+
               E.vaginatum:Ev_tussock_count+
              C.vulgaris+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.M)#AIC 432.8
check_collinearity(Mod.M)

Mod.N<-glmer(C.t_PA ~ 
               Bryophytes_mixed+
               E.angustifolium+
               #M.caerulea+
               E.tetralix+
               E.vaginatum:Ev_tussock_count+
               (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.N)#AIC 433.6
check_collinearity(Mod.N)

Mod.O<-glmer(C.t_PA ~ 
               E.vaginatum:Ev_tussock_count+
               Bryophytes_mixed+
               (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.O)#AIC 435.3
check_collinearity(Mod.O)
r.squaredGLMM(Mod.O)

Mod.P<-glmer(C.t_PA ~ 
               M.caerulea+
               Ev_tussock_count+
               E.angustifolium+
               Sphagnum_total+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.P)#AIC 421.0 #Table 3.
check_collinearity(Mod.P)
r.squaredGLMM(Mod.P)
Mod.P.hp <-glmm.hp(Mod.P)
Mod.P.hp
plot.glmmhp(Mod.P.hp)
Mod.Q<-glmer(C.t_PA ~ 
              E.tetralix+
              Ev_tussock_count+
              Sphagnum_total+#ORP+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.Q)#AIC 422.0
check_collinearity(Mod.Q)

Mod.R<-glmer(C.t_PA ~ 
              Sphagnum_total+
               E.tetralix:E.vaginatum+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.R)#AIC 424.9
check_collinearity(Mod.R)

Mod.S<-glmer(C.t_PA ~ 
             E.vaginatum+
               E.tetralix+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.S)#AIC 428.5
#check_collinearity(Mod.S)

Mod.T<-glmer(C.t_PA ~ 
              E.tetralix+
              Ev_tussock_count+
              C.vulgaris+
               #EC+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.T)#AIC 427.7
check_collinearity(Mod.T)

Mod.U<-glmer(C.t_PA ~ 
              E.tetralix+
              E.vaginatum+
              Ev_tussock_count+
              C.vulgaris+
              Sphagnum_total+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.U)#AIC 424.1
check_collinearity(Mod.U)
glmm.hp(Mod.U)

Mod.V<-glmer(C.t_PA ~ 
              E.tetralix:E.vaginatum+
              Ev_tussock_count+
              Sphagnum_total+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.V)#AIC 422.0
check_collinearity(Mod.V)

glmm.hp(Mod.V)#returns an error message, glmm.hd requires an interaction term to be 
#generated from the product of E.tetralix * E.vaginatum
H_D$E.t_E.v <-H_D$E.tetralix*H_D$E.vaginatum
#The E.t_E.v term behaves the same as E.tetralix:E.vaginatum when it is run in Mod.Y
Mod.V2<-glmer(C.t_PA ~ 
                E.t_E.v+
                Ev_tussock_count+
                Sphagnum_total+
                (1|Flt_Season),data = H_D, family=binomial,
              glmerControl(optimizer = "bobyqa", 
                           optCtrl = list(maxfun = 100000)))
summary(Mod.V2)
Mod.V2.hp <-glmm.hp(Mod.V2)
Mod.V2.hp
plot.glmmhp(Mod.V2.hp)

Mod.W<-glmer(C.t_PA ~ 
              E.tetralix+
              E.vaginatum+
              C.vulgaris+
              M.caerulea+
               E.angustifolium+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.W)#AIC 433.5
check_collinearity(Mod.W)


Mod.X<-glmer(C.t_PA ~ 
               E.tetralix:E.vaginatum+
              M.caerulea+
              Ev_tussock_count+
              E.angustifolium+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.X)#AIC 432.0
check_collinearity(Mod.X)

Mod.Y<-glmer(C.t_PA ~ 
              #E.t_E.v+
              E.tetralix:E.vaginatum+
              Ev_tussock_count+
              Sphagnum_total+
              M.caerulea+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.Y)#AIC 421.0 #Table 3.
check_collinearity(Mod.Y)
r.squaredGLMM(Mod.Y)
glmm.hp(Mod.Y)#returns an error message, glmm.hd requires an interaction term to be 
#generated from the product of E.tetralix * E.vaginatum
H_D$E.t_E.v <-H_D$E.tetralix*H_D$E.vaginatum
#The E.t_E.v term behaves the same as E.tetralix:E.vaginatum when it is run in Mod.Y
Mod.Y2<-glmer(C.t_PA ~ 
               E.t_E.v+
               Ev_tussock_count+
               Sphagnum_total+
               M.caerulea+
               (1|Flt_Season),data = H_D, family=binomial,
             glmerControl(optimizer = "bobyqa", 
                          optCtrl = list(maxfun = 100000)))
summary(Mod.Y2)
Mod.Y2.hp <-glmm.hp(Mod.Y2)
Mod.Y2.hp
plot.glmmhp(Mod.Y2.hp)

Mod.Z<-glmer(C.t_PA ~ 
              E.tetralix:Ev_tussock_count+
               #E.tetralix:E.vaginatum+
              Bryophytes_mixed+
              E.angustifolium+
              (1|Flt_Season),data = H_D, family=binomial,
            glmerControl(optimizer = "bobyqa", 
                         optCtrl = list(maxfun = 100000)))
summary(Mod.Z )#AIC 441.4
check_collinearity(Mod.Z)

 library(MuMIn)
models.SB<-model.sel(Mod.A,Mod.B,Mod.C,Mod.D,Mod.E,Mod.F,Mod.G,Mod.H, 
                     Mod.I,Mod.J,Mod.K,Mod.L,Mod.M,Mod.N,Mod.O,Mod.P,
                     Mod.Q,Mod.R,Mod.S,Mod.T,Mod.U,Mod.V,Mod.W,Mod.X,Mod.Y,
                     Mod.Z,
                  rank=AICc) 
#Table 3.
models.SB
#Table 4.
importance(models.SB)

#####extract p.values#######
#reformat the coefficients into data.frames from #Table 4.
{M.A <-data.frame(summary(Mod.A)$coefficients[,4])
names(M.A) <-c("p.value")
M.A["Predictor"]<-row.names(M.A)
M.A

M.B <-data.frame(summary(Mod.B)$coefficients[,4])
names(M.B) <-c("p.value")
M.B["Predictor"]<-row.names(M.B)
M.B

M.C <-data.frame(summary(Mod.C)$coefficients[,4])
names(M.C) <-c("p.value")
M.C["Predictor"]<-row.names(M.C)
M.C

M.D <-data.frame(summary(Mod.D)$coefficients[,4])
names(M.D) <-c("p.value")
M.D["Predictor"]<-row.names(M.D)
M.D

M.E <-data.frame(summary(Mod.E)$coefficients[,4])
names(M.E) <-c("p.value")
M.E["Predictor"]<-row.names(M.E)
M.E

M.F <-data.frame(summary(Mod.F)$coefficients[,4])
names(M.F) <-c("p.value")
M.F["Predictor"]<-row.names(M.F)
M.F

M.G <-data.frame(summary(Mod.G)$coefficients[,4])
names(M.G) <-c("p.value")
M.G["Predictor"]<-row.names(M.G)
M.G

M.H <-data.frame(summary(Mod.H)$coefficients[,4])
names(M.H) <-c("p.value")
M.H["Predictor"]<-row.names(M.H)
M.H

M.I <-data.frame(summary(Mod.I)$coefficients[,4])
names(M.I) <-c("p.value")
M.I["Predictor"]<-row.names(M.I)
M.I

M.J <-data.frame(summary(Mod.J)$coefficients[,4])
names(M.J) <-c("p.value")
M.J["Predictor"]<-row.names(M.J)
M.J

M.K <-data.frame(summary(Mod.K)$coefficients[,4])
names(M.K) <-c("p.value")
M.K["Predictor"]<-row.names(M.K)
M.K

M.L <-data.frame(summary(Mod.L)$coefficients[,4])
names(M.L) <-c("p.value")
M.L["Predictor"]<-row.names(M.L)
M.L

M.M <-data.frame(summary(Mod.M)$coefficients[,4])
names(M.M) <-c("p.value")
M.M["Predictor"]<-row.names(M.M)
M.M

M.N <-data.frame(summary(Mod.N)$coefficients[,4])
names(M.N) <-c("p.value")
M.N["Predictor"]<-row.names(M.N)
M.N

M.O <-data.frame(summary(Mod.O)$coefficients[,4])
names(M.O) <-c("p.value")
M.O["Predictor"]<-row.names(M.O)
M.O

M.P <-data.frame(summary(Mod.P)$coefficients[,4])
names(M.P) <-c("p.value")
M.P["Predictor"]<-row.names(M.P)
M.P

M.Q <-data.frame(summary(Mod.Q)$coefficients[,4])
names(M.Q) <-c("p.value")
M.Q["Predictor"]<-row.names(M.Q)
M.Q

M.R <-data.frame(summary(Mod.R)$coefficients[,4])
names(M.R) <-c("p.value")
M.R["Predictor"]<-row.names(M.R)
M.R

M.S <-data.frame(summary(Mod.S)$coefficients[,4])
names(M.S) <-c("p.value")
M.S["Predictor"]<-row.names(M.S)
M.S

M.T <-data.frame(summary(Mod.T)$coefficients[,4])
names(M.T) <-c("p.value")
M.T["Predictor"]<-row.names(M.T)
M.T

M.U <-data.frame(summary(Mod.U)$coefficients[,4])
names(M.U) <-c("p.value")
M.U["Predictor"]<-row.names(M.U)
M.U

M.V <-data.frame(summary(Mod.V)$coefficients[,4])
names(M.V) <-c("p.value")
M.V["Predictor"]<-row.names(M.V)
M.V

M.W <-data.frame(summary(Mod.W)$coefficients[,4])
names(M.W) <-c("p.value")
M.W["Predictor"]<-row.names(M.W)
M.W

M.X <-data.frame(summary(Mod.X)$coefficients[,4])
names(M.X) <-c("p.value")
M.X["Predictor"]<-row.names(M.X)
M.X

M.Y <-data.frame(summary(Mod.Y)$coefficients[,4])
names(M.Y) <-c("p.value")
M.Y["Predictor"]<-row.names(M.Y)
M.Y

M.Z <-data.frame(summary(Mod.Z)$coefficients[,4])
names(M.Z) <-c("p.value")
M.Z["Predictor"]<-row.names(M.Z)
M.Z
}

#Bind all of the lists of coefficients
M26.SB <-rbind(M.A,M.B,M.C,M.D,M.E,M.F,M.G,M.H,M.I,M.J,
            M.K,M.L,M.M,M.N,M.O,M.P,M.Q,M.R,M.S,M.T,
            M.U,M.V,M.W,M.Y,M.X,M.Z)
names(M26.SB)
options(scipen=999)
#Calculate the mean of each Predictor variable (Table 4)
with(M26.SB, tapply(p.value,Predictor,mean))
#Calculate the min of each Predictor variable (Table 4)
with(M26.SB, tapply(p.value,Predictor,min))


##########################################
#### Correlations of numeric variables (Appendix/S2)####
##########################################

# Load in the Habitat Dispersal data 
#Data must be nummeric for corr to work
#H_D contains data which has been scaled and clipped to polygons only
#Need to remove duplicates
H_D.xd <-H_D[!duplicated(H_D[c("Number")]), ]
# First step is to check all numeric variables for correlation
# We can then exclude appropriate variables and avoid multicollinearity
names(H_D.xd)
#Convert factors to numeric
H_D.xd$C.t_PA <-as.numeric(H_D.xd$C.t_PA)
H_D.xd$CGB_1_0 <-as.numeric(H_D.xd$CGB_1_0)
# Create a new dataframe with just the numeric variables
Vars <- (H_D.xd[c(4,17,7:16,21,23)])
head(Vars)
names(Vars)
# Coerce into a data matrix
Var <- as.matrix(Vars)
#Carry out correlation
M <- cor(Var, method = "spearman")
round(M, digits=2)

library(corrplot)
#Corrplot 1.
corrplot(M, method="circle",type="lower", tl.col="black", tl.srt=45)

#calc p-values
cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}


# matrix of the p-value of the correlation
M.p <- cor.mtest(Var)


#Corrplot 2.
corrplot(M, p.mat =M.p, sig.level=0.05, insig = "blank", type="lower", tl.col="black", tl.srt=45)

#Cutom correlogram with correlation coefficients. 
#Corrplot 3.
#?corrplot 
{col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
  cex.before <- par("cex")
  par(cex=0.6)  #sets font size for correlation coef
  
  corrplot(M, method="color", col=col(200),  
           type="lower", 
           p.mat = M.p,sig.level = 0.05, insig = "blank", # p-value 
           # hide correlation coefficient on the principal diagonal
           addgrid="grey",
           addCoef.col = "black",# Add coefficient of correlation,
           cl.cex=1/par("cex"),
           tl.col="black", tl.srt=45, #Text label color and rotation
           # Combine with significance
           
           diag=TRUE
  )
  par(cex=cex.before)
}
#####################

###############################################
####Binomial Models R2 Analysis####
###############################################
set.seed(137)
names(H_D)
##Global Model##

Mod.global<-glmer(C.t_PA ~ 
                    
                    Sphagnum_total+
                    Ev_tussock_count+
                    M.caerulea+
                    E.tetralix+
                    E.tetralix:E.vaginatum+
                    E.vaginatum+
                    E.angustifolium+
                    C.vulgaris+
                    E.vaginatum:Ev_tussock_count+
                    Bryophytes_mixed+
                    E.tetralix:Ev_tussock_count+
                    
                    (1|Flt_Season),data = H_D, family=binomial,
                  glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun = 100000)))
summary(Mod.global)#AIC 430.8
check_collinearity(Mod.global) #The interaction terms are highly collinear
r.squaredGLMM(Mod.global)
library(glmm.hp)
help(glmm.hp)
glmm.hp(Mod.global)#returns an error message, glmm.hd requires interactions term to be 
#generated from the product of E.tetralix * E.vaginatum etc
H_D$E.t_E.v <-H_D$E.tetralix*H_D$E.vaginatum
#The E.t_E.v term behaves the same as E.tetralix:E.vaginatum when it is run in Mod.Y
H_D$E.v_EvT <-H_D$E.vaginatum*H_D$Ev_tussock_count
H_D$E.t_EvT <-H_D$E.tetralix*H_D$Ev_tussock_count

Mod.global2<-glmer(C.t_PA ~ 
                     
                     Sphagnum_total+
                     Ev_tussock_count+
                     M.caerulea+
                     E.tetralix+
                     E.t_E.v+
                     E.vaginatum+
                     E.angustifolium+
                     C.vulgaris+
                     E.v_EvT+
                     Bryophytes_mixed+
                     E.t_EvT+
                     
                     (1|Flt_Season),data = H_D, family=binomial,
                   glmerControl(optimizer = "bobyqa", 
                                optCtrl = list(maxfun = 100000)))
summary(Mod.global2)
check_collinearity(Mod.global2) #The interaction terms are highly collinear

#The last term is dropped from the result matrix, but is passed through into the
#Individual Part R2 calcultation
set.seed(137)
Mod.global2.hp <-glmm.hp(Mod.global2)
Mod.global2.hp
plot.glmmhp(Mod.global2.hp)

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