## This program spatio-temporally thins iNaturalist milkweed data from specified 
## dates to within 30 km for unique week/year combinations for phenological analysis
library(maptools)
library(rgdal)
library(sp)
library(raster)
library(geosphere)
library(spThin)
library(maptools)
library(plyr)
library(foreign)
library(dplyr)
library(ggpubr)
library(data.table)
library(electivity)
library(reshape2)
library(ggcorrplot)
library(psych)
library(lubridate)
#
FunctDirect <- "C:/Users/jamesltracy61/Documents/R/win-library/"
setwd(FunctDirect)
#
# Read in User Defined Functions
source("EvalStatVars_Summary_Function2.R")
source("Round2_Function.R")
#
InDirect <- "K:/Monarch"
setwd(InDirect)
Species <- "Danaus_plexippus"
SpatFiltBuff <- 10 # km thinned data
Units <- "km"  # units SpatFiltBuff
StudyAreaName <- paste0("SCStatesEcoregionsDissGeo")
PopGridName <- paste0("sc2010pop1k")
## Read in population density grid
PopDenRaster <- raster(PopGridName)
RandomPointsName <- paste0("NectarPlants_RandomPoints_Thinned10km_100000")
TargetGroupDates <- "1Jan2011_31Dec2020" # date of downloading iNaturalist data
TargetGroupList <- c("Asteraceae", "Apocynaceae", "Lamiaceae", "Verbenaceae")
HostList1 <- c("Asclepias viridis", "Asclepias asperula", "Asclepias oenotheroides", "Asclepias latifolia", "Cynanchum laeve", "Asclepias tuberosa","Asclepias syriaca", "Asclepias perennis", "Asclepias texana", "Asclepias variegata", "Asclepias purpurascens", "Asclepias subverticillata", "Asclepias linearis", "Asclepias incarnata")
HostList2 <- c("Asclepias amplexicaulis", "Asclepias hirtella", "Asclepias speciosa", "Asclepias sullivantii", "Asclepias longifolia", "Asclepias viridiflora", "Asclepias engelmanniana", "Asclepias arenaria", "Asclepias brachystephana", "Asclepias quadrifolia", "Matelea reticulata")
HostList <- c(HostList1, HostList2)
#
# Read in data
HostPlantDat.df <- data.frame(read.csv(paste0(Species, "_PhenologyThinned10kmSC2011_20_AllNativeMilkweeds.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(HostPlantDat.df)
head(HostPlantDat.df)
tail(HostPlantDat.df)
nrow(HostPlantDat.df)
# Clip data to study area
# Read in study area
StudyAreaShp <- readOGR(".", StudyAreaName)
crs(StudyAreaShp)
# Convert data to spatial points data frame
xy <- HostPlantDat.df[,c("lon","lat")]
HostPlantDat.spdf <- SpatialPointsDataFrame(coords = xy, data = HostPlantDat.df, proj4string = crs(StudyAreaShp))
head(HostPlantDat.spdf)
nrow(HostPlantDat.spdf)
#
dev.new()
plot(StudyAreaShp, col='light green')
plot(HostPlantDat.spdf, add=TRUE, axes=TRUE, col='red')
# Clip pres.spdf by StudyAreaShp
HostPlantDat.spdfc <- HostPlantDat.spdf[StudyAreaShp,]
plot(HostPlantDat.spdfc, add=TRUE, axes=TRUE, col= 'blue')
#
# Add population data
plot(PopDenRaster, add=TRUE, axes=TRUE)
HostPlantDat.spdfc1 <- raster::extract(PopDenRaster, HostPlantDat.spdfc, sp=TRUE)
head(HostPlantDat.spdfc1)
nrow(HostPlantDat.spdfc1)
# Convert back to data frame
HostPlantDat.df2 <- as.data.frame(HostPlantDat.spdfc1, stringsAsFactors=FALSE)
head(HostPlantDat.df2)
# Delete last two columns
HostPlantDat.df2 <- HostPlantDat.df2[,-c((ncol(HostPlantDat.df2)-1):(ncol(HostPlantDat.df2)))]
nrow(HostPlantDat.df2)
# Rename last column to "PopDen"
colnames(HostPlantDat.df2)[ncol(HostPlantDat.df2)] <- "PopDen"
nrow(HostPlantDat.df2)
# Omit any rows with NA values for PopDen
# Make function to delete rows with NAs in a specified column
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}
HostPlantDat.df2 <- completeFun(HostPlantDat.df2, "PopDen")
nrow(HostPlantDat.df2)
head(HostPlantDat.df2)
#
HostPlantDat.df2$RawCount <- 1
head(HostPlantDat.df2)
nrow(HostPlantDat.df2)
#
write.table(HostPlantDat.df2, file=paste0(Species, "_PhenologyThinned10kmSCClipped2011_20_AllNativeMilkweeds.csv"), sep=",", row.names=FALSE)
##
#Find clipped points and save in separate file
HostPlantDat.df3 <- HostPlantDat.df2[,c(1:(ncol(HostPlantDat.df2)-2))]
nrow(HostPlantDat.df3)
head(HostPlantDat.df3)
nrow(HostPlantDat.df)
HostPlantDatCombined.df <- rbind(HostPlantDat.df, HostPlantDat.df3)
nrow(HostPlantDatCombined.df)
OmittedDat.df <- setdiff(HostPlantDat.df, HostPlantDat.df3)
nrow(OmittedDat.df)
#
write.table(OmittedDat.df, file=paste0(Species, "_PhenologyThinned10kmSCOmittedinClip2011_20_AllNativeMilkweeds.csv"), sep=",", row.names=FALSE)
#
#######################################
## Read in and process Target Group data for sample bias correction
##
GroupDataList <- list()
Count <- 0
for(Group in TargetGroupList) {
  #Group <- "Asteraceae"
  GroupDataRaw.df <- data.frame(read.csv(paste0(Group, "_iNaturalist_", TargetGroupDates, "_TXOKARLARaw.csv"), stringsAsFactors=FALSE))
  head(GroupDataRaw.df)
  GroupDataRaw.df$Source <- "iNat"
  #str(GroupDataRaw.df)
  GroupDataRaw.df$Year <-  substr(GroupDataRaw.df$observed_on, nchar(GroupDataRaw.df$observed_on)-4+1, nchar(GroupDataRaw.df$observed_on))
  GroupDataRaw.df$Year <- as.numeric(GroupDataRaw.df$Year)
  #
  GroupDataRaw.df$family <- Group
  GroupData.df1 <- GroupDataRaw.df[,c("Source", "id", "latitude", "longitude", "positional_accuracy", "user_id", "observed_on", "Year", "scientific_name", "family", "place_guess", "url", "geoprivacy")]
  head(GroupData.df1)
  nrow(GroupData.df1)
  # Keep only records with no geoprivacy
  GroupData.df2 <- GroupData.df1[ which(GroupData.df1$geoprivacy==""),]
  head(GroupData.df2)
  nrow(GroupData.df2)
  # Delete records with any missing geocoordinates
  GroupData.df2 <- GroupData.df2[complete.cases(GroupData.df2[,c("latitude","longitude")]),]
  #str(GroupData.df2)
  Count <- Count + 1
  GroupDataList[[Count]] <- GroupData.df2
}
GroupDataAll.df <- do.call(rbind, GroupDataList)
head(GroupDataAll.df)
tail(GroupDataAll.df)
nrow(GroupDataAll.df)
colnames(GroupDataAll.df) <- gsub("latitude", "lat", colnames(GroupDataAll.df))
colnames(GroupDataAll.df) <- gsub("longitude", "lon", colnames(GroupDataAll.df))
colnames(GroupDataAll.df) <- gsub("id", "ID", colnames(GroupDataAll.df))
#
# Convert to spatial points data frame
xy <- GroupDataAll.df[,c("lon","lat")]
GroupData.spdf <- SpatialPointsDataFrame(coords = xy, data = GroupDataAll.df, proj4string = crs(StudyAreaShp))
head(GroupData.spdf)
nrow(GroupData.spdf)
# Clip data to study area
GroupData.spdfc <- GroupData.spdf[StudyAreaShp,]
#
dev.new()
plot(StudyAreaShp, col='light green')
plot(GroupData.spdfc, add=TRUE, axes=TRUE, col='red')
# Add population data
## Read in population density grid
GroupData.spdfc1 <- raster::extract(PopDenRaster, GroupData.spdfc, sp=TRUE)
head(GroupData.spdfc1)
nrow(GroupData.spdfc1)
# Convert back to data frame
GroupData.df2 <- as.data.frame(GroupData.spdfc1, stringsAsFactors=FALSE)
head(GroupData.df2)
# Delete last two columns
GroupData.df2 <- GroupData.df2[,-c((ncol(GroupData.df2)-1):(ncol(GroupData.df2)))]
nrow(GroupData.df2)
# Rename last column to "PopDen"
colnames(GroupData.df2)[ncol(GroupData.df2)] <- "PopDen"
nrow(GroupData.df2)
# Omit any rows with NA values for PopDen
# Make function to delete rows with NAs in a specified column
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}
GroupData.df2 <- completeFun(GroupData.df2, "PopDen")
nrow(GroupData.df2)
head(GroupData.df2)
#
GroupData.df2$RawCount <- 1
head(GroupData.df2)
nrow(GroupData.df2)
##########################
# Add needed columns to match larval host data
GroupData.df2$State.Prov <- "."
GroupData.df2$Comments <- "."
GroupData.df2$Date <- GroupData.df2$observed_on
GroupData.df2$Host <- GroupData.df2$scientific_name
WeekNum <- function(Date) {
  WeekNum <- as.numeric(format(as.Date(Date, format="%m/%d/%Y"),"%U"))+1
  return(WeekNum)
}
GroupData.df2$WeekNum <- as.numeric(lapply(GroupData.df2$Date, FUN= WeekNum))
GroupData.df2$Town <- GroupData.df2$place_guess
GroupData.df2$Phenological.Stage <- "."
GroupData.df2$Accuracy <- GroupData.df2$positional_accuracy
GroupData.df2$Number <- 0
GroupData.df2$NumberMilkweeds <- "."
GroupData.df2$MLMPID <- "."
head(GroupData.df2)
nrow(GroupData.df2)
col_order <- colnames(HostPlantDat.df2)
GroupData.df3 <- GroupData.df2[, col_order]
head(GroupData.df3)
nrow(GroupData.df3)
#################
## Extract known host species from data
HostPlantAvail.df <- GroupData.df3[GroupData.df3$Host %in% HostList, ]
head(HostPlantAvail.df)
sort(unique(HostPlantAvail.df$Host))
nrow(HostPlantAvail.df)
# Combine with larval host data
HostPlantAvail.df2 <- rbind(HostPlantAvail.df, HostPlantDat.df2)
nrow(HostPlantAvail.df2)
# Replace Asclepias with A.
HostPlantAvail.df2$Host <- sub("Asclepias", "A.", HostPlantAvail.df2$Host)
head(HostPlantAvail.df2)
HostPlantAvail.df2$PopDen
unique(HostPlantAvail.df2$Host)
# Write data
write.table(HostPlantAvail.df2, file=paste0(Species, "_HostPlantAvail_UnThinned_.csv"), sep=",", row.names=FALSE)

##################
### Thin data by 10 km by species
# Randomize order of data
set.seed(42)
rows <- sample(nrow(HostPlantAvail.df2))
HostPlantAvail.df3 <- HostPlantAvail.df2[rows, ]
#
## Obtain frequency of plants by species
HostFreq.df <- as.data.frame(table(HostPlantAvail.df2$Host), stringsAsFactors=FALSE)
colnames(HostFreq.df)[1] <- c("Host")
HostFreq.df <- HostFreq.df[order(-HostFreq.df$Freq),]
#
# Thin data by 10 km within s
HostPlantAvailList <- list()
SpeciesList <- unique(HostPlantAvail.df2$Host)
length(SpeciesList)
#
Count <- 0
#
for(SpeciesName in SpeciesList) {
  #SpeciesName <- "A. variegata"
  HostPlantAvail.df2a <- HostPlantAvail.df2[which(HostPlantAvail.df2$Host==SpeciesName),]
  nrow(HostPlantAvail.df2a)
  # Randomize order of rows
  set.seed(42)
  rows <- sample(nrow(HostPlantAvail.df2a))
  HostPlantAvail.df3 <- HostPlantAvail.df2a[rows, ]
  # Remove duplicate values for Lat/Lon
  HostPlantAvail.df4 <- HostPlantAvail.df3[!duplicated(HostPlantAvail.df3[c("lat", "lon")]),]
  nrow(HostPlantAvail.df4)
  if(nrow(HostPlantAvail.df4)< 3) {
    if(nrow(HostPlantAvail.df4)==2) {
      # calculate distance between the two points
      PointDist <- distm(c(HostPlantAvail.df4$lon[1], HostPlantAvail.df4$lat[1]), c(HostPlantAvail.df4$lon[2], HostPlantAvail.df4$lat[2]), fun = distHaversine)
      # Delete second record if closer than 10 km (10,000 m)
      if(PointDist < 10000) {
         HostPlantAvail.df5 <- HostPlantAvail.df4[1,]
      } else {
         HostPlantAvail.df5 <- HostPlantAvail.df4
      }
    } else {
      HostPlantAvail.df5 <- HostPlantAvail.df4
    }
  } else {
    # Spatially thin to 10 km
    # Select only lat lon and id columns
    pres.df <- HostPlantAvail.df4[,c("lat","lon", "ID")]
    head(pres.df)
    nrow(pres.df)
    # Omit any rows with NA values
    pres.df <- na.omit(pres.df)
    pres.df$Species <- Species
    ## Spatially thin presence points to 10 km
    presthin.df <- as.data.frame(thin(pres.df, lat.col="lat", long.col="lon", spec.col = "Species", thin.par=10, reps=1, write.file=FALSE, locs.thinned.list.return = TRUE))
    colnames(presthin.df) <- c("lon", "lat")
    # Merge data with original roadkill by lat/long
    HostPlantAvail.df5 <- merge(presthin.df, HostPlantAvail.df4, by=c("lon", "lat"))
    head(HostPlantAvail.df5)
    nrow(HostPlantAvail.df5)
  }
  # Save data to list
  Count <- Count + 1
  HostPlantAvailList[[Count]] <- HostPlantAvail.df5
}  
##########################
HostPlantAvailAll.df <- do.call(rbind, HostPlantAvailList)
colnames(HostPlantAvailAll.df)
head(HostPlantAvailAll.df)
tail(HostPlantAvailAll.df)
unique(HostPlantAvailAll.df$Host)
# Re-order columns
HostPlantAvailAll.df <- HostPlantAvailAll.df[, c(3,4,1,2,5:ncol(HostPlantAvailAll.df))]
#Delete column ID.1
#HostPlantAvailAll.df <- within(HostPlantAvailAll.df, rm("ID.1"))
head(HostPlantAvailAll.df)
tail(HostPlantAvailAll.df)
nrow(HostPlantAvailAll.df)
#
## Obtain frequency of plants by species
HostFreq.df <- as.data.frame(table(HostPlantAvailAll.df$Host), stringsAsFactors=FALSE)
colnames(HostFreq.df)[1] <- c("Host")
HostFreq.df <- HostFreq.df[order(-HostFreq.df$Freq),]
#
write.table(HostPlantAvailAll.df, file=paste0(Species, "_HostPlantAvail_Thinned10km_.csv"), sep=",", row.names=FALSE)
##
###################################################
## Look through and tag "Ornamental" in Comments for common garden species of 
## A. speciosa, A. syriaca, A. tuberosa, and A. incarnata
##
####################################################
## Write target group data to shapefiles
# First, convert to spatial points data frame
xy <- GroupData.df3[,c("lon","lat")]
GroupData.spdf <- SpatialPointsDataFrame(coords = xy, data = GroupData.df3, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
head(GroupData.spdf)
nrow(GroupData.spdf)
##
# Save unthinned target data points that are clipped to study area
writeOGR(GroupData.spdf, ".", paste0(Species, "_HostTargetGroup_Clipped_UnThinned"), driver="ESRI Shapefile", overwrite=TRUE)
#
####################################################
### In ArcGIS, use script MonarchLarvaeLandscapeHostPlantSelectivityAssessment_Processing.py
### to spatially thin target group points with buffers to 10 km (too many points for spthin).
### Also use CreateRandomPoints tool to create 100,000 random points spaced 10 km apart in study area
#####################################################
#
#####################################################
## Read in data for population weighting
HostPlantDat.df2 <- data.frame(read.csv(paste0(Species, "_PhenologyThinned10kmSCClipped2011_20_AllNativeMilkweeds.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(HostPlantDat.df2)
head(HostPlantDat.df2)
tail(HostPlantDat.df2)
nrow(HostPlantDat.df2)
#
HostPlantAvailAll.df <- data.frame(read.csv(paste0(Species, "_HostPlantAvail_Thinned10km_.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(HostPlantAvailAll.df)
head(HostPlantAvailAll.df)
tail(HostPlantAvailAll.df)
nrow(HostPlantAvailAll.df)
# Delete Cynanchum unifarium
HostPlantAvailAll.df <- HostPlantAvailAll.df[!(HostPlantAvailAll.df$Host=="Cynanchum unifarium"), ]
#
## Remove "Ornamental"
# Delete records with comment of "Ornamental"
HostPlantAvailAll.df <- HostPlantAvailAll.df[!grepl("Ornamental", HostPlantAvailAll.df$Comments), , drop = FALSE] 
head(HostPlantAvailAll.df)
nrow(HostPlantAvailAll.df)
#
### Read in thinned target group points and random points shapefiles and convert to data frames
## 
RandomDat.shp <- readOGR(".", "NectarPlants_RandomPoints_Thinned10km_100000")
# Extract PopDen data
RandomDat.spdf1 <- raster::extract(PopDenRaster, RandomDat.shp, sp=TRUE)
head(RandomDat.spdf1)
#str(RandomDat.spdf1)
# Convert back to data frame
RandomDat.df1 <- as.data.frame(RandomDat.spdf1, stringsAsFactors=FALSE)
head(RandomDat.df1)
nrow(RandomDat.df1)
RandomDat.df1 <- na.omit(RandomDat.df1)
#
RandomDat.df1$ID <- rownames(RandomDat.df1)
colnames(RandomDat.df1) <- c("CID", "PopDen", "lon", "lat", "ID")
RandomDat.df2 <- RandomDat.df1[,c("ID", "lon", "lat", "PopDen")]
RandomDat.df2$RawCount <- 1
head(RandomDat.df2)
#
TargetGroupDat.shp <- readOGR(".", paste0(Species, "_HostTargetGroup_Clipped_Thinned_10Kilometers"))
TargetGroupDat.df2 <- as.data.frame(TargetGroupDat.shp)
head(TargetGroupDat.df2)
# Delete last three columns
TargetGroupDat.df3 <- TargetGroupDat.df2[, -c((ncol(TargetGroupDat.df2)-3):ncol(TargetGroupDat.df2))]
head(TargetGroupDat.df3)
colnames(TargetGroupDat.df3) <- sub("RawCont", "RawCount", colnames(TargetGroupDat.df3))
nrow(TargetGroupDat.df3)
TargetGroupDat.df3$PopDen
#
nrow(HostPlantDat.df2)
HostPlantDat.df2$PopDen

#
##################################################################################################
### Histogram Bin Ratio Weighting (HBRWt) for correcting sample bias for thinned host plant records
### Calculate histograms of PopDen for random, target group, and species with same bins. Divide target group
### percentage into random percentage for each bin to obtain a correction factor to weight
### the species count in the bin. Use zero weight, 50% weight and 100% weight to species and compare 
### results
################################################################################################### 
##  First normalize all values of PopDen across all three datasets
# Find the minimum and maximum values from all PopDen data to use 
# for normalization from zero to one
PopDenMin <- min(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, HostPlantDat.df2$PopDen)
PopDenMax <- max(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, HostPlantDat.df2$PopDen)
# Calculate Normalized PopDen, PopDenN
RandomDat.df2$PopDenN <- (RandomDat.df2$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
head(RandomDat.df2)
TargetGroupDat.df3$PopDenN <- (TargetGroupDat.df3$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
HostPlantDat.df2$PopDenN <- (HostPlantDat.df2$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
## Prepare histogram for random points PopDen
png(file = paste0(Species, "_RandomPoints_PopDenNHistogram.png"))
randomhistres <- hist(RandomDat.df2$PopDenN, breaks=seq(0,1,by=0.001))
plot(randomhistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
## Prepare histogram for target group PopDen
png(file = paste0(Species, "_TargetGroup_PopDenNHistogram.png"))
targethistres <- hist(TargetGroupDat.df3$PopDenN, breaks=seq(0,1,by=0.001))
plot(targethistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
## Prepare histogram for presence group PopDen
png(file = paste0(Species, "_Presence_PopDenNHistogram.png"))
presencehistres <- hist(HostPlantDat.df2$PopDenN, breaks=seq(0,1,by=0.001))
plot(presencehistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
############
length(targethistres$count)
randomhistres$count
presencehistres$count

## Assemble a dataframe with all histogram results
##
PopDenHBRWt.df <- data.frame(cbind(randomhistres$count, targethistres$count, presencehistres$count))
head(PopDenHBRWt.df)
tail(PopDenHBRWt.df)
colnames(PopDenHBRWt.df) <- c("Random", "TargetGroup", "Presence")
# Calculate frequency percent for each column (column will sum to 1)
PopDenHBRWt.df$RandomFreqPerc <- PopDenHBRWt.df$Random/sum(PopDenHBRWt.df$Random) 
PopDenHBRWt.df$TargetGroupFreqPerc <- PopDenHBRWt.df$TargetGroup/sum(PopDenHBRWt.df$TargetGroup) 
PopDenHBRWt.df$PresenceFreqPerc <- PopDenHBRWt.df$Presence/sum(PopDenHBRWt.df$Presence)
# Replace all values of zero with 50% of minimum value above zero
PopDenHBRWt.df$RandomFreqPerc[PopDenHBRWt.df$RandomFreqPerc==0] <- unique(sort(PopDenHBRWt.df$RandomFreqPerc))[2]*0.5
PopDenHBRWt.df$TargetGroupFreqPerc[PopDenHBRWt.df$TargetGroupFreqPerc==0] <- unique(sort(PopDenHBRWt.df$TargetGroupFreqPerc))[2]*0.5
PopDenHBRWt.df$PresenceFreqPerc[PopDenHBRWt.df$PresenceFreqPerc==0] <- unique(sort(PopDenHBRWt.df$PresenceFreqPerc))[2]*0.5
# Calculate the PopDenHBRWt
PopDenHBRWt.df$PopDenHBRWt <- PopDenHBRWt.df$RandomFreqPerc/PopDenHBRWt.df$TargetGroupFreqPerc
# Replace all values of zero with 50% of minimum value above zero
PopDenHBRWt.df$PopDenHBRWt[PopDenHBRWt.df$PopDenHBRWt==0] <- unique(sort(PopDenHBRWt.df$PopDenHBRWt))[2]*0.5
# Provide bin values
PopDenHBRWt.df$Bin <- seq(0.001,1,by=0.001)
head(PopDenHBRWt.df)
## Assign Presence data PopDenN values to bins corresponding to PopDenHBRWt.df$Bin
HostPlantDat.df2$PopDenNBinCode <- .bincode(HostPlantDat.df2$PopDenN, seq(0,1,by=0.001), include.lowest=TRUE)
head(HostPlantDat.df2)
# Retrieve HBRWt for PopDenN
HostPlantDat.df2$PopDenNHBRWt <- PopDenHBRWt.df[HostPlantDat.df2$PopDenNBinCode,7]
head(HostPlantDat.df2)
PopDenHBRWt.df[1:20,]
############
# Normalize all weighting factors to sum to total raw count for calculating expected counts
TotRawCount <- sum(HostPlantDat.df2$RawCount)
TotPopDenNHBRWt <- sum(HostPlantDat.df2$PopDenNHBRWt)
# Normalize weights to sum to TotRawCount
HostPlantDat.df2$PopDenNHBRWtN <- HostPlantDat.df2$PopDenNHBRWt/(TotPopDenNHBRWt/TotRawCount)
sum(HostPlantDat.df2$PopDenNHBRWtN)
# Write data
write.table(HostPlantDat.df2, file=paste0(Species, "_HostPlantThinned10kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt.csv"), sep=",", row.names=FALSE)
##
###########################
## Obtain frequency of plants by species using PopDenNHBRWtN
EvaluationStats <- HostPlantDat.df2[,c("Host", "PopDenNHBRWtN")]
head(EvaluationStats)
OutName <- paste0(Species, "_HostFrequencyforThinned10kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt.csv")
SortGroups <- c("Host")
ncol(EvaluationStats)
## All statistics, and only statistics, should be at and after the column specified below
StatVarFirstColumn <- 2
OutDirectIn <- InDirect
EvaluationStatsSummary.df <- EvalStatVars.Summary(EvaluationStats, SortGroups, StatVarFirstColumn, OutName, OutDirectIn)
head(EvaluationStatsSummary.df)
##################################################################################


##################################################################################################
### Histogram Bin Ratio Weighting (HBRWt) for correcting sample bias for available host plants
### Calculate histograms of PopDen for random, target group, and available host species with same bins. Divide target group
### percentage into random percentage for each bin to obtain a correction factor to weight
### the species count in the bin. Use zero weight, 50% weight and 100% weight to species and compare
### results
###################################################################################################
##  First normalize all values of PopDen across all three datasets
# Find the minimum and maximum values from all PopDen data to use
# for normalization from zero to one
# Read back in HostPlantAvail Data
PopDenMin <- min(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen,HostPlantAvailAll.df$PopDen)
PopDenMax <- max(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, HostPlantAvailAll.df$PopDen)
# Calculate Normalized PopDen, PopDenN
RandomDat.df2$PopDenN <- (RandomDat.df2$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
head(RandomDat.df2)
TargetGroupDat.df3$PopDenN <- (TargetGroupDat.df3$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
HostPlantAvailAll.df$PopDenN <- (HostPlantAvailAll.df$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
## Prepare histogram for random points PopDen
png(file = paste0(Species, "_RandomPoints_PopDenNHistogram.png"))
randomhistres <- hist(RandomDat.df2$PopDenN, breaks=seq(0,1,by=0.001))
plot(randomhistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
## Prepare histogram for target group PopDen
png(file = paste0(Species, "_TargetGroup_PopDenNHistogram.png"))
targethistres <- hist(TargetGroupDat.df3$PopDenN, breaks=seq(0,1,by=0.001))
plot(targethistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
## Prepare histogram for presence group PopDen
png(file = paste0(Species, "_Presence_PopDenNHistogram.png"))
presencehistres <- hist(HostPlantAvailAll.df$PopDenN, breaks=seq(0,1,by=0.001))
plot(presencehistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
############
length(targethistres$count)
randomhistres$count
presencehistres$count

## Assemble a dataframe with all histogram results
##
PopDenHBRWt.df <- data.frame(cbind(randomhistres$count, targethistres$count, presencehistres$count))
head(PopDenHBRWt.df)
tail(PopDenHBRWt.df)
colnames(PopDenHBRWt.df) <- c("Random", "TargetGroup", "Presence")
# Calculate frequency percent for each column (column will sum to 1)
PopDenHBRWt.df$RandomFreqPerc <- PopDenHBRWt.df$Random/sum(PopDenHBRWt.df$Random)
PopDenHBRWt.df$TargetGroupFreqPerc <- PopDenHBRWt.df$TargetGroup/sum(PopDenHBRWt.df$TargetGroup)
PopDenHBRWt.df$PresenceFreqPerc <- PopDenHBRWt.df$Presence/sum(PopDenHBRWt.df$Presence)
# Replace all values of zero with 50% of minimum value above zero
PopDenHBRWt.df$RandomFreqPerc[PopDenHBRWt.df$RandomFreqPerc==0] <- unique(sort(PopDenHBRWt.df$RandomFreqPerc))[2]*0.5
PopDenHBRWt.df$TargetGroupFreqPerc[PopDenHBRWt.df$TargetGroupFreqPerc==0] <- unique(sort(PopDenHBRWt.df$TargetGroupFreqPerc))[2]*0.5
PopDenHBRWt.df$PresenceFreqPerc[PopDenHBRWt.df$PresenceFreqPerc==0] <- unique(sort(PopDenHBRWt.df$PresenceFreqPerc))[2]*0.5
# Calculate the PopDenHBRWt
PopDenHBRWt.df$PopDenHBRWt <- PopDenHBRWt.df$RandomFreqPerc/PopDenHBRWt.df$TargetGroupFreqPerc
# Replace all values of zero with 50% of minimum value above zero
PopDenHBRWt.df$PopDenHBRWt[PopDenHBRWt.df$PopDenHBRWt==0] <- unique(sort(PopDenHBRWt.df$PopDenHBRWt))[2]*0.5
# Provide bin values
PopDenHBRWt.df$Bin <- seq(0.001,1,by=0.001)
head(PopDenHBRWt.df)
## Assign Presence data PopDenN values to bins corresponding to PopDenHBRWt.df$Bin
HostPlantAvailAll.df$PopDenNBinCode <- .bincode(HostPlantAvailAll.df$PopDenN, seq(0,1,by=0.001), include.lowest=TRUE)
head(HostPlantAvailAll.df)
# Retrieve HBRWt for PopDenN
HostPlantAvailAll.df$PopDenNHBRWt <- PopDenHBRWt.df[HostPlantAvailAll.df$PopDenNBinCode,7]
head(HostPlantAvailAll.df)
PopDenHBRWt.df[1:20,]
############
# Normalize all weighting factors to sum to total raw count for calculating expected counts
TotPopDenNHBRWtAvail <- sum(HostPlantAvailAll.df$PopDenNHBRWt)
# Normalize weights to sum to TotRawCount
HostPlantAvailAll.df$PopDenNHBRWtN <- HostPlantAvailAll.df$PopDenNHBRWt/(TotPopDenNHBRWtAvail/TotRawCount)
sum(HostPlantAvailAll.df$PopDenNHBRWtN)
##
# Write data
write.table(HostPlantAvailAll.df, file=paste0(Species, "_HostPlantAvail_Thinned10km_PopDenHBRWt.csv"), sep=",", row.names=FALSE)
###########################
## Obtain frequency of plants by species using PopDenNHBRWtN
EvaluationStats <- HostPlantAvailAll.df[,c("Host", "PopDenNHBRWtN")]
head(EvaluationStats)
OutName <- paste0(Species, "_HostAvailableFrequencyforThinned10kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt.csv")
SortGroups <- c("Host")
ncol(EvaluationStats)
## All statistics, and only statistics, should be at and after the column specified below
StatVarFirstColumn <- 2
OutDirectIn <- InDirect
EvaluationStatsSummary.df <- EvalStatVars.Summary(EvaluationStats, SortGroups, StatVarFirstColumn, OutName, OutDirectIn)
head(EvaluationStatsSummary.df)
##################################################################################

########################################################################
### Summarize population weighted data 
################################################
## Obtain counts by Host
ActualDatTab.dfa <- aggregate(RawCount ~ Host, data=HostPlantDat.df2, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfb <- aggregate(PopDenNHBRWtN ~ Host, data=HostPlantDat.df2, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfc <- Reduce(function(x,y) merge(x,y, by="Host", all=TRUE), list(ActualDatTab.dfa, ActualDatTab.dfb))
#
sum(ActualDatTab.dfc$RawCount)
sum(ActualDatTab.dfc$PopDenNHBRWtN)
##
#########
ActualDatTab.dfc$RawPercUsed <- ActualDatTab.dfc$RawCount/(sum(ActualDatTab.dfc$RawCount))
ActualDatTab.dfc$PopDenNHBRWtNPercUsed <- ActualDatTab.dfc$PopDenNHBRWtN/(sum(ActualDatTab.dfc$RawCount))
#
###########
## Obtain counts by Host Availability
ActualDatTab.dfd <- aggregate(RawCount ~ Host, data=HostPlantAvailAll.df, FUN=function(x) c(Sum = sum(x)))
colnames(ActualDatTab.dfd)[2] <- "AvailCount"
ActualDatTab.dfe <- aggregate(PopDenNHBRWtN ~ Host, data=HostPlantAvailAll.df, FUN=function(x) c(Sum = sum(x)))
colnames(ActualDatTab.dfe)[2] <- "AvailPopDenNHBRWtN"
ActualDatTab.dff <- Reduce(function(x,y) merge(x,y, by="Host", all=TRUE), list(ActualDatTab.dfd, ActualDatTab.dfe))
#
sum(ActualDatTab.dff$AvailCount)
sum(ActualDatTab.dff$AvailPopDenNHBRWtN)
#
ActualDatTab.dff$RawPercAvail <- ActualDatTab.dff$AvailCount/(sum(ActualDatTab.dff$AvailCount))
ActualDatTab.dff$PopDenNHBRWtNPercAvail <- ActualDatTab.dff$AvailPopDenNHBRWtN/TotRawCount
##
ActualDatTab.df <- Reduce(function(x,y) merge(x,y, by="Host", all=TRUE), list(ActualDatTab.dfc, ActualDatTab.dff))
#########################
#
# Calculate expected count, which is RandomPercAvail * Total Actual PopDenInd across all classes except water
ActualDatTab.df$ExpectedCount <-  ActualDatTab.df$PopDenNHBRWtNPercAvail*TotRawCount
TotExpectedCount <- sum(ActualDatTab.df$ExpectedCount)
#
## Identify and combine classes that have expected counts below 5 (c.f. Bryers et al. 1984, Manly et al. 2002)
CombineDropClasses.df <- ActualDatTab.df[ which(ActualDatTab.df$ExpectedCount<5),]
print(CombineDropClasses.df)
#str(DropClasses.df)
# Save Dropped Classes
write.table(CombineDropClasses.df, sep= ",", col.names=NA, file=paste0(Species, "CombinedHostsforMonarchHostSelectivity.csv"))
#
# Save above list of dropped classes to text file
out1 <- paste0("List of Hosts combined that had expected count less than 5")
cat(out1, file=paste0(Species, "ListCombinedHostsforMonarchHostSelectivity.txt"), sep="\n", append=TRUE)
out2 <- capture.output(print(CombineDropClasses.df))
cat(out2, file=paste0(Species, "ListCombinedHostsforMonarchHostSelectivity.txt"), sep="\n", append=TRUE)
#
CombineDropClasses <- CombineDropClasses.df$Host
## Examine data to determine which similar classes to combine and change names for combined classes
## Combine Hosts for ActualDatTab.df using new HostsN
ActualDatTab.df$HostN <- ActualDatTab.df$Host
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. syriaca"] <- "A. speciosa Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. speciosa"] <- "A. speciosa Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. sullivantii"] <- "A. speciosa Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. arenaria"] <- "A. arenaria Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. brachystephana"] <- "A. arenaria Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. subverticillata"] <- "A. arenaria Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. incarnata"] <- "A. linearis Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. purpurascens"] <- "A. quadrifolia Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. quadrifolia"] <- "A. quadrifolia Grp"
ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. linearis"] <- "A. linearis Grp"
#
write.table(ActualDatTab.df, sep = ",", col.names=NA, file=paste0(Species, "_Thinned", SpatFiltBuff, Units,"_OriginalHostDataCounts.csv"))
#
ActualDatTab.df2 <- ActualDatTab.df
ActualDatTab.df2$Host <- ActualDatTab.df2$HostN
ActualDatTab.df2 <- ActualDatTab.df2[, !(names(ActualDatTab.df2) %in% c("HostN"))]
# Resummarize by Host
Columns <- colnames(ActualDatTab.df2)[2:ncol(ActualDatTab.df2)]
ActualDatTab.df3 <- aggregate(. ~ Host, data=ActualDatTab.df2, FUN=function(x) c(Sum = sum(x)))
#
# Save data
write.table(ActualDatTab.df3, sep = ",", col.names=NA, file=paste0(Species, "_Thinned", SpatFiltBuff, Units,"_GroupedHostDataCounts.csv"))
#

###############################################################################################################
## Calculate Bonferroni Confidence Intervals for whether RandomPercAvail falls within BCI for PercUsed
## based on raw data and bias adjustments for human populations
## following Bryers et al. 1984 and Manly et al. 2002
###############################################################################################################
# Read back in data
ActualDatTab.df3 <- data.frame(read.csv(paste0(Species, "_Thinned", SpatFiltBuff, Units,"_GroupedHostDataCounts.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
## Delete row names column
ActualDatTab.df3 <- ActualDatTab.df3[,-1]
#
sum(ActualDatTab.df3$ExpectedCount)
sum(ActualDatTab.df3$PopDenNHBRWtN)
sum(ActualDatTab.df3$PopDenNHBRWtNPercAvail)
NumberComparisons <- nrow(ActualDatTab.df3)
#str(ActualDatTab.df3)
#
## Loop through the test for the raw and six bias adjustments for percent land cover used
## and calculate electivity indices
UsedColumnList <- seq(2, 3, 1)
UsedCategoryList <- colnames(ActualDatTab.df3)[UsedColumnList]
PercUsedColumnList <- seq(4,5,1)
PercUsedCategoryList <- colnames(ActualDatTab.df3)[PercUsedColumnList]
count=0
for(UsedCol in UsedColumnList) {
  #UsedCol=2
  PercUsedCol <- UsedCol + 2
  count=count+1
  ChiSqGFT_Results <- chisq.test(ActualDatTab.df3[,UsedCol], p=ActualDatTab.df3$PopDenNHBRWtNPercAvail)
  #str(ChiSqGFT_Results)
  #ChiSqGFT_Results$expected
  # Save above test results
  out1 <- paste0("Chi-Square Goodness of Fit Test for Differeces in Expected (Random) Versus ", UsedCategoryList[count]," Utilization of Milkweed Hosts")
  cat(out1, file=paste0(Species, "2011_2020MilkweedHostChiSquareGoodnessFitTestfor", UsedCategoryList[count],".txt"), sep="\n")
  out2 <- capture.output(print(ChiSqGFT_Results))
  cat(out2, file=paste0(Species, "2011_2020MilkweedHostChiSquareGoodnessFitTestfor", UsedCategoryList[count],".txt"), sep="\n", append=TRUE)
  #
  # Assuming above ChiSquare Test is significant, find bonferoni critical z-value of 0.05 with NumberComparisons categories for two-tailed 95% CI
  pc <- 0.05/(NumberComparisons)
  za <- qnorm(1-(pc/2))
  ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCI")] <- za*sqrt((ActualDatTab.df3[,PercUsedCol]*(1-ActualDatTab.df3[,PercUsedCol]))/sum(ActualDatTab.df3$ExpectedCount))
  ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"LowerBCI")] <- ActualDatTab.df3[,PercUsedCol] - ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCI")]
  ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"UpperBCI")] <- ActualDatTab.df3[,PercUsedCol] + ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCI")]
  # If RandomPercAvail outside BCI, identify as Significant difference
  if(ChiSqGFT_Results$p.value<0.05) {
    ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCISig")] <- ifelse(ActualDatTab.df3$PopDenNHBRWtNPercAvail <= ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"UpperBCI")] & ActualDatTab.df3$PopDenNHBRWtNPercAvail >=  ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"LowerBCI")],  "NonSig",  "Significant")
    } else {
    ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCISig")] <- "NonSig"
    #
  }
  ActualDatTab.df3[, paste0("ForateRat", PercUsedCategoryList[count])] <- ivlev_forage(ActualDatTab.df3[,PercUsedCol], ActualDatTab.df3$PopDenNHBRWtNPercAvail, log10=FALSE)
  ActualDatTab.df3[, paste0("VS_Elect", PercUsedCategoryList[count])] <- vs_electivity(ActualDatTab.df3[,PercUsedCol], ActualDatTab.df3$PopDenNHBRWtNPercAvail)
  ActualDatTab.df3[, paste0("VS_Elect", PercUsedCategoryList[count], "_10K")] <- Round2(ActualDatTab.df3[, paste0("VS_Elect", PercUsedCategoryList[count])]*10000,0)
}
## Calculated Normalized VS_Elect
VS_ElectMin <- min(ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsed)
VS_ElectMax <- max(ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsed)
# Calculate Normalized PopDen, PopDenN
ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedN_10K <- 10000*((ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsed - VS_ElectMin)/(VS_ElectMax - VS_ElectMin))
head(ActualDatTab.df3)
# Calculate Percent of Total Normalized VS_Elect
ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedNP <- ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedN_10K/sum(ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedN_10K)
#
## Save data
write.table(ActualDatTab.df3, sep = ",", col.names=NA, file=paste0(Species, "_MilkweedGroupHostActualPercentUsedStatsIndices.csv"))

##################################################################################################################################
### Analyze for fall and spring records. Will need to filter available and used hosts by WeekNum
###################################################################################################################################
##
#####################################################
## Read in data for population weighting
HostPlantDat.df2 <- data.frame(read.csv(paste0(Species, "_PhenologyThinned10kmSCClipped2011_20_AllNativeMilkweeds.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(HostPlantDat.df2)
head(HostPlantDat.df2)
tail(HostPlantDat.df2)
nrow(HostPlantDat.df2)
#
HostPlantAvailAll.df <- data.frame(read.csv(paste0(Species, "_HostPlantAvail_Thinned10km_.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(HostPlantAvailAll.df)
head(HostPlantAvailAll.df)
tail(HostPlantAvailAll.df)
nrow(HostPlantAvailAll.df)
# Delete Cynanchum unifarium
HostPlantAvailAll.df <- HostPlantAvailAll.df[!(HostPlantAvailAll.df$Host=="Cynanchum unifarium"), ]
#
## Remove "Ornamental"
# Delete records with comment of "Ornamental"
HostPlantAvailAll.df <- HostPlantAvailAll.df[!grepl("Ornamental", HostPlantAvailAll.df$Comments), , drop = FALSE] 
head(HostPlantAvailAll.df)
nrow(HostPlantAvailAll.df)
#
### Read in thinned target group points and random points shapefiles and convert to data frames
## 
RandomDat.shp <- readOGR(".", "NectarPlants_RandomPoints_Thinned10km_100000")
# Extract PopDen data
RandomDat.spdf1 <- raster::extract(PopDenRaster, RandomDat.shp, sp=TRUE)
head(RandomDat.spdf1)
#str(RandomDat.spdf1)
# Convert back to data frame
RandomDat.df1 <- as.data.frame(RandomDat.spdf1, stringsAsFactors=FALSE)
head(RandomDat.df1)
nrow(RandomDat.df1)
RandomDat.df1 <- na.omit(RandomDat.df1)
#
RandomDat.df1$ID <- rownames(RandomDat.df1)
colnames(RandomDat.df1) <- c("CID", "PopDen", "lon", "lat", "ID")
RandomDat.df2 <- RandomDat.df1[,c("ID", "lon", "lat", "PopDen")]
RandomDat.df2$RawCount <- 1
head(RandomDat.df2)
#
TargetGroupDat.shp <- readOGR(".", paste0(Species, "_HostTargetGroup_Clipped_Thinned_10Kilometers"))
TargetGroupDat.df2 <- as.data.frame(TargetGroupDat.shp)
head(TargetGroupDat.df2)
# Delete last three columns
TargetGroupDat.df3 <- TargetGroupDat.df2[, -c((ncol(TargetGroupDat.df2)-3):ncol(TargetGroupDat.df2))]
head(TargetGroupDat.df3)
colnames(TargetGroupDat.df3) <- sub("RawCont", "RawCount", colnames(TargetGroupDat.df3))
nrow(TargetGroupDat.df3)
TargetGroupDat.df3$PopDen
#
nrow(HostPlantDat.df2)
HostPlantDat.df2$PopDen

################################################################################################
## Loop through data to filter by early and late season
##
WeekNumMinList <- c(11, 32)
WeekNumMaxList <- c(28, 47)
SeasonList <- c("Spring", "Fall")
for(i in 1:2) {
  #i=1
  WeekNumMin <- WeekNumMinList[i]
  WeekNumMax <- WeekNumMaxList[i]
  Season <- SeasonList[i]
  #
  HostPlantDat.df2a <- HostPlantDat.df2[ which(HostPlantDat.df2$WeekNum>=WeekNumMin & HostPlantDat.df2$WeekNum<=WeekNumMax),]
  nrow(HostPlantDat.df2a)
  HostPlantDat.df2a$WeekNum
  #
  HostPlantAvailAll.dfa <- HostPlantAvailAll.df[ which(HostPlantAvailAll.df$WeekNum>=WeekNumMin & HostPlantAvailAll.df$WeekNum<=WeekNumMax),]
  nrow(HostPlantAvailAll.dfa)
  HostPlantAvailAll.dfa$WeekNum
  #
  ##################################################################################################
  ### Histogram Bin Ratio Weighting (HBRWt) for correcting sample bias for thinned host plant records
  ### Calculate histograms of PopDen for random, target group, and species with same bins. Divide target group
  ### percentage into random percentage for each bin to obtain a correction factor to weight
  ### the species count in the bin. Use zero weight, 50% weight and 100% weight to species and compare
  ### results
  ###################################################################################################
  ##  First normalize all values of PopDen across all three datasets
  # Find the minimum and maximum values from all PopDen data to use
  # for normalization from zero to one
  PopDenMin <- min(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, HostPlantDat.df2a$PopDen)
  PopDenMax <- max(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, HostPlantDat.df2a$PopDen)
  # Calculate Normalized PopDen, PopDenN
  RandomDat.df2$PopDenN <- (RandomDat.df2$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
  head(RandomDat.df2)
  TargetGroupDat.df3$PopDenN <- (TargetGroupDat.df3$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
  HostPlantDat.df2a$PopDenN <- (HostPlantDat.df2a$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
  ## Prepare histogram for random points PopDen
  png(file = paste0(Species, "_", Season, "_RandomPoints_PopDenNHistogram.png"))
  randomhistres <- hist(RandomDat.df2$PopDenN, breaks=seq(0,1,by=0.001))
  plot(randomhistres, col="blue", border= "blue", xlab = "Normalized Population Density")
  dev.off()
  ## Prepare histogram for target group PopDen
  png(file = paste0(Species, "_", Season, "_TargetGroup_PopDenNHistogram.png"))
  targethistres <- hist(TargetGroupDat.df3$PopDenN, breaks=seq(0,1,by=0.001))
  plot(targethistres, col="blue", border= "blue", xlab = "Normalized Population Density")
  dev.off()
  ## Prepare histogram for presence group PopDen
  png(file = paste0(Species, "_", Season, "_Presence_PopDenNHistogram.png"))
  presencehistres <- hist(HostPlantDat.df2a$PopDenN, breaks=seq(0,1,by=0.001))
  plot(presencehistres, col="blue", border= "blue", xlab = "Normalized Population Density")
  dev.off()
  ############
  length(targethistres$count)
  randomhistres$count
  presencehistres$count

  ## Assemble a dataframe with all histogram results
  ##
  PopDenHBRWt.df <- data.frame(cbind(randomhistres$count, targethistres$count, presencehistres$count))
  head(PopDenHBRWt.df)
  tail(PopDenHBRWt.df)
  colnames(PopDenHBRWt.df) <- c("Random", "TargetGroup", "Presence")
  # Calculate frequency percent for each column (column will sum to 1)
  PopDenHBRWt.df$RandomFreqPerc <- PopDenHBRWt.df$Random/sum(PopDenHBRWt.df$Random)
  PopDenHBRWt.df$TargetGroupFreqPerc <- PopDenHBRWt.df$TargetGroup/sum(PopDenHBRWt.df$TargetGroup)
  PopDenHBRWt.df$PresenceFreqPerc <- PopDenHBRWt.df$Presence/sum(PopDenHBRWt.df$Presence)
  # Replace all values of zero with 50% of minimum value above zero
  PopDenHBRWt.df$RandomFreqPerc[PopDenHBRWt.df$RandomFreqPerc==0] <- unique(sort(PopDenHBRWt.df$RandomFreqPerc))[2]*0.5
  PopDenHBRWt.df$TargetGroupFreqPerc[PopDenHBRWt.df$TargetGroupFreqPerc==0] <- unique(sort(PopDenHBRWt.df$TargetGroupFreqPerc))[2]*0.5
  PopDenHBRWt.df$PresenceFreqPerc[PopDenHBRWt.df$PresenceFreqPerc==0] <- unique(sort(PopDenHBRWt.df$PresenceFreqPerc))[2]*0.5
  # Calculate the PopDenHBRWt
  PopDenHBRWt.df$PopDenHBRWt <- PopDenHBRWt.df$RandomFreqPerc/PopDenHBRWt.df$TargetGroupFreqPerc
  # Replace all values of zero with 50% of minimum value above zero
  PopDenHBRWt.df$PopDenHBRWt[PopDenHBRWt.df$PopDenHBRWt==0] <- unique(sort(PopDenHBRWt.df$PopDenHBRWt))[2]*0.5
  # Provide bin values
  PopDenHBRWt.df$Bin <- seq(0.001,1,by=0.001)
  head(PopDenHBRWt.df)
  ## Assign Presence data PopDenN values to bins corresponding to PopDenHBRWt.df$Bin
  HostPlantDat.df2a$PopDenNBinCode <- .bincode(HostPlantDat.df2a$PopDenN, seq(0,1,by=0.001), include.lowest=TRUE)
  head(HostPlantDat.df2a)
  # Retrieve HBRWt for PopDenN
  HostPlantDat.df2a$PopDenNHBRWt <- PopDenHBRWt.df[HostPlantDat.df2a$PopDenNBinCode,7]
  head(HostPlantDat.df2a)
  PopDenHBRWt.df[1:20,]
  ############
  # Normalize all weighting factors to sum to total raw count for calculating expected counts
  TotRawCount <- sum(HostPlantDat.df2a$RawCount)
  TotPopDenNHBRWt <- sum(HostPlantDat.df2a$PopDenNHBRWt)
  # Normalize weights to sum to TotRawCount
  HostPlantDat.df2a$PopDenNHBRWtN <- HostPlantDat.df2a$PopDenNHBRWt/(TotPopDenNHBRWt/TotRawCount)
  sum(HostPlantDat.df2a$PopDenNHBRWtN)
  # Write data
  write.table(HostPlantDat.df2a, file=paste0(Species, "_", Season, "_HostPlantThinned10kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt.csv"), sep=",", row.names=FALSE)
  ##
  ###########################
  ## Obtain frequency of plants by species using PopDenNHBRWtN
  EvaluationStats <- HostPlantDat.df2a[,c("Host", "PopDenNHBRWtN")]
  head(EvaluationStats)
  OutName <- paste0(Species, "_", Season, "_HostFrequencyforThinned10kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt.csv")
  SortGroups <- c("Host")
  ncol(EvaluationStats)
  ## All statistics, and only statistics, should be at and after the column specified below
  StatVarFirstColumn <- 2
  OutDirectIn <- InDirect
  EvaluationStatsSummary.df <- EvalStatVars.Summary(EvaluationStats, SortGroups, StatVarFirstColumn, OutName, OutDirectIn)
  head(EvaluationStatsSummary.df)
  ##################################################################################


  ##################################################################################################
  ### Histogram Bin Ratio Weighting (HBRWt) for correcting sample bias for available host plants
  ### Calculate histograms of PopDen for random, target group, and available host species with same bins. Divide target group
  ### percentage into random percentage for each bin to obtain a correction factor to weight
  ### the species count in the bin. Use zero weight, 50% weight and 100% weight to species and compare
  ### results
  ###################################################################################################
  ##  First normalize all values of PopDen across all three datasets
  # Find the minimum and maximum values from all PopDen data to use
  # for normalization from zero to one
  # Read back in HostPlantAvail Data
  PopDenMin <- min(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen,HostPlantAvailAll.dfa$PopDen)
  PopDenMax <- max(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, HostPlantAvailAll.dfa$PopDen)
  # Calculate Normalized PopDen, PopDenN
  RandomDat.df2$PopDenN <- (RandomDat.df2$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
  head(RandomDat.df2)
  TargetGroupDat.df3$PopDenN <- (TargetGroupDat.df3$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
  HostPlantAvailAll.dfa$PopDenN <- (HostPlantAvailAll.dfa$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
  ## Prepare histogram for random points PopDen
  png(file = paste0(Species, "_", Season, "_RandomPoints_PopDenNHistogram.png"))
  randomhistres <- hist(RandomDat.df2$PopDenN, breaks=seq(0,1,by=0.001))
  plot(randomhistres, col="blue", border= "blue", xlab = "Normalized Population Density")
  dev.off()
  ## Prepare histogram for target group PopDen
  png(file = paste0(Species, "_", Season, "_TargetGroup_PopDenNHistogram.png"))
  targethistres <- hist(TargetGroupDat.df3$PopDenN, breaks=seq(0,1,by=0.001))
  plot(targethistres, col="blue", border= "blue", xlab = "Normalized Population Density")
  dev.off()
  ## Prepare histogram for presence group PopDen
  png(file = paste0(Species, "_", Season, "_Presence_PopDenNHistogram.png"))
  presencehistres <- hist(HostPlantAvailAll.dfa$PopDenN, breaks=seq(0,1,by=0.001))
  plot(presencehistres, col="blue", border= "blue", xlab = "Normalized Population Density")
  dev.off()
  ############
  length(targethistres$count)
  randomhistres$count
  presencehistres$count

  ## Assemble a dataframe with all histogram results
  ##
  PopDenHBRWt.df <- data.frame(cbind(randomhistres$count, targethistres$count, presencehistres$count))
  head(PopDenHBRWt.df)
  tail(PopDenHBRWt.df)
  colnames(PopDenHBRWt.df) <- c("Random", "TargetGroup", "Presence")
  # Calculate frequency percent for each column (column will sum to 1)
  PopDenHBRWt.df$RandomFreqPerc <- PopDenHBRWt.df$Random/sum(PopDenHBRWt.df$Random)
  PopDenHBRWt.df$TargetGroupFreqPerc <- PopDenHBRWt.df$TargetGroup/sum(PopDenHBRWt.df$TargetGroup)
  PopDenHBRWt.df$PresenceFreqPerc <- PopDenHBRWt.df$Presence/sum(PopDenHBRWt.df$Presence)
  # Replace all values of zero with 50% of minimum value above zero
  PopDenHBRWt.df$RandomFreqPerc[PopDenHBRWt.df$RandomFreqPerc==0] <- unique(sort(PopDenHBRWt.df$RandomFreqPerc))[2]*0.5
  PopDenHBRWt.df$TargetGroupFreqPerc[PopDenHBRWt.df$TargetGroupFreqPerc==0] <- unique(sort(PopDenHBRWt.df$TargetGroupFreqPerc))[2]*0.5
  PopDenHBRWt.df$PresenceFreqPerc[PopDenHBRWt.df$PresenceFreqPerc==0] <- unique(sort(PopDenHBRWt.df$PresenceFreqPerc))[2]*0.5
  # Calculate the PopDenHBRWt
  PopDenHBRWt.df$PopDenHBRWt <- PopDenHBRWt.df$RandomFreqPerc/PopDenHBRWt.df$TargetGroupFreqPerc
  # Replace all values of zero with 50% of minimum value above zero
  PopDenHBRWt.df$PopDenHBRWt[PopDenHBRWt.df$PopDenHBRWt==0] <- unique(sort(PopDenHBRWt.df$PopDenHBRWt))[2]*0.5
  # Provide bin values
  PopDenHBRWt.df$Bin <- seq(0.001,1,by=0.001)
  head(PopDenHBRWt.df)
  ## Assign Presence data PopDenN values to bins corresponding to PopDenHBRWt.df$Bin
  HostPlantAvailAll.dfa$PopDenNBinCode <- .bincode(HostPlantAvailAll.dfa$PopDenN, seq(0,1,by=0.001), include.lowest=TRUE)
  head(HostPlantAvailAll.dfa)
  # Retrieve HBRWt for PopDenN
  HostPlantAvailAll.dfa$PopDenNHBRWt <- PopDenHBRWt.df[HostPlantAvailAll.dfa$PopDenNBinCode,7]
  head(HostPlantAvailAll.dfa)
  PopDenHBRWt.df[1:20,]
  ############
  # Normalize all weighting factors to sum to total raw count for calculating expected counts
  TotPopDenNHBRWtAvail <- sum(HostPlantAvailAll.dfa$PopDenNHBRWt)
  # Normalize weights to sum to TotRawCount
  HostPlantAvailAll.dfa$PopDenNHBRWtN <- HostPlantAvailAll.dfa$PopDenNHBRWt/(TotPopDenNHBRWtAvail/TotRawCount)
  sum(HostPlantAvailAll.dfa$PopDenNHBRWtN)
  ##
  # Write data
  write.table(HostPlantAvailAll.dfa, file=paste0(Species, "_", Season, "_HostPlantAvail_Thinned10km_PopDenHBRWt.csv"), sep=",", row.names=FALSE)
  ###########################
  ## Obtain frequency of plants by species using PopDenNHBRWtN
  EvaluationStats <- HostPlantAvailAll.dfa[,c("Host", "PopDenNHBRWtN")]
  head(EvaluationStats)
  OutName <- paste0(Species, "_", Season, "_HostAvailableFrequencyforThinned10kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt.csv")
  SortGroups <- c("Host")
  ncol(EvaluationStats)
  ## All statistics, and only statistics, should be at and after the column specified below
  StatVarFirstColumn <- 2
  OutDirectIn <- InDirect
  EvaluationStatsSummary.df <- EvalStatVars.Summary(EvaluationStats, SortGroups, StatVarFirstColumn, OutName, OutDirectIn)
  head(EvaluationStatsSummary.df)
  ##################################################################################

  ########################################################################
  ### Summarize population weighted data
  ################################################
  ## Obtain counts by Host
  ActualDatTab.dfa <- aggregate(RawCount ~ Host, data=HostPlantDat.df2a, FUN=function(x) c(Sum = sum(x)))
  ActualDatTab.dfb <- aggregate(PopDenNHBRWtN ~ Host, data=HostPlantDat.df2a, FUN=function(x) c(Sum = sum(x)))
  ActualDatTab.dfc <- Reduce(function(x,y) merge(x,y, by="Host", all=TRUE), list(ActualDatTab.dfa, ActualDatTab.dfb))
  #
  sum(ActualDatTab.dfc$RawCount)
  sum(ActualDatTab.dfc$PopDenNHBRWtN)
  ##
  #########
  ActualDatTab.dfc$RawPercUsed <- ActualDatTab.dfc$RawCount/(sum(ActualDatTab.dfc$RawCount))
  ActualDatTab.dfc$PopDenNHBRWtNPercUsed <- ActualDatTab.dfc$PopDenNHBRWtN/(sum(ActualDatTab.dfc$RawCount))
  #
  ###########
  ## Obtain counts by Host Availability
  ActualDatTab.dfd <- aggregate(RawCount ~ Host, data=HostPlantAvailAll.dfa, FUN=function(x) c(Sum = sum(x)))
  colnames(ActualDatTab.dfd)[2] <- "AvailCount"
  ActualDatTab.dfe <- aggregate(PopDenNHBRWtN ~ Host, data=HostPlantAvailAll.dfa, FUN=function(x) c(Sum = sum(x)))
  colnames(ActualDatTab.dfe)[2] <- "AvailPopDenNHBRWtN"
  ActualDatTab.dff <- Reduce(function(x,y) merge(x,y, by="Host", all=TRUE), list(ActualDatTab.dfd, ActualDatTab.dfe))
  #
  sum(ActualDatTab.dff$AvailCount)
  sum(ActualDatTab.dff$AvailPopDenNHBRWtN)
  #
  ActualDatTab.dff$RawPercAvail <- ActualDatTab.dff$AvailCount/(sum(ActualDatTab.dff$AvailCount))
  sum(ActualDatTab.dff$AvailPopDenNHBRWtN)
  ActualDatTab.dff$PopDenNHBRWtNPercAvail <- ActualDatTab.dff$AvailPopDenNHBRWtN/TotRawCount
  sum(ActualDatTab.dff$PopDenNHBRWtNPercAvail)
  ##
  ActualDatTab.df <- Reduce(function(x,y) merge(x,y, by="Host", all=TRUE), list(ActualDatTab.dfc, ActualDatTab.dff))
  sum(ActualDatTab.df$PopDenNHBRWtNPercAvail)
  # Write taxa that were missing in spring host data
  MissingHosts.df <- as.data.frame(ActualDatTab.df[is.na(ActualDatTab.df$RawCount),1])
  colnames(MissingHosts.df) <- "Host"
  # Save Missing Hosts
  write.table(MissingHosts.df, sep= ",", col.names=NA, file=paste0(Species, "_", Season, "MissingHostsforMonarchHostSelectivity.csv"))
  #
  # Omit any rows with missing data
  ActualDatTab.df <- na.omit(ActualDatTab.df)
  sum(ActualDatTab.df$RawCount)
  sum(ActualDatTab.df$PopDenNHBRWtNPercAvail)
  # Recalculate RawPercAvail and PopDenNHBRWtNPercAvail based on remaining seasonal hosts
  sum(ActualDatTab.df$AvailCount)
  ActualDatTab.df$RawPercAvail <- ActualDatTab.df$AvailCount/(sum(ActualDatTab.df$AvailCount))
  # Renormalize PopDenNHBRWtNPercAvail by calculating the product of each value and 1 divided by the sum of values
  ActualDatTab.df$PopDenNHBRWtNPercAvailN <- ActualDatTab.df$PopDenNHBRWtNPercAvail * (1/sum(ActualDatTab.df$PopDenNHBRWtNPercAvail))
  sum(ActualDatTab.df$RawPercAvail)
  sum(ActualDatTab.df$PopDenNHBRWtNPercAvailN)
  ActualDatTab.df$PopDenNHBRWtNPercAvail <- ActualDatTab.df$PopDenNHBRWtNPercAvailN
  sum(ActualDatTab.df$PopDenNHBRWtNPercAvail)
  ActualDatTab.df <- ActualDatTab.df[, !(names(ActualDatTab.df) %in% "PopDenNHBRWtNPercAvailN")] 
  # 
  #########################
  #
  # Calculate expected count, which is RandomPercAvail * Total Actual PopDenInd across all classes except water
  ActualDatTab.df$ExpectedCount <-  ActualDatTab.df$PopDenNHBRWtNPercAvail*TotRawCount
  TotExpectedCount <- sum(ActualDatTab.df$ExpectedCount)
  #
  write.table(ActualDatTab.df, sep = ",", col.names=NA, file=paste0(Species, "_", Season, "_Thinned", SpatFiltBuff, Units,"_HostExpectedDataCounts.csv"))  
  #
  ## Identify and combine classes that have expected counts below 5 (c.f. Bryers et al. 1984, Manly et al. 2002)
  #
  CombineDropClasses.df <- ActualDatTab.df[ which(ActualDatTab.df$ExpectedCount<5),]
  print(CombineDropClasses.df)
  #str(DropClasses.df)
  # Save Dropped Classes
  write.table(CombineDropClasses.df, sep= ",", col.names=NA, file=paste0(Species, "_", Season, "CombinedHostsforMonarchHostSelectivity.csv"))
  #
  # Save above list of dropped classes to text file
  out1 <- paste0("List of Hosts combined that had expected count less than 5")
  cat(out1, file=paste0(Species, "_", Season, "ListCombinedHostsforMonarchHostSelectivity.txt"), sep="\n", append=TRUE)
  out2 <- capture.output(print(CombineDropClasses.df))
  cat(out2, file=paste0(Species, "_", Season, "ListCombinedHostsforMonarchHostSelectivity.txt"), sep="\n", append=TRUE)
  #
  CombineDropClasses <- CombineDropClasses.df$Host
  ## Examine data to determine which similar classes to combine and change names for combined classes
  ## Combine Hosts for ActualDatTab.df using new HostsN
  #
  if(i==1) {
    ActualDatTab.df$HostN <- ActualDatTab.df$Host
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. syriaca"] <- "A. speciosa Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. speciosa"] <- "A. speciosa Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. sullivantii"] <- "A. speciosa Grp"
    #ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. arenaria"] <- "A. arenaria Grp"
    #ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. brachystephana"] <- "A. arenaria Grp"
    #ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. subverticillata"] <- "A. arenaria Grp"
    #ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. incarnata"] <- "A. linearis Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. purpurascens"] <- "A. quadrifolia Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. quadrifolia"] <- "A. quadrifolia Grp"
    #ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. linearis"] <- "A. linearis Grp"
    #
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. arenaria"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. hirtella"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. texana"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="Cynanchum laeve"] <- "Other Grp"
  } else {
    ActualDatTab.df$HostN <- ActualDatTab.df$Host
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. brachystephana"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. engelmanniana"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. incarnata"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. linearis"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. perennis"] <- "Other Grp"    
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. purpurascens"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. subverticillata"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. syriaca"] <- "Other Grp"    
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="A. texana"] <- "Other Grp"
    ActualDatTab.df$HostN[ActualDatTab.df$HostN=="Matelea reticulata"] <- "Other Grp"
  }
  write.table(ActualDatTab.df, sep = ",", col.names=NA, file=paste0(Species, "_", Season, "_Thinned", SpatFiltBuff, Units,"_OriginalHostDataCounts.csv"))
  #
  ActualDatTab.df2 <- ActualDatTab.df
  ActualDatTab.df2$Host <- ActualDatTab.df2$HostN
  ActualDatTab.df2 <- ActualDatTab.df2[, !(names(ActualDatTab.df2) %in% c("HostN"))]
  # Resummarize by Host
  Columns <- colnames(ActualDatTab.df2)[2:ncol(ActualDatTab.df2)]
  ActualDatTab.df3 <- aggregate(. ~ Host, data=ActualDatTab.df2, FUN=function(x) c(Sum = sum(x)))
  #
  # Save data
  write.table(ActualDatTab.df3, sep = ",", col.names=NA, file=paste0(Species, "_", Season, "_Thinned", SpatFiltBuff, Units,"_GroupedHostDataCounts.csv"))
  #

  ###############################################################################################################
  ## Calculate Bonferroni Confidence Intervals for whether RandomPercAvail falls within BCI for PercUsed
  ## based on raw data and bias adjustments for human populations
  ## following Bryers et al. 1984 and Manly et al. 2002
  ###############################################################################################################
  # Read back in data
  ActualDatTab.df3 <- data.frame(read.csv(paste0(Species, "_", Season, "_Thinned", SpatFiltBuff, Units,"_GroupedHostDataCounts.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
  ##
  ## Delete row names column
  ActualDatTab.df3 <- ActualDatTab.df3[,-1]
  ##
  sum(ActualDatTab.df3$ExpectedCount)
  sum(ActualDatTab.df3$PopDenNHBRWtN)
  sum(ActualDatTab.df3$PopDenNHBRWtNPercAvail)
  NumberComparisons <- nrow(ActualDatTab.df3)
  #
  ## Loop through the test for the raw and six bias adjustments for percent land cover used
  ## and calculate electivity indices
  UsedColumnList <- seq(2, 3, 1)
  UsedCategoryList <- colnames(ActualDatTab.df3)[UsedColumnList]
  PercUsedColumnList <- seq(4,5,1)
  PercUsedCategoryList <- colnames(ActualDatTab.df3)[PercUsedColumnList]
  count=0
  for(UsedCol in UsedColumnList) {
    #UsedCol=2
    PercUsedCol <- UsedCol + 2
    count=count+1
    ChiSqGFT_Results <- chisq.test(ActualDatTab.df3[,UsedCol], p=ActualDatTab.df3$PopDenNHBRWtNPercAvail)
    #str(ChiSqGFT_Results)
    #ChiSqGFT_Results$expected
    # Save above test results
    out1 <- paste0("Chi-Square Goodness of Fit Test for Differeces in Expected (Random) Versus ", UsedCategoryList[count]," Utilization of Milkweed Hosts")
    cat(out1, file=paste0(Species, "_", Season, "2011_2020MilkweedHostChiSquareGoodnessFitTestfor", UsedCategoryList[count],".txt"), sep="\n")
    out2 <- capture.output(print(ChiSqGFT_Results))
    cat(out2, file=paste0(Species, "_", Season, "2011_2020MilkweedHostChiSquareGoodnessFitTestfor", UsedCategoryList[count],".txt"), sep="\n", append=TRUE)
    #
    # Assuming above ChiSquare Test is significant, find bonferoni critical z-value of 0.05 with NumberComparisons categories for two-tailed 95% CI
    pc <- 0.05/(NumberComparisons)
    za <- qnorm(1-(pc/2))
    ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCI")] <- za*sqrt((ActualDatTab.df3[,PercUsedCol]*(1-ActualDatTab.df3[,PercUsedCol]))/sum(ActualDatTab.df3$ExpectedCount))
    ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"LowerBCI")] <- ActualDatTab.df3[,PercUsedCol] - ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCI")]
    ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"UpperBCI")] <- ActualDatTab.df3[,PercUsedCol] + ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCI")]
    # If RandomPercAvail outside BCI, identify as Significant difference
    if(ChiSqGFT_Results$p.value<0.05) {
      ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCISig")] <- ifelse(ActualDatTab.df3$PopDenNHBRWtNPercAvail <= ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"UpperBCI")] & ActualDatTab.df3$PopDenNHBRWtNPercAvail >=  ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"LowerBCI")],  "NonSig",  "Significant")
      } else {
      ActualDatTab.df3[, paste0(PercUsedCategoryList[count],"BCISig")] <- "NonSig"
      #
    }
    ActualDatTab.df3[, paste0("ForateRat", PercUsedCategoryList[count])] <- ivlev_forage(ActualDatTab.df3[,PercUsedCol], ActualDatTab.df3$PopDenNHBRWtNPercAvail, log10=FALSE)
    ActualDatTab.df3[, paste0("VS_Elect", PercUsedCategoryList[count])] <- vs_electivity(ActualDatTab.df3[,PercUsedCol], ActualDatTab.df3$PopDenNHBRWtNPercAvail)
    ActualDatTab.df3[, paste0("VS_Elect", PercUsedCategoryList[count], "_10K")] <- Round2(ActualDatTab.df3[, paste0("VS_Elect", PercUsedCategoryList[count])]*10000,0)
  }
  ## Calculated Normalized VS_Elect
  VS_ElectMin <- min(ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsed)
  VS_ElectMax <- max(ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsed)
  # Calculate Normalized PopDen, PopDenN
  ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedN_10K <- 10000*((ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsed - VS_ElectMin)/(VS_ElectMax - VS_ElectMin))
  head(ActualDatTab.df3)
  # Calculate Percent of Total Normalized VS_Elect
  ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedNP <- ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedN_10K/sum(ActualDatTab.df3$VS_ElectPopDenNHBRWtNPercUsedN_10K)
  #
  ## Save data
  write.table(ActualDatTab.df3, sep = ",", col.names=NA, file=paste0(Species, "_", Season, "_MilkweedGroupHostActualPercentUsedStatsIndices.csv"))
}