## Conduct accuracy assessment on 2016 NLCD data using Google Earth verified
## iNaturalist records
# Specify  directory for saved functions
#FunctDirect <- "D:/My Drive/win-library/"
FunctDirect <- "C:/Users/jamesltracy61/Documents/R/win-library/"
setwd(FunctDirect)
# Read in User Defined Functions
source("ClassContingencyTable2_Function.R")
source("Kappa.Calc_Function.R")  # used in ClassContingencyTable_Function.R
source("Round2_Function.R")
##
# Load needed packages
#install.packages(c("maptools", "rgdal", "sp", "raster", "plyr","foreign","dplyr","ggpubr","data.table","electivity","reshape","ggcorrplot","psych")) 
library(plyr)
library(foreign)
library(dplyr)
library(ggpubr)
library(data.table)
library(reshape2)
library(ggcorrplot)
library(psych)
library(tidyr)
library(rcompanion)
###
InDirect <- "K:/Monarch/"
Units <- "km"
SpatFiltBuff <- 1 # units in km for spatial filter
#
####################################################################################
### Statistically compare percent utilizaton of habitats among species
setwd(InDirect)
#
SpeciesList <- c("A_cap","A_vir","A_oeno","A_lat")
DataList <- list()
Count <- 0
for(Species in SpeciesList) {
  PresenceDatComb.df4a <- data.frame(read.csv(paste0(Species, "_Thinned", SpatFiltBuff, Units,"_presNLCD_", "1kPopDenIndexdat.csv"), stringsAsFactors=FALSE))
  #
  head(PresenceDatComb.df4a)
  Count <- Count + 1
  if(Count==1) {
    PresenceDatComb.df4a <- PresenceDatComb.df4a[,-1] 
    colnames(PresenceDatComb.df4a)[1] <- "X" 
  }
  PresenceDatComb.df4a$Species <- Species
  DataList[[Count]] <- PresenceDatComb.df4a
  #
}
#head(DataList[[1]])
#ncol(DataList[[4]])
#colnames(DataList[[1]])
##
### Find columns in first but not second data frame
#NewCols.df <- as.data.frame(setdiff(colnames(DataList[[2]]), colnames(DataList[[3]])), stringsAsFactors=FALSE)
##
PresenceDatComb.df4 <- do.call(rbind, DataList)
#
unique(PresenceDatComb.df4$Habitat)
# Combine Trail, Utility, and RR ROWs
PresenceDatComb.df4$Habitat[grepl("Trail|Utility", PresenceDatComb.df4$Habitat)] <- "Trail/Utility/RR ROW"
#
head(PresenceDatComb.df4)
############
# Construct data of contingency table for species versus habitat
HabitatTable <- aggregate(PopRoadHBRWtN ~ Species+Habitat, PresenceDatComb.df4, FUN=function(x) c(CountSum = sum(x)))
#str(HabitatTable)
#
# Transform long to wide table using tidyr
CrossTab.df <- as.data.frame(pivot_wider(HabitatTable, names_from = Habitat, values_from = PopRoadHBRWtN)) 
# Replace na with zero
CrossTab.df[is.na(CrossTab.df)]=0
#str(CrossTab.df)
CrossTab.df2 <- CrossTab.df[,2:(ncol(CrossTab.df))]
rownames(CrossTab.df2) <- CrossTab.df$Species
rowSums(CrossTab.df2)
#
# Run Pearson's Chi-squared test on table
ChiSqResults <- chisq.test(CrossTab.df2)
ChiSqResults$expected
#
CrossTabFreq <- as.data.frame(prop.table(as.matrix(CrossTab.df2), 1))
rowSums(CrossTabFreq)
#
########
#############################################################
## Combine Habitat into Road ROW versus non-Road and analyze
PresenceDatComb.df5 <- PresenceDatComb.df4
# Combine all Habitats with "Road"
PresenceDatComb.df5$Habitat[grepl("Road", PresenceDatComb.df5$Habitat)] <- "Road ROW"
# Combine all Habitats without "ROW" and "Yard"
PresenceDatComb.df5$Habitat[!grepl("Road", PresenceDatComb.df5$Habitat)] <- "non-Road"
#
head(PresenceDatComb.df5)
#################
# Construct data of contingency table for species versus habitat
HabitatTable <- aggregate(PopRoadHBRWtN ~ Species+Habitat, PresenceDatComb.df5, FUN=function(x) c(CountSum = sum(x)))
#str(HabitatTable)
#
# Transform long to wide table using tidyr
CrossTab.df <- as.data.frame(pivot_wider(HabitatTable, names_from = Habitat, values_from = PopRoadHBRWtN)) 
# Replace na with zero
CrossTab.df[is.na(CrossTab.df)]=0
#str(CrossTab.df)
CrossTab.df2 <- CrossTab.df[,2:(ncol(CrossTab.df))]
rownames(CrossTab.df2) <- CrossTab.df$Species
rowSums(CrossTab.df2)
#
########
# Run Pearson's Chi-squared test on table
ChiSqResults <- chisq.test(CrossTab.df2)
ChiSqResults$expected
###
### Conduct Post Hoc tests
CrossTab.mat <- as.matrix(CrossTab.df2)
PostHocChiSqResults <- pairwiseNominalIndependence(CrossTab.mat,
                            fisher = FALSE,
                            gtest  = FALSE,
                            chisq  = TRUE,
                            method = "fdr")
#
FUN = function(i,j){    
      chisq.test(matrix(c(CrossTab.mat[i,1], CrossTab.mat[i,2],
                          CrossTab.mat[j,1], CrossTab.mat[j,2]),
                 nrow=2,
                 byrow=TRUE))$ p.value
                }
#     
PairwisePostHocChiSqResults <- pairwise.table(FUN,
               rownames(CrossTab.mat),
               p.adjust.method="holm")                            
#
# Save above correlation output
out <- paste0("   ")
cat(out, file=paste0("MilkweedSpeciesRoadROWVsNon-Road_ChiSquareIndependence.txt"), sep="\n")
out1 <- paste0("\n", "p-value with Holm Correction for Post-hoc Chi-Square Independence: ")
cat(out1, file=paste0("MilkweedSpeciesRoadROWVsNon-Road_ChiSquareIndependence.txt"), sep="\n", append=TRUE)
out2 <-capture.output(print(PairwisePostHocChiSqResults ))
cat(out2, file=paste0("MilkweedSpeciesRoadROWVsNon-Road_ChiSquareIndependence.txt"), sep="\n", append=TRUE)
out3 <- paste0("   ")
cat(out3, file=paste0("MilkweedSpeciesRoadROWVsNon-Road_ChiSquareIndependence.txt"), sep="\n", append=TRUE)
#
CrossTabFreq <- as.data.frame(prop.table(as.matrix(CrossTab.df2), 1))
rowSums(CrossTabFreq)
# Save CrossTabFreq
write.table(CrossTabFreq, sep= ",", col.names=NA, file=paste0("MilkweedSpeciesRoadROWVsNon-RoadFrequencyPercent.csv"))



