####Chi_sq test of Rest Points vs Perching Substrate Cover (Table 2) ######
#### Basic Stats on Habitat Data (Table 1)####
## A. Osborne  January 2023 - August 2023 ##

#Display decimals (for p-values)
options(scipen=999)

##Load 2020-2023 Flight Point Data##
setwd("~/Desktop/PC_Folders/LHB/Manuscripts/BC_JIC")
list.files()
Y1_3 <-read.csv("2020-23 Flight Point Data.csv")
names(Y1_3)
nrow(Y1_3)
#Remove Outliers (Dispersal out of the intensively surveyed reintroduction area)
Y1_3.XO <-subset(Y1_3, Dispersal!="Outlier")
nrow(Y1_3.XO)

###################################
#Extract Flight/ Rest Point Statistics
FP.Tot <-nrow(Y1_3.XO); FP.Tot
#Perching on each plant species
FP.unknown <-nrow(subset(Y1_3.XO,Behaviour==""));FP.unknown
FP.Flight <-nrow(subset(Y1_3.XO,Behaviour=="Flight"));FP.Flight
FP.Pair <-nrow(subset(Y1_3.XO,Behaviour=="Pair"));FP.Pair
FP.Mc <-nrow(subset(Y1_3.XO,Behaviour=="Rest_M.caerulea"));FP.Mc
FP.Cv <-nrow(subset(Y1_3.XO,Behaviour=="Rest_C.vulgaris"));FP.Cv
FP.Ev <-nrow(subset(Y1_3.XO,Behaviour=="Rest_E.vaginatum"));FP.Ev
FP.Et <-nrow(subset(Y1_3.XO,Behaviour=="Rest_E.tetralix")); FP.Et
FP.N <-nrow(subset(Y1_3.XO,Behaviour=="Nectaring"));FP.N
FP.Ea <-nrow(subset(Y1_3.XO,Behaviour=="Rest_E.angustifolium"));FP.Ea

## Basic Flight Point Calcs
#Perching on E. tetralix shrub or inflorescence
#Total is FP.XLH (cross-leaved heath)
FP.XLH <-FP.Et + FP.N ; FP.XLH
FP.N/FP.XLH # Fraction nectaring on inflorescence
##Significance test for nectaring, ie Perching on Inflorescence more often than expected##
Inf.EXP_ON = 15 ; Inf.EXP_ON #Expected to be on habitat plant species
Inf.EXP_OFF = 15 ; Inf.EXP_OFF   # Expected to be off habitat plant species
Inf.OBS_ON = FP.N ; Inf.OBS_ON #Observed on habitat plant species
Inf.OBS_OFF = FP.Et ; Inf.OBS_OFF #Observed off habitat plant species
#set up chi table
Inf.Mc.chi<-as.table(rbind(c(Inf.OBS_ON,Inf.EXP_ON),c(Inf.OBS_OFF,Inf.EXP_OFF)))
dimnames(Inf.Mc.chi) <- list(Habitat = c("Inf.ON","Inf.OFF"), Presence = c("Observed", "Expected"))   #labelling Chisq axis
#chi table
Inf.Mc.chi
#perform chi-squared test
Inf.Mc.chi.scores<-chisq.test(Inf.Mc.chi)
#test result
Inf.Mc.chi.scores
Inf.OBS_ON/(Inf.OBS_ON+Inf.OBS_OFF) #proportion observed ON



FP.Rests <-FP.Mc+FP.Cv+FP.Ea+FP.Ev+FP.XLH ; FP.Rests
FP.Rests+FP.Flight+FP.Pair+(FP.unknown+1) #Breakdown of flights 1 pair not followed ie unknown
FP.Tot
FP.XLH; FP.Rests; FP.XLH/FP.Rests

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

####Habitat Survey (2021 initial 144 ARG quadrats, 
  #plus extra quadrats)  combined with FP data (Table 2) ####
  #Testing the null hypothesis that
    #the probability of Rest Flight Points occurring on each specific 
    #habitat plant  species is not dependent on the 
    # % cover of habitat plant species

##########################################
#### Load the Habitat Data ####
##########################################

# Load in the Habitat Dispersal data 
H_D <-read.csv("Astley_Habitat_Dispersal_Data.csv")
names(H_D)
nrow(H_D)

#Remove quadrats/compartments that aren't included in any polygon
H_D.P <-subset(H_D,H_D$Polygon==1)
nrow(H_D.P)
#Remove duplicates that are part of more than one of the 2020, 2021, 2022 polygon
H_D.P <-H_D.P[!duplicated(H_D.P[c("Number")]), ]
#Survey quadrats/compartments included
nrow(H_D.P)
head(H_D.P)
H_D.P$Survey_quadrat

#####################
#Extract %cover data and calculate mean %cover
HD.Mc <-mean(H_D.P$M.caerulea); HD.Mc
HD.Cv <-mean(H_D.P$C.vulgaris); HD.Cv
HD.Ev <-mean(H_D.P$E.vaginatum); HD.Ev
HD.Et <-mean(H_D.P$E.tetralix); HD.Et
HD.Ea <-mean(H_D.P$E.angustifolium); HD.Ea
#Sum of %cover of all plant species used for perching
HD.sum <-sum(HD.Mc,HD.Cv,HD.Ev,HD.Et,HD.Ea);HD.sum #99.49 

#### Table 2;  chi-squared tests for significance of perching plant species ##
#Stats based on plant species % cover  
  #in quadrats within C.tullia dispersal polygons.

# M.caerulea chi-squared
#Set expected values according to %cover 
FP.Mc#Resting Events on M.caerulea
Mc.EXP_ON = FP.Rests*(HD.Mc/HD.sum) ; Mc.EXP_ON #Expected to be on habitat plant species
Mc.EXP_OFF = FP.Rests*(1-(HD.Mc/HD.sum)) ; Mc.EXP_OFF   # Expected to be off habitat plant species
Mc.OBS_ON = FP.Mc ; Mc.OBS_ON #Observed on habitat plant species
Mc.OBS_OFF = FP.Rests-FP.Mc ; Mc.OBS_OFF #Observed off habitat plant species
Mc.OBS_ON+Mc.OBS_OFF
Mc.EXP_ON+Mc.EXP_OFF
#set up chi table
Mc.chi<-as.table(rbind(c(Mc.OBS_ON,Mc.EXP_ON),c(Mc.OBS_OFF,Mc.EXP_OFF)))
dimnames(Mc.chi) <- list(Habitat = c("Mc.ON","Mc.OFF"), Presence = c("Observed", "Expected"))   #labelling Chisq axis
#chi table
Mc.chi
#perform chi-squared test
Mc.chi.scores<-chisq.test(Mc.chi)
#test result
Mc.chi.scores #p<0.001
Mc.OBS_ON/(Mc.OBS_ON+Mc.OBS_OFF) #proportion observed ON
#Mc.OBS_OFF/(Mc.OBS_ON+Mc.OBS_OFF) #proportion observed OFF
HD.Mc/HD.sum #proportion expected ON

# E.tetralix chi-squared
#Set expected values according to %cover 
FP.XLH#Resting Events Total on E.tetralix
Xlh.EXP_ON = FP.Rests*(HD.Et/HD.sum) ; Xlh.EXP_ON #Expected to be on habitat plant species
Xlh.EXP_OFF = FP.Rests*(1-(HD.Et/HD.sum)) ; Xlh.EXP_OFF   # Expected to be off habitat plant species
Xlh.OBS_ON = FP.XLH ; Xlh.OBS_ON #Observed on habitat plant species
Xlh.OBS_OFF = FP.Rests-FP.XLH ; Xlh.OBS_OFF #Observed off habitat plant species
Xlh.OBS_ON+Xlh.OBS_OFF
Xlh.EXP_ON+Xlh.EXP_OFF
#set up chi table
Xlh.chi<-as.table(rbind(c(Xlh.OBS_ON,Xlh.EXP_ON),c(Xlh.OBS_OFF,Xlh.EXP_OFF)))
dimnames(Xlh.chi) <- list(Habitat = c("Xlh.ON","Xlh.OFF"), Presence = c("Observed", "Expected"))   #labelling Chisq axis
#chi table
Xlh.chi
#perform chi-squared test
Xlh.chi.scores<-chisq.test(Xlh.chi)
#test result
Xlh.chi.scores #p<0.0001
Xlh.OBS_ON/(Xlh.OBS_ON+Xlh.OBS_OFF) #proportion observed ON
#Xlh.OBS_OFF/(Xlh.OBS_ON+Xlh.OBS_OFF)#proportion observed OFF
HD.Et/HD.sum #proportion expected ON

# E.vaginatum chi-squared  
#Set expected values according to %cover 
FP.Ev#Resting Events on E.vaginatum 
Ev.EXP_ON = FP.Rests*(HD.Ev/HD.sum) ; Ev.EXP_ON #Expected to be on habitat plant species
Ev.EXP_OFF = FP.Rests*(1-(HD.Ev/HD.sum)) ; Ev.EXP_OFF   # Expected to be off habitat plant species
Ev.OBS_ON = FP.Ev ; Ev.OBS_ON #Observed on habitat plant species
Ev.OBS_OFF = FP.Rests-FP.Ev ; Ev.OBS_OFF #Observed off habitat plant species
Ev.OBS_ON+Ev.OBS_OFF
Ev.EXP_ON+Ev.EXP_OFF
#set up chi table
Ev.chi<-as.table(rbind(c(Ev.OBS_ON,Ev.EXP_ON),c(Ev.OBS_OFF,Ev.EXP_OFF)))
dimnames(Ev.chi) <- list(Habitat = c("Ev.ON","Ev.OFF"), Presence = c("Observed", "Expected"))   #labelling Chisq axis
#chi table
Ev.chi
#perform chi-squared test
Ev.chi.scores<-chisq.test(Ev.chi)
#test result
Ev.chi.scores #p=0.029
Ev.OBS_ON/(Ev.OBS_ON+Ev.OBS_OFF) #proportion observed ON
#Ev.OBS_OFF/(Ev.OBS_ON+Ev.OBS_OFF) #proportion observed OFF
HD.Ev/HD.sum #proportion expected ON

# C.vulgaris chi-squared
FP.Cv#Resting Events on C. vulgaris 
#Set expected values according to %cover 
Cv.EXP_ON = FP.Rests*(HD.Cv/HD.sum) ; Cv.EXP_ON #Expected to be on habitat plant species
Cv.EXP_OFF = FP.Rests*(1-(HD.Cv/HD.sum)) ; Cv.EXP_OFF   # Expected to be off habitat plant species
Cv.OBS_ON = FP.Cv ; Cv.OBS_ON #Observed on habitat plant species
Cv.OBS_OFF = FP.Rests-FP.Cv ; Cv.OBS_OFF #Observed off habitat plant species
Cv.OBS_ON+Cv.OBS_OFF
Cv.EXP_ON+Cv.EXP_OFF
#set up chi table
Cv.chi<-as.table(rbind(c(Cv.OBS_ON,Cv.EXP_ON),c(Cv.OBS_OFF,Cv.EXP_OFF)))
dimnames(Cv.chi) <- list(Habitat = c("Cv.ON","Cv.OFF"), Presence = c("Observed", "Expected"))   #labelling Chisq axis
#chi table
Cv.chi
#perform chi-squared test
Cv.chi.scores<-chisq.test(Cv.chi)
#test result
Cv.chi.scores #p<0.250
Cv.OBS_ON/(Cv.OBS_ON+Cv.OBS_OFF) #proportion observed ON
#Cv.OBS_OFF/(Cv.OBS_ON+Cv.OBS_OFF) #proportion observed OFF
HD.Cv/HD.sum #proportion expected ON

# E.angustifolium chi-squared
FP.Ea#Resting Events on E.angustifolium
#Set expected values according to %cover 
Ea.EXP_ON = FP.Rests*(HD.Ea/HD.sum) ; Ea.EXP_ON #Expected to be on habitat plant species
Ea.EXP_OFF = FP.Rests*(1-(HD.Ea/HD.sum)) ; Ea.EXP_OFF   # Expected to be off habitat plant species
Ea.OBS_ON = FP.Ea ; Ea.OBS_ON #Observed on habitat plant species
Ea.OBS_OFF = FP.Rests-FP.Ea ; Ea.OBS_OFF #Observed off habitat plant species
Ea.OBS_ON+Ea.OBS_OFF
Ea.EXP_ON+Ea.EXP_OFF
#set up chi table
Ea.chi<-as.table(rbind(c(Ea.OBS_ON,Ea.EXP_ON),c(Ea.OBS_OFF,Ea.EXP_OFF)))
dimnames(Ea.chi) <- list(Habitat = c("Ea.ON","Ea.OFF"), Presence = c("Observed", "Expected"))   #labelling Chisq axis
#chi table
Ea.chi
#perform chi-squared test
Ea.chi.scores<-chisq.test(Ea.chi)
#test result
Ea.chi.scores #p<0.001
Ea.OBS_ON/(Ea.OBS_ON+Ea.OBS_OFF) #proportion observed ON
#Ea.OBS_OFF/(Ea.OBS_ON+Ea.OBS_OFF) #proportion observed OFF
HD.Ea/HD.sum #proportion expected ON



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

#####################################
####Table 1; Environmental Survey Summary ####
#####################################
#Stats based on all quadrats surveyed 
  #(many outside C.tullia dispersal polygons)
H_D <-read.csv("Astley_Habitat_Dispersal_Data.csv")
names(H_D)
nrow(H_D)
#Remove duplicates 
H_D.S <-H_D[!duplicated(H_D[c("Number")]), ]
nrow(H_D.S)
#H_D.S - survey quadrats
nrow(H_D.S)
names(H_D.S)
summary(H_D.S$E.tetralix)
summary(H_D.S$E.vaginatum)
summary(H_D.S$Ev_tussock_count)
summary(H_D.S$Bryophytes_mixed)
summary(H_D.S$C.vulgaris)
summary(H_D.S$E.angustifolium)
summary(H_D.S$Sphagnum_total)
summary(H_D$M.caerulea)
summary(H_D.S$EC)
summary(H_D.S$ORP)
##Subset for CGB_1
H_D1.S <-subset(H_D.S,H_D.S$CGB_1_0==1)
nrow(H_D1.S)
summary(H_D1.S$E.tetralix)
summary(H_D1.S$E.vaginatum)
summary(H_D1.S$Ev_tussock_count)
summary(H_D1.S$Bryophytes_mixed)
summary(H_D1.S$C.vulgaris)
summary(H_D1.S$E.angustifolium)
summary(H_D1.S$Sphagnum_total)
summary(H_D1.S$M.caerulea)
summary(H_D1.S$EC)
summary(H_D1.S$ORP)
##Subset for CGB_0
H_D0.S <-subset(H_D.S,H_D.S$CGB_1_0==0)
nrow(H_D0.S)
summary(H_D0.S$E.tetralix)
summary(H_D0.S$E.vaginatum)
summary(H_D0.S$Ev_tussock_count)
summary(H_D0.S$Bryophytes_mixed)
summary(H_D0.S$C.vulgaris)
summary(H_D0.S$E.angustifolium)
summary(H_D0.S$Sphagnum_total)
summary(H_D0.S$M.caerulea)
summary(H_D0.S$EC)
summary(H_D0.S$ORP)

##Test for difference##
#test for normality
shapiro.test(H_D0.S$Sphagnum_total)
shapiro.test(H_D1.S$Sphagnum_total)
#data not normally distributed

#mean Sphagnum cover quoted in discussion
mean(H_D0.S$Sphagnum_total)


#Global test for difference in plant species composition
  #between cottongrass beds and Molinea tussock areas
#PERMANOVA
library(vegan)
names(H_D.S)
adonis(H_D.S[,7:13]~H_D.S$CGB_1_0, data=H_D.S, method='bray')

#Pairwise comparison of individual variables
wilcox.test(H_D0.S$Sphagnum_total,H_D1.S$Sphagnum_total)
wilcox.test(H_D0.S$Bryophytes_mixed,H_D1.S$Bryophytes_mixed)
wilcox.test(H_D0.S$E.tetralix,H_D1.S$E.tetralix)
wilcox.test(H_D0.S$E.vaginatum,H_D1.S$E.vaginatum)
wilcox.test(H_D0.S$E.angustifolium,H_D1.S$E.angustifolium)
wilcox.test(H_D0.S$M.caerulea,H_D1.S$M.caerulea)
wilcox.test(H_D0.S$C.vulgaris,H_D1.S$C.vulgaris)
wilcox.test(H_D0.S$Ev_tussock_count,H_D1.S$Ev_tussock_count)
wilcox.test(H_D0.S$EC,H_D1.S$EC)
wilcox.test(H_D0.S$ORP,H_D1.S$ORP)