## 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(plotrix)
#
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"
PopGridName <- paste0("sc2010pop1k")
# Read in data
PhenologyDat.df <- data.frame(read.csv(paste0(Species, "_ImmaturePhenologyUnthinnedWeekYearSCStatesRandom.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(PhenologyDat.df)
head(PhenologyDat.df)
tail(PhenologyDat.df)
nrow(PhenologyDat.df)
unique(PhenologyDat.df$Phenological.Stage)
#str(PhenologyDat.df)
# Only keep data after 2010
PhenologyDat.df1 <- PhenologyDat.df[which(PhenologyDat.df$Year>2010),]
nrow(PhenologyDat.df1)
# Drop exotic Apocynaceae
Drop <- c("Gomphocarpus", "curassavica")
PhenologyDat.df1 <- PhenologyDat.df1[!grepl(paste(Drop, collapse= "|"), PhenologyDat.df1$Host), ]
nrow(PhenologyDat.df1)
head(PhenologyDat.df1)
unique(PhenologyDat.df1$Host)
# Drop blanks for Host
PhenologyDat.df1 <- PhenologyDat.df1[!(PhenologyDat.df1$Host=="."), ]
nrow(PhenologyDat.df1)
#
# Delete records with comment of "Ornamental"
PhenologyDat.df1 <- PhenologyDat.df1[!grepl("Ornamental", PhenologyDat.df1$Comments), , drop = FALSE] 
head(PhenologyDat.df1)
nrow(PhenologyDat.df1)
#
# Delete records from David Berman's late season study
PhenologyDat.df1 <- PhenologyDat.df1[!grepl("Berman", PhenologyDat.df1$Comments), , drop = FALSE] 
head(PhenologyDat.df1)
nrow(PhenologyDat.df1)
#
# Randomize order of data
set.seed(42)
rows <- sample(nrow(PhenologyDat.df1))
PhenologyDat.df1 <- PhenologyDat.df1[rows, ]
#
# Thin data by 30 km within same week and year
YearList <- unique(PhenologyDat.df1$Year)
WeekList <- unique(PhenologyDat.df1$WeekNum)
PhenologyDatList1 <- list()
PhenologyDatList2 <- list()
SpeciesList <- unique(PhenologyDat.df1$Host)
length(SpeciesList)
#
Count2 <- 0
#
for(SpeciesName in SpeciesList) {
  #SpeciesName <- "A. asperula"
  PhenologyDat.df1a <- PhenologyDat.df1[which(PhenologyDat.df1$Host==SpeciesName),]
  Count1 <- 0
  for(Year in YearList) {
    #Year <- 2020
    PhenologyDat.df2 <-  PhenologyDat.df1a[which(PhenologyDat.df1a$Year==Year),]
    for(Week in WeekList) {
      #Week <- 15
      PhenologyDat.df3 <-  PhenologyDat.df2[which(PhenologyDat.df2$WeekNum==Week),]
      head(PhenologyDat.df3)
      nrow(PhenologyDat.df3)
      # Randomize order of rows
      set.seed(42)
      rows <- sample(nrow(PhenologyDat.df3))
      PhenologyDat.df3 <- PhenologyDat.df3[rows, ]
      # Remove duplicate values for Lat/Lon
      PhenologyDat.df3 <- PhenologyDat.df3[!duplicated(PhenologyDat.df3[c("lat", "lon")]),]
      if(nrow(PhenologyDat.df3)< 3) {
        if(nrow(PhenologyDat.df3)==2) {
          # calculate distance between the two points
          PointDist <- distm(c(PhenologyDat.df3$lon[1], PhenologyDat.df3$lat[1]), c(PhenologyDat.df3$lon[2], PhenologyDat.df3$lat[2]), fun = distHaversine)
          # Delete second record if closer than 30 km (30,000 m)
          if(PointDist < 30000) {
             PhenologyDat.df4 <- PhenologyDat.df3[1,]
          } else {
             PhenologyDat.df4 <- PhenologyDat.df3
          }
        } else {
          PhenologyDat.df4 <- PhenologyDat.df3
        }
      } else {
      # Spatially thin to 30 km
      # Select only lat lon and id columns
      pres.df <- PhenologyDat.df3[,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 30 km
      presthin.df <- as.data.frame(thin(pres.df, lat.col="lat", long.col="lon", spec.col = "Species", thin.par=30, reps=1, write.file=FALSE, locs.thinned.list.return = TRUE))
      colnames(presthin.df) <- c("lon", "lat")
      # Merge data with original roadkill by lat/long
      PhenologyDat.df4 <- merge(presthin.df, PhenologyDat.df3, by=c("lon", "lat"))
      head(PhenologyDat.df4)
      nrow(PhenologyDat.df4)
      PhenologyDat.df4 <- PhenologyDat.df4[,c(3:4, 1:2, 5:ncol(PhenologyDat.df4))]
      }
      # Save data to list
      Count1 <- Count1 + 1
      PhenologyDatList1[[Count1]] <- PhenologyDat.df4
    }
  }
  Count2 <- Count2 + 1
  PhenologyDatList2[[Count2]] <- do.call(rbind, PhenologyDatList1)
}  
##########################
PhenologyDatAll.df <- do.call(rbind, PhenologyDatList2)
colnames(PhenologyDatAll.df)
head(PhenologyDatAll.df)
tail(PhenologyDatAll.df)
nrow(PhenologyDatAll.df)
write.table(PhenologyDatAll.df, file=paste0(Species, "_PhenologyThinnedWeekYear30kmSC2011_20_AllNativeMilkweeds.csv"), sep=",", row.names=FALSE)
############################################
###
##################################################################################################
### Histogram Bin Ratio Weighting (HBRWt) for correcting sample bias for thinned phenology 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
################################################################################################### 
#
# Read back in spatiotemporally thinned data
PhenologyDatAll.df <- data.frame(read.csv(paste0(Species, "_PhenologyThinnedWeekYear30kmSC2011_20_AllNativeMilkweeds.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(PhenologyDatAll.df)
head(PhenologyDatAll.df)
tail(PhenologyDatAll.df)
nrow(PhenologyDatAll.df)
#
## Extract population densities for data
# Convert data to spatial points data frame
xy <- PhenologyDatAll.df[,c("lon","lat")]
PhenologyDatAll.spdf <- SpatialPointsDataFrame(coords = xy, data = PhenologyDatAll.df, proj4string = crs(StudyAreaShp))
head(PhenologyDatAll.spdf)
nrow(PhenologyDatAll.spdf)
#
# Add population data
## Read in population density grid
PopDenRaster <- raster(PopGridName)
#plot(PopDenRaster, add=TRUE, axes=TRUE)
PhenologyDatAll.spdf1 <- raster::extract(PopDenRaster, PhenologyDatAll.spdf, sp=TRUE)
head(PhenologyDatAll.spdf1)
nrow(PhenologyDatAll.spdf1)
# Convert back to data frame
PhenologyDatAll.df2 <- as.data.frame(PhenologyDatAll.spdf1, stringsAsFactors=FALSE)
head(PhenologyDatAll.df2)
# Delete last two columns
PhenologyDatAll.df2 <- PhenologyDatAll.df2[,-c((ncol(PhenologyDatAll.df2)-1):(ncol(PhenologyDatAll.df2)))]
nrow(PhenologyDatAll.df2)
# Rename last column to "PopDen"
colnames(PhenologyDatAll.df2)[ncol(PhenologyDatAll.df2)] <- "PopDen"
nrow(PhenologyDatAll.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, ])
}
PhenologyDatAll.df2 <- completeFun(PhenologyDatAll.df2, "PopDen")
nrow(PhenologyDatAll.df2)
head(PhenologyDatAll.df2)
#
PhenologyDatAll.df2$RawCount <- 1
head(PhenologyDatAll.df2)
nrow(PhenologyDatAll.df2)
#
### Read in thinned target group points and random points shapefiles and convert to data frames
### Obtained from program 'Danus__plexippus_Immature_BHRWeight_HostRecords.R'
## 
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
#
#
##################################################################
##  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, PhenologyDatAll.df2$PopDen)
PopDenMax <- max(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, PhenologyDatAll.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)
PhenologyDatAll.df2$PopDenN <- (PhenologyDatAll.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(PhenologyDatAll.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
PhenologyDatAll.df2$PopDenNBinCode <- .bincode(PhenologyDatAll.df2$PopDenN, seq(0,1,by=0.001), include.lowest=TRUE)
head(PhenologyDatAll.df2)
# Retrieve HBRWt for PopDenN
PhenologyDatAll.df2$PopDenNHBRWt <- PopDenHBRWt.df[PhenologyDatAll.df2$PopDenNBinCode,7]
head(PhenologyDatAll.df2)
PopDenHBRWt.df[1:20,]
############
# Normalize all weighting factors to sum to total raw count for calculating expected counts
TotRawCount <- sum(PhenologyDatAll.df2$RawCount)
TotPopDenNHBRWt <- sum(PhenologyDatAll.df2$PopDenNHBRWt)
# Normalize weights to sum to TotRawCount
PhenologyDatAll.df2$PopDenNHBRWtN <- PhenologyDatAll.df2$PopDenNHBRWt/(TotPopDenNHBRWt/TotRawCount)
sum(PhenologyDatAll.df2$PopDenNHBRWtN)
head(PhenologyDatAll.df2)
# Write data
write.table(PhenologyDatAll.df2, file=paste0(Species, "_PhenologyThinnedWeekYear30kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt.csv"), sep=",", row.names=FALSE)
##
###########################
## Obtain frequency of plants by species using PopDenNHBRWtN
EvaluationStats <- PhenologyDatAll.df2[,c("Host", "PopDenNHBRWtN")]
head(EvaluationStats)
OutName <- paste0(Species, "_HostFrequencyforThinnedWeekYear30kmSC2011_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)
##################################################################################
### Read in seasonal species frequency by PopDenHBRWt and anlayze for differences
##################################################################################
##
#
PhenologyDatFreq.df <- data.frame(read.csv(paste0(Species, "_PhenologyThinnedWeekYear30kmSC2011_20_AllNativeMilkweeds_PopDenHBRWt_WeightedHistogramSeasonalFrequencies.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
#
head(PhenologyDatFreq.df)
###############################################################################################################
## 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
###############################################################################################################
##
####################################################
## Chi-square Goodness of Fit Tests for early and late season proportions
#
SpeciesList <- colnames(PhenologyDatFreq.df)[2:ncol(PhenologyDatFreq.df)]
SpeciesFreq <- list()
Count <- 0
#
for(SpeciesF in SpeciesList) {
  #SpeciesF <- "A.viridis"
  SpringTotal <- PhenologyDatFreq.df[1, SpeciesF]    
  FallTotal <- PhenologyDatFreq.df[2, SpeciesF]  
  ExpectedCount <- (PhenologyDatFreq.df[3, SpeciesF])/2
  #
  SeasonCounts <- c(SpringTotal, FallTotal)
  ChiSqGFTres <- chisq.test(SeasonCounts, p = c(1/2,1/2))
  ChiSqGFTres$expected
  #
  # Save above test output
  out1 <- paste0("Chi-square Goodneess of Fit test for ", SpeciesF, " Spring vs Fall Monarch Larval Frequency")
  cat(out1, file=paste0(Species, "ChisquareGFT_", SpeciesF, "_SpringVsFall.txt"), sep="\n", append=TRUE)
  out2 <- paste0("Spring Count: ", SpringTotal)
  cat(out2, file=paste0(Species, "ChisquareGFT_", SpeciesF, "_SpringVsFall.txt"), sep="\n", append=TRUE)
  out3 <- paste0("Fall Count: ", FallTotal)
  cat(out3, file=paste0(Species, "ChisquareGFT_", SpeciesF, "_SpringVsFall.txt"), sep="\n", append=TRUE)
  out4 <- paste0("Expected Count: ", ExpectedCount)
  cat(out4, file=paste0(Species, "ChisquareGFT_", SpeciesF, "_SpringVsFall.txt"), sep="\n", append=TRUE)
  out5 <- capture.output(print(ChiSqGFTres))
  cat(out5, file=paste0(Species, "ChisquareGFT_", SpeciesF, "_SpringVsFall.txt"), sep="\n", append=TRUE)
  #
  # Find bonferoni critical z-value of 0.05 with NumberComparisons categories for two-tailed 95% CI
  SeasonPercent.df <- data.frame(c("MidMar-MidJul", "MidAug-LateNov"))
  colnames(SeasonPercent.df)[1] <- "Season"
  SeasonPercent.df$Species <- SpeciesF
  SeasonPercent.df$SeasonCounts <- SeasonCounts
  SeasonPercent.df$SeasonPercent <-SeasonCounts/sum(SeasonCounts)
  SeasonPercent.df$RandomPercent <- c(rep(0.5,2))
  SeasonPercent.df$ExpectedCount <- c(rep(ExpectedCount,2))
  #
  pc <- 0.05/2
  za <- qnorm(1-(pc/2))
  SeasonPercent.df$SeasonPercentBCI <- za*sqrt((SeasonPercent.df$SeasonPercent*(1-SeasonPercent.df$SeasonPercent))/sum(SeasonPercent.df$ExpectedCount))
  SeasonPercent.df$PercUsedLowerBCI <- SeasonPercent.df$SeasonPercent - SeasonPercent.df$SeasonPercentBCI
  SeasonPercent.df$PercUsedUpperBCI <- SeasonPercent.df$SeasonPercent + SeasonPercent.df$SeasonPercentBCI
  # If RandomPercAvail outside BCI, place asterisk under BCISig
  SeasonPercent.df$BCISig <- ifelse(SeasonPercent.df$RandomPercent <= SeasonPercent.df$PercUsedUpperBCI & SeasonPercent.df$RandomPercent >= SeasonPercent.df$PercUsedLowerBCI,  "NonSig",  "Significant")
  Count <- Count + 1
  SpeciesFreq[[Count]] <- SeasonPercent.df
}
SeasonPercentAll.df <- do.call(rbind, SpeciesFreq)
## Save Data
write.table(SeasonPercentAll.df, sep = ",", col.names=NA, file=paste0(Species, "_HostSeasonFrequencyPercentStats.csv"))
#
######### Following code not used, but made pivot tables in excel to build stacked histograms by PopDenHBRWt per stage and species
## Create weighted histograms for plotting in excel
PhenologyHistogram <- weighted.hist(PhenologyDatAll.df2$WeekNum, PhenologyDatAll.df2$PopDenNHBRWtN, breaks=1:52)
#str(PhenologyHistogram)
PhenologyHistogram.df <- cbind(lapply(PhenologyHistogram$mids, function(x) x + 0.5), PhenologyHistogram$counts)
head(PhenologyHistogram.df)
colnames(PhenologyHistogram.df) <- c("WeekNum", "PopDenNHBRWtN_Counts")
# Write data
write.table(PhenologyHistogram.df, file=paste0(Species, "_WeightedPhenologyHistograms_PopDenHBRWt_ThinnedWeekYear30kmSC2011_20_AllNativeMilkweeds.csv"), sep=",", row.names=FALSE)

