#
Species <- "Asclepias_latifolia"
InDirect <- "K:/Monarch/"
#
# Read csv file of weekly histogram
setwd(InDirect)
SpeciesData.df <- data.frame(read.csv(paste0(Species, "_histogrampoints.csv"), stringsAsFactors=FALSE))
head(SpeciesData.df)
plot(SpeciesData.df$Count)
#
## define the integrated function for early season
integrand <- function(x) {(6.943701847130800E-12*6.355663334470514E-01^x)*(x^1.229030569676149E+01)}  # Hoerl Model
# add plot of curve
SpeciesData.df$hoerlcurve <- sapply(SpeciesData.df$WeekNum, integrand)## integrate the function from 0 to infinity
lines(SpeciesData.df$hoerlcurve, col="blue")
# integrate function
integrate(integrand, lower = 0, upper = 31)
integrate(integrand, lower = 11, upper = 28)
#
## define the integrated function for late season
integrand <- function(x) {(1.007291622917991E+00 + -2.027875328233310E-02*x)/(1 + -4.979268107935871E-02*x + 6.334700520665867E-04*x^2)}  # Rational Model
# add plot of curve
SpeciesData.df$rationalcurve2 <- sapply(SpeciesData.df$WeekNum, integrand)## integrate the function from 0 to infinity
lines(SpeciesData.df$rationalcurve2, col="red")
## integrate the function from 0 to infinity
integrate(integrand, lower = 32, upper = 52)
integrate(integrand, lower = 32, upper = 47)
#
####################################################
## Chi-square Goodness of Fit Tests for early and late season proportions
#
Total <- sum(SpeciesData.df$Count)
SpringCount <- SpeciesData.df[ which(SpeciesData.df$WeekNum>=11 & SpeciesData.df$WeekNum <= 28),]
SpringTotal <- sum(SpringCount$Count)
FallCount <- SpeciesData.df[ which(SpeciesData.df$WeekNum>=32 & SpeciesData.df$WeekNum <= 47),]
FallTotal <- sum(FallCount$Count)
ExpectedCount <- (SpringTotal + FallTotal)/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 Spring vs Fall Milkweed Frequency")
cat(out1, file=paste0(Species, "ChisquareGFTSpringVsFall.txt"), sep="\n", append=TRUE)
out2 <- paste0("Spring Count: ", SpringTotal)
cat(out2, file=paste0(Species, "ChisquareGFTSpringVsFall.txt"), sep="\n", append=TRUE)
out3 <- paste0("Fall Count: ", FallTotal)
cat(out3, file=paste0(Species, "ChisquareGFTSpringVsFall.txt"), sep="\n", append=TRUE)
out4 <- paste0("Expected Count: ", ExpectedCount)
cat(out4, file=paste0(Species, "ChisquareGFTSpringVsFall.txt"), sep="\n", append=TRUE)
out5 <- capture.output(print(ChiSqGFTres))
cat(out5, file=paste0(Species, "ChisquareGFTSpringVsFall.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$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")
## Save data
write.table(SeasonPercent.df, sep = ",", col.names=NA, file=paste0(Species, "_SeasonPercentStats.csv"))


## Save data
write.table(SeasonPercent.df, sep = ",", col.names=NA, file=paste0(Species, "_SeasonPercentStats.csv"))
