####Code for developing and mixing distributions of overwintering monarch butterfly densities and associated implications for population size
###W. E. Thogmartin, wthogmartin@usgs.gov

# This software has been approved for release by the U.S. Geological Survey (USGS). Although the software has been subjected to rigorous review, the USGS reserves the right to update the software as needed pursuant to further analysis and review. No warranty, expressed or implied, is made by the USGS or the U.S. Government as to the functionality of the software and related material nor shall the fact of release constitute any such warranty. Furthermore, the software is released on condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from its authorized or unauthorized use.

library(plotrix)
library(rriskDistributions)
library(TeachingDemos)
library(MASS)
library(outliers)
library(Unicode)
library(reshape2)
library(evd)
library(fitdistrplus)
library(lattice)
library(ggplot2)
library(AICcmodavg)
library(effects)
library(ggplot2)
library(matrixStats)


#############Density estimates from previously published studies

uplimP=100    #Upper 97.5 limit reported for the Petersen estimate
lowlimP=40.6  #Lower 2.5% limit reported for the Petersen estimate
meanP=60.9    #Mean reported for the Petersen estimate
uplimP-meanP  #Identifying extent of asymmetry in the confidence interval
meanP-lowlimP

uplimJS=68.7  #Upper 97.5 limit reported for the Jolly-Seber estimate
lowlimJS=21.1 #Lower 2.5 limit reported for the Jolly-Seber estimate
meanJS=33.8   #Mean reported for the Jolly-Seber estimate
uplimJS-meanJS #Identifying extent of asymmetry in the confidence interval
meanJS-lowlimJS

uplimDec=9.6    #Upper 97.5 limit reported for the Petersen estimate
lowlimDec=5.6  #Lower 2.5% limit reported for the Petersen estimate
meanDec=6.9    #Mean reported for the Petersen estimate
uplimDec-meanDec  #Identifying extent of asymmetry in the confidence interval
meanDec-lowlimDec

######Petersen January
log(meanP)
lmeanP <- log(meanP)

(log(meanP)-log(lowlimP))/2.11  #0.1921636 for t = 2.11; 0.20687 for z = 1.96
(log(uplimP)-log(meanP))/2.11  #0.2350412; 0.2530291 for z = 1.96

(((log(meanP)-log(lowlimP))/2.11) + ((log(uplimP)-log(meanP))/2.11))/2   #0.2136024; 0.2299495
SDP <- (((log(meanP)-log(lowlimP))/2.11) + ((log(uplimP)-log(meanP))/2.11))/2   #0.2136024; 0.2299495


thetaP <- rlnorm(1000000, lmeanP, SDP)  #Re-constructing a lognormal distribution for the Petersen estimate
#summary(thetaP)
quantile(thetaP, probs=c(0.025, 0.5, 0.975))



#######Jolly-Seber
log(meanJS)
lmeanJS <- log(meanJS)

(log(meanJS)-log(lowlimJS))/2.11 
(log(uplimJS)-log(meanJS))/2.11  

(((log(meanJS)-log(lowlimJS))/2.11) + ((log(uplimJS)-log(meanJS))/2.11))/2  
SDJS <- (((log(meanJS)-log(lowlimJS))/2.11) + ((log(uplimJS)-log(meanJS))/2.11))/2  


thetaJS <- rlnorm(1000000, lmeanJS, SDJS)
summary(thetaJS)
quantile(thetaJS, probs=c(0.025, 0.5, 0.975))


######Petersen December
log(meanDec)
lmeanDec <- log(meanDec)

(log(meanDec)-log(lowlimDec))/2.11  #0.1921636 for t = 2.11; 0.20687 for z = 1.96
(log(uplimDec)-log(meanDec))/2.11  #0.2350412; 0.2530291 for z = 1.96

(((log(meanDec)-log(lowlimDec))/2.11) + ((log(uplimP)-log(meanP))/2.11))/2   #0.2136024; 0.2299495
SDDec <- (((log(meanDec)-log(lowlimDec))/2.11) + ((log(uplimP)-log(meanP))/2.11))/2   #0.2136024; 0.2299495


thetaDec <- rlnorm(1000000, lmeanDec, SDDec)  #Re-constructing a lognormal distribution for the Petersen estimate
#summary(thetaP)
quantile(thetaDec, probs=c(0.025, 0.5, 0.975))

exp(SDDec)


############Branch
#n = 17, sample size of Branches in Calvert (2004)
sdP <- sqrt(17)*(uplimP - lowlimP)/4.22     #Re-constructing the standard deviation from the confidence interval
sdJS <- sqrt(17)*(uplimJS - lowlimJS)/4.22

sdP
log(sdP)
sdJS

branchCIupP <- 10.28 + 2.11*sqrt(sdP/17)    #Constructing an interval for Branch, from sd for Petersen
branchCIlowP <- 10.28 - 2.11*sqrt(sdP/17)

branchCIupJS <- 10.28 + 2.11*sqrt(sdJS/17)    #Constructing an interval for Branch, from sd for Jolly-Seber
branchCIlowJS <- 10.28 - 2.11*sqrt(sdJS/17)


meanB = 10.3
lmeanB = log(meanB)  #logged

sdPscaled <- (((log(meanB)-log(branchCIlowP))/2.11) + ((log(branchCIupP)-log(meanB))/2.11))/2  
exp(sdPscaled)

###E[x]
exp(lmeanB)


#SDB = (SDP + SDJS)/2   calculated an average of SDs from the CMR results
thetaB <- rlnorm(1000000, lmeanB, sdPscaled) #note use of sdPscaled
summary(thetaB)
quantile(thetaB, probs=c(0.025, 0.5, 0.975))


##########Alternative Branch Variance
#Page 124 of Calvert (2004) identifies 5 steps in the calculation of butterfly abundance on crowns, and 3 on trunks
#Thus, according to Machtans and Thogmartin (2014), the range of variation in this Fermi approximation is on the order of

#if we assume each step is correct within a factor of 2
###Tree Crowns
uppercrown <- 2^sqrt(5)
lowercrown <- 1/(2^sqrt(5))

meancrown <- 24961256
lowerlimcrown <- meancrown*lowercrown
upperlimcrown <- meancrown*uppercrown

###Tree Trunks
uppertrunk <- 2^sqrt(3)
lowertrunk <- 1/(2^sqrt(3))

meantrunk <- 2688947
lowerlimtrunk <- meantrunk*lowertrunk
upperlimtrunk <- meantrunk*uppertrunk

meanbranch <- meancrown + meantrunk
lowerbranch <- lowerlimcrown + lowerlimtrunk
upperbranch <- upperlimcrown + upperlimtrunk

meanbranchdensity <- meanbranch/2.69
lowerbranchdensity <- lowerbranch/2.69
upperbranchdensity <- upperbranch/2.69

meanbranchdens <- meanbranchdensity/1000000
lowerbranchdens <- lowerbranchdensity/1000000
upperbranchdens <- upperbranchdensity/1000000

log(meanbranchdens)
lmeanbranch <- log(meanbranchdens)

log(lowerbranchdens)
llowerbranch <- log(lowerbranchdens)

log(upperbranchdens)
lupperbranch <- log(upperbranchdens)

(lmeanbranch - llowerbranch)/2.11  
(lupperbranch - lmeanbranch)/2.11  

(((lmeanbranch - llowerbranch)/2.11) + ((lupperbranch - lmeanbranch)/2.11))/2   
SDbranch <- (((lmeanbranch - llowerbranch)/2.11) + ((lupperbranch - lmeanbranch)/2.11))/2


thetabranch <- rlnorm(1000000, lmeanbranch, SDbranch)  #Re-constructing a lognormal distribution for the Petersen estimate
#summary(thetaP)
quantile(thetabranch, probs=c(0.025, 0.5, 0.975))


exp(SDbranch)





################Plot

x <- seq(0,150, length=1000000)

hist(thetaP)
dP <- density(thetaP)
plot(dP, xlim=c(0,150), ylim=c(0,0.045), lwd=2, main="", sub="")

dJS <- density(thetaJS)
par(new=TRUE)
plot(dJS, xlim=c(0,150), ylim=c(0,0.045), lwd=2, main="", sub="")

#polygon(dP$x, pmin(dP$y, dJS$y), col="gray")
abline(v=60.9, lwd=2, lty=4)
abline(v=33.8, lwd=2, lty=3)
abline(v=10.3, lwd=2, lty=2)

dB <- density(thetabranch)
par(new=TRUE)
plot(dB, xlim=c(0,150), ylim=c(0,0.18), main="", lwd=2, sub="")


#




plot(0, xlim=c(0,110), ylim=c(0,0.35), main=NA, bty="n", pch="", xlab="Population Density (millions/ha)", ylab="Density")
lines(density(thetaP), xlim=c(0,110), ylim=c(0,0.35), main=NA, bty="n", lwd=2, lty=4)
#par(new=TRUE)
lines(density(thetaJS), xlim=c(0,110), ylim=c(0,0.35), main=NA, lwd=2, lty=3, ylab="", yaxt="n", xaxt="n", bty="n")
#par(new=TRUE)
lines(density(thetabranch), xlim=c(0,110), ylim=c(0,0.35), main=NA, lwd=2, lty=2, ylab="", yaxt="n", xaxt="n", bty="n")
lines(density(thetaDec), xlim=c(0,110), ylim=c(0,0.35), main=NA, lwd=2, lty=2, ylab="", yaxt="n", xaxt="n", bty="n")
ablineclip(v=60.9,y1=0.01, y2=0.05, lwd=2, lty=1, col="dark gray")
ablineclip(v=33.8, y1=0.01, y2=0.05, lwd=2, lty=1, col="dark gray")
ablineclip(v=10.3, y1=0.01, y2=0.05, lwd=2, lty=1, col="dark gray")
ablineclip(v=6.9, y1=0.01, y2=0.05, lwd=2, lty=1, col="dark gray")
text(18,0.15, "Branch", cex=1.2)
text(46,0.05, "Jolly-Seber", cex=1.2)
text(78,0.044, "Petersen (Jan)", cex=1.2)
text(26,0.074, "Petersen (Dec)", cex=1.2)





#####December

plot(0, xlim=c(0,20), ylim=c(0,0.35), main=NA, bty="n", pch="", xlab="Population Density (millions/ha)", ylab="Density")
lines(density(thetaDec), xlim=c(0,20), ylim=c(0,0.35), main=NA, lwd=2, lty=5, ylab="", yaxt="n", xaxt="n", bty="n")
ablineclip(v=6.9, y1=0.01, y2=0.31, lwd=2, lty=1, col="dark gray")
text(13,0.24, "Petersen (Dec)", cex=1.2)



######January

plot(0, xlim=c(0,110), ylim=c(0,0.08), main=NA, bty="n", pch="", xlab="Population Density (millions/ha)", ylab="Density")
lines(density(thetabranch), xlim=c(0,110), ylim=c(0,0.08), main=NA, lwd=2, lty=2, ylab="", yaxt="n", xaxt="n", bty="n")
lines(density(thetaP), xlim=c(0,110), ylim=c(0,0.35), main=NA, bty="n", lwd=2, lty=4)
lines(density(thetaJS), xlim=c(0,110), ylim=c(0,0.35), main=NA, lwd=2, lty=3, ylab="", yaxt="n", xaxt="n", bty="n")
ablineclip(v=10.3, y1=0.001, y2=0.04, lwd=2, lty=1, col="dark gray")
ablineclip(v=60.9,y1=0.001, y2=0.04, lwd=2, lty=1, col="dark gray")
ablineclip(v=33.8, y1=0.001, y2=0.04, lwd=2, lty=1, col="dark gray")
text(18,0.07, "Branch", cex=1.2)
text(46,0.05, "Jolly-Seber", cex=1.2)
text(78,0.044, "Petersen (Jan)", cex=1.2)









dev.off()

########Brower as a Lognormal##########

# library(rriskDistributions)
# library(TeachingDemos)
# library(MASS)
# library(outliers)

conejosmd <- c(1,
               104,
               287,
               54,
               11,
               86,
               12,
               81,
               1097,
               43,
               2354,
               195,
               77,
               21,
               189,
               26,
               125,
               18,
               136,
               14,
               52,
               148,
               132,
               72,
               763,
               462,
               75,
               589,
               1199)

zapateromd <- c(109,
                158,
                80,
                9,
                118,
                41,
                45,
                52,
                372,
                158,
                17,
                26,
                0,
                0,
                0,
                22,
                209,
                238,
                55,
                166,
                152,
                275,
                45,
                121,
                92,
                83,
                266,
                100,
                40)

#library(Unicode)
summary(zapateromd)
var(zapateromd)
sd(zapateromd)
summary(conejosmd)
var(conejosmd)
sd(conejosmd)
truehist(zapateromd, xlim=c(0,2500))
truehist(conejosmd)
plot(conejosmd, ylab=expression(paste('Count per 0.2 \U00D7 0.2 m'^'2')))
points(zapateromd, pch=16, add=TRUE)





#library(reshape2)
zapcon <- melt(data.frame(zapateromd,conejosmd))
colnames(zapcon) <- c("Group", "Count")
zapcon$Group <- as.character(zapcon$Group)
zapcon$Group <- replace(zapcon$Group, zapcon$Group=="zapateromd", "Zapatero")
zapcon$Group <- replace(zapcon$Group, zapcon$Group=="conejosmd", "Conejos")
print(zapcon)
tapply(zapcon$Count, zapcon$Group, summary)



boxplot(Count ~ Group, data=zapcon, outpch = NA, ylab=expression(paste('Count per 0.2 \U00D7 0.2 m'^'2'))) 
stripchart(Count ~ Group, data=zapcon,
           vertical = TRUE, method = "jitter", 
           pch = 21, col = "maroon", bg = "bisque", 
           add = TRUE) 



#library(evd)
fitzapev <- fgev(zapateromd)

par(mfrow=c(2,2))
plot(fitzapev)

par(mfrow=)
hist(zapateromd,freq=F,col="gray",
     xlab=expression(paste('Count per 0.2 \U00D7 0.2 m'^'2')), ylab="Density", main="")
xval1 <- seq(from=0, to=max(zapateromd), length.out=200)

param2z <- fitzapev$estimate
locz <- param2z[["loc"]]
scalz <- param2z[["scale"]]
shapez <- param2z[["shape"]]
lines(xval1, dgev(xval1, loc=locz, scale=scalz, shape=shapez), col="blue", lwd=2)

thetaZapateroEV <- rgev(1000000, loc=locz, shape=shapez, scale=scalz)
mean(thetaZapateroEV)
sd(thetaZapateroEV)
thetaZapateroEV[thetaZapateroEV < 0] <- 0
thetaZapateroEVHA <- thetaZapateroEV*10000/(0.2*0.2)
thetaZapateroEVHA <- thetaZapateroEVHA[thetaZapateroEVHA >= 0 & thetaZapateroEVHA <= 100000000]
#thetaZapateroEVHA[thetaZapateroEVHA>100000000] <- 100000000
thetaZapateroEVHAm <- thetaZapateroEVHA/1000000
summary(thetaZapateroEVHAm)
dZapateroEVHAm <- density(thetaZapateroEVHAm, from=0)

summary(thetaZapateroEVHAm)
log(sd(thetaZapateroEVHAm))
log(mean(thetaZapateroEVHAm))
summary(zapateromd)
plot(dZapateroEVHAm)


#outlier analysis
plot(ecdf(conejosmd))
plot(ecdf(zapateromd))

chisq.out.test(conejosmd, variance=var(conejosmd), opposite = FALSE) #checking to see if highest value is an outlier
chisq.out.test(zapateromd, variance=var(zapateromd), opposite = FALSE)
dixon.test(conejosmd, type = 0, opposite = FALSE, two.sided = FALSE)
dixon.test(zapateromd, type = 0, opposite = FALSE, two.sided = FALSE)
grubbs.test(conejosmd, type = 20, opposite = FALSE, two.sided = FALSE)
grubbs.test(zapateromd, type = 20, opposite = FALSE, two.sided = FALSE)
qgrubbs(c(0.025,0.25,0.5,0.75,0.975), 29, type = 20, rev = FALSE)

qqnorm(log(zapateromd+0.001))
qqnorm(log(conejosmd))
qqnorm(zapateromd)
qqnorm(conejosmd)

qex <- function(x) qexp((rank(x)-.375)/(length(x)+.25))
plot(qex(zapateromd+0.001),log(zapateromd+0.001))
plot(qex(conejosmd),log(conejosmd))
conejosmd1 <- rm.outlier(conejosmd, fill=FALSE)
points(qex(conejosmd1),log(conejosmd1), pch=0, add=TRUE)
conejosmd2 <- rm.outlier(conejosmd1, fill=FALSE)
points(qex(conejosmd2),log(conejosmd2), pch=3, add=TRUE)
qqnorm(log(conejosmd2))
zapateromd1 <- rm.outlier(zapateromd, fill=FALSE)
zapateromd2 <- rm.outlier(zapateromd1, fill=FALSE)


#library(fitdistrplus)
zapateromd <- zapateromd+0.0001
# hist(zapateromd, freq=F, ylim=c(0, 0.075))
# fitzap<-fitdistr(zapateromd,"log-normal")$estimate
# lines(dlnorm(0:max(zapateromd),fitzap[1],fitzap[2]), lwd=3)
# fitzapg<-fitdistr(zapateromd,"gamma")$estimate
# lines(dgamma(0:max(zapateromd),fitzapg[1],fitzapg[2]), lty=2, lwd=3)
# 
# fitdistr(zapateromd,densfun=dgamma,start=list(scale=1,shape=2))


plotdist(zapateromd, histo = TRUE, demp = TRUE)
plotdist(conejosmd, histo=TRUE, demp=TRUE)

descdist(zapateromd, discrete=FALSE, boot=500)
descdist(conejosmd, discrete=FALSE, boot=500)

fitzap_w  <- fitdist(zapateromd, "weibull")
fitzap_g  <- fitdist(zapateromd, "gamma")
fitzap_ln <- fitdist(zapateromd, "lnorm")
summary(fitzap_ln)
summary(fitzap_w)
summary(fitzap_g)

fitzap_g$aic
fitzap_ln$aic
fitzap_w$aic

par(mfrow=c(2,2))
plot.legend <- c("Weibull", "lognormal", "gamma")
denscomp(list(fitzap_w, fitzap_g, fitzap_ln), legendtext = plot.legend)
cdfcomp (list(fitzap_w, fitzap_g, fitzap_ln), legendtext = plot.legend)
qqcomp  (list(fitzap_w, fitzap_g, fitzap_ln), legendtext = plot.legend)
ppcomp  (list(fitzap_w, fitzap_g, fitzap_ln), legendtext = plot.legend)


shapeZapatero <- fitzap_g$estimate[1]
rateZapatero <- fitzap_g$estimate[2]
thetaZapatero <- rgamma(1000000, shape = shapeZapatero, rate = rateZapatero)

dZapatero <- density(thetaZapatero, from=0)


fitcon_w  <- fitdist(conejosmd, "weibull")
#fitcon_g  <- fitdist(conejosmd, "gamma")
fitcon_ln <- fitdist(conejosmd, "lnorm")
fitcon_lnr <- fitdistr(conejosmd, "lognormal")
summary(fitcon_ln)
fitcon_lnr
summary(fitcon_w)

fitcon_ln$aic
fitcon_w$aic

par(mfrow=c(2,2))
plot.legend <- c("Weibull", "lognormal")
denscomp(list(fitcon_w, fitcon_ln), legendtext = plot.legend)
cdfcomp (list(fitcon_w, fitcon_ln), legendtext = plot.legend)
qqcomp  (list(fitcon_w, fitcon_ln), legendtext = plot.legend)
ppcomp  (list(fitcon_w, fitcon_ln), legendtext = plot.legend)

plot.legend <- c("lognormal")
denscomp(list(fitcon_ln), legendtext = plot.legend)
cdfcomp (list(fitcon_ln), legendtext = plot.legend)
qqcomp  (list(fitcon_ln), legendtext = plot.legend)
ppcomp  (list(fitcon_ln), legendtext = plot.legend)


lmeanConejos <- fitcon_ln$estimate[1]
SDConejos <- fitcon_ln$estimate[2]
shapeConejos <- fitcon_w$estimate[1]
scaleConejos <- fitcon_w$estimate[2]

thetaConejos <- rlnorm(1000000, lmeanConejos, SDConejos)   #Draw from a lognormal with estimated logmean and SD
weibullConejos <- rweibull(1000000,shapeConejos,scaleConejos)

thetaConejosHA <- thetaConejos*10000/(0.2*0.2)   #0.2 * 0.2 m^2 equals 0.000116 ha
weibullConejosHA <- weibullConejos*10000/(0.2*0.2) 
  
thetaConejosHA <- thetaConejosHA[thetaConejosHA >= 0 & thetaConejosHA <= 100000000]
weibullConejosHA <-  weibullConejosHA[weibullConejosHA >= 0 & weibullConejosHA <= 100000000] 
  
#thetaConejosHA[thetaConejosHA>100000000] <- 100000000      #Remove extremes exceeding 100 million
thetaConejosHAm <- thetaConejosHA/1000000                  #Express in numbers per million
weibullConejosHAm <- weibullConejosHA/1000000
#thetaConejosHAm <- thetaConejosHAm * 2  #because only 50% of Conejos sample was dead or moribund, other half was alive
#thetaConejosHAm <- thetaConejosHAm      #PeerJ review by Brower indicates that the 50% is incorrect
dConejosHAm <- density(thetaConejosHAm, from=0)
dweibullConejosHAm <- density(weibullConejosHAm, from=0)

summary(thetaConejosHA)
summary(weibullConejosHA)
log(sd(thetaConejosHAm))
sd(thetaConejosHAm)
sd(weibullConejosHAm)
log(mean(thetaConejosHAm))

summary(weibullConejosHA)

fitdist(thetaConejosHAm, "lnorm")











#Figure 1
png(filename="//../Figure1.png", 
    units="in", 
    width=7, 
    height=4, 
    pointsize=12, 
    res=900)

old1.par <- par(mfrow = c(3,1),
                oma = c(2,1,1,0) + 0.1,
                mar = c(4,4.1,0,0) + 0.1)


#####Structural

plot(0, xlim=c(0,100), ylim=c(0,0.08), main=NA, bty="n", cex.lab=1.6, cex.axis=1.4, pch="", xlab="", ylab="Density")
lines(density(thetabranch), xlim=c(0,100), ylim=c(0,0.08), main=NA, lwd=2, lty=1, ylab="", yaxt="n", xaxt="n", bty="n")
ablineclip(v=10.3, y1=0.001, y2=0.065, lwd=2, lty=1, col="dark gray")
text(16,0.075, "Branch", cex=1.4)
#text(50, 0.06, "Structural", cex=1.6)
mtext("A", side=3, line=-0.8, adj=-0.14, cex=1.8)

######Capture Mark Recapture

plot(0, xlim=c(0,100), ylim=c(0,0.35), main=NA, cex.lab=1.6, cex.axis=1.4, bty="n", pch="", xlab="", ylab="Density")
lines(density(thetaDec), xlim=c(0,100), ylim=c(0,0.35), main=NA, lwd=2, lty=2, ylab="", yaxt="n", xaxt="n", bty="n")
lines(density(thetaP), xlim=c(0,100), ylim=c(0,0.35), main=NA, bty="n", lwd=2, lty=4)
lines(density(thetaJS), xlim=c(0,100), ylim=c(0,0.35), main=NA, lwd=2, lty=1, ylab="", yaxt="n", xaxt="n", bty="n")
ablineclip(v=60.9,y1=0.003, y2=0.05, lwd=2, lty=1, col="dark gray")
ablineclip(v=33.8, y1=0.003, y2=0.05, lwd=2, lty=1, col="dark gray")
ablineclip(v=6.9, y1=0.003, y2=0.3, lwd=2, lty=1, col="dark gray")
text(18,0.25, "Petersen (Dec)", cex=1.4)
text(40,0.09, "Jolly-Seber", cex=1.4)
text(69,0.09, "Petersen (Jan)", cex=1.4)
#text(50, 0.30, "Capture-Mark-Recapture", cex=1.6)
mtext("B", side=3, line=-0.8, adj=-0.14, cex=1.8)

#######January storm-mortality est.

plot(0, xlim=c(0,100), ylim=c(0,0.05), main=NA, bty="n", cex.lab=1.6, cex.axis=1.4, pch="", xlab="Population Density (millions/ha)", ylab="Density")
lines(dZapateroEVHAm, lwd=2, main="", sub="")
text(24,0.035, "Zapatero", cex=1.4)
ablineclip(v=median(thetaZapateroEVHAm), y1=0.003, y2=0.027, lwd=2, col="dark gray")
lines(dConejosHAm, lwd=2, lty=2, main="", sub="")
text(10,0.013, "Conejos", cex=1.4)
ablineclip(v=median(thetaConejosHAm), y1=0.003, y2=0.019, lwd=2, col="dark gray")
#text(50, 0.036, "Storm Mortality", cex=1.6)
mtext("C", side=3, line=-0.8, adj=-0.14, cex=1.8)


dev.off() 



######
#The number of samples from the mixture distribution
N = 1000000    

#Sample N random uniforms U
U = runif(N)

#Variable to store the samples from the mixture distribution                                             
rand.samples = rep(NA,N)

#Sampling from the mixture
for(i in 1:N){
  if(U[i]<0.1666667){
    rand.samples[i] = rlnorm(1, lmeanP, SDP)
  }else if(U[i]<0.333333){
    rand.samples[i] = rlnorm(1, lmeanDec, SDDec)
  }else if(U[i]<0.50){
    rand.samples[i] = rlnorm(1, lmeanbranch, SDbranch)
  }else if(U[i]<0.66666667){
    rand.samples[i] = rlnorm(1, lmeanJS, SDJS)
  }else if(U[i]<0.83333333){
    rand.samples[i] = rlnorm(1, lmeanConejos, SDConejos)
    rand.samples[i] <- rand.samples[i]*10000/(0.2*0.2)
    rand.samples[i] <- rand.samples[i]/1000000
  }else {
    rand.samples[i] = rgev(1, loc=locz, shape=shapez, scale=scalz)
    rand.samples[i] <- rand.samples[i]*10000/(0.2*0.2)
    rand.samples[i] <- rand.samples[i]/1000000
  }}

max(rand.samples)
min(rand.samples)
rand.samples <- rand.samples[rand.samples >= 0 & rand.samples <= 100]

summary(rand.samples)
quantile(rand.samples, c(0.025, 0.5, 0.975))

fitrandsamp_ln <- fitdist(rand.samples, "lnorm")
fitrandsamp_lnr <- fitdistr(rand.samples, "lognormal")  #Equivalent to above procedure

summary(fitrandsamp_ln)
fitrandsamp_lnr
exp(fitrandsamp_ln$estimate[1])
plot(fitrandsamp_ln)

#Figure 2
png(filename="//.../Figure2.png", 
    units="in", 
    width=7, 
    height=4, 
    pointsize=12, 
    res=900)

par(mfrow=c(1,1))

#Density plot of the random samples   xlim=c(0,100), ylim=c(0,0.1), 
plot(0, main=NA, bty="n",   xlim=c(0,100), ylim=c(0,0.07), pch="", xlab="Population Density (millions/ha)", ylab="Density")
lines(density(rand.samples),    xlim=c(0,100), ylim=c(0,0.07), main=NA, bty="n", lwd=3, lty=1)
segments(60.9, 0, 60.9, 0.012, lwd=2, lty=1, col="blue")  #Petersen (Jan)
segments(33.8, 0, 33.8, 0.020, lwd=2, lty=1, col="blue")  #Jolly Seber
text(42,0.025, "Jolly-Seber", cex=1.2)
text(71,0.017, "Petersen (Jan)", cex=1.2)
segments(mean(thetaConejosHAm), 0, mean(thetaConejosHAm), 0.048, lwd=2, lty=1, col="blue")  #Conejos
segments(mean(thetaZapateroEVHAm), 0, mean(thetaZapateroEVHAm), 0.035, lwd=2, lty=1, col="blue")  #Zapatero
text(31,0.053, "Conejos", cex=1.2)
text(20,0.040, "Zapatero", cex=1.2)
segments(mean(rand.samples), 0, mean(rand.samples), 0.028, lwd=2, lty=2, col="black")
#segments(quantile(rand.samples, probs=c(0.5)), 0, quantile(rand.samples, probs=c(0.5)), 0.048, lwd=2, lty=1, col="black")
shadowtext(30,0.031, "Mean", col="black", bg="white", cex=1.2)
# segments(50, 0, 50, 0.025, lwd=2, lty=1, col="blue")  #Brower
# text(53,0.0281, "Brower", cex=1.2)
segments(10.3, 0, 10.3, 0.045, lwd=2, lty=1, col="blue")  #Branch
segments(6.9, 0, 6.9, 0.062, lwd=2, lty=1, col="blue")  #Petersen (Dec)
text(18,0.067, "Petersen (Dec)", cex=1.2)
text(16,0.050, "Branch", cex=1.2)

#Demarcation between distributions
#segments(17.9, 0, 17.9, 0.01, lwd=3, lty=2, col="dark gray")

dev.off()


mean(rand.samples)
summary(rand.samples)
quantile(rand.samples, probs=c(0.025, 0.5, 0.975))
# sum(rand.samples)
# head(rand.samples)


#sum(rand.samples <= 17.8)/(sum(rand.samples <= 17.8)+sum(rand.samples > 17.8))

# rand.low <- rand.samples[rand.samples <= 17.8]
# mean(rand.low)
# quantile(rand.low, probs=c(0.025, 0.5, 0.975))
# rand.high <- rand.samples[rand.samples > 17.8]
# mean(rand.high)
# quantile(rand.high, probs=c(0.025, 0.5, 0.975))
# 
# summary(rand.low)
# summary(rand.high)










#######Factors affecting density
###################
means <- c(6.9, 60.9, 33.8, 24.945, 23.6147, 10.3)
#y6 <- c( (6.9-min(y)), (60.9-min(y)), (33.8-min(y)), (25.9589-min(y)), (23.58251-min(y)), (10.3-min(y)) )
days <- c(8, 42, 33, 37, 37, 52)             #x equals days of year, with 12 Dec = 1, 20 Dec = 8, 22 Jan = 42, 1 Feb = 52
temperature <- c(49.8, 49.813, 48.368, 48.365, 48.365, 50.561)
dewp <- c(35.87391, 33.82174, 33.08636, 34.8, 34.8, 35.73043)



options(digits=4)
options(scipen=999)

#library(lattice)

degrees.C= function(x) {
  tc=(x-32)*5/9
  return(tc)
}

degrees.F =function(y) {
  tf=9*y/5+32
  return(tf)
}

dewp.C <- degrees.C(dewp)
dataOW <- as.data.frame(cbind(means, days, temperature, dewp, dewp.C))

Cand.models <- list( )
Cand.models[[1]] <- lm(means~days) #Day of Year
# 
# fit.eff <- allEffects(fit)
# plot(fit.eff)
Cand.models[[2]]  <- lm(means~days+temperature+dewp.C) #Day of Year + Env
Cand.models[[3]]  <- lm(means~temperature+dewp.C) #Env
Cand.models[[4]]  <- lm(means~days+temperature) #Day of Year + Temperature
Cand.models[[5]]  <- lm(means~days+dewp.C) #Day of Year + Dew Point
Cand.models[[6]]  <- lm(means~dewp.C) #Dew Point
Cand.models[[7]]  <- lm(means~temperature) #Temperature


#library(AICcmodavg)

aictab(cand.set = Cand.models, modnames = c("Day", "Day+Temp+Dew", "Temp+Dew", "Day+Temp", "Day+Dew", "Dew", "Temp"), sort = TRUE, c.hat = 1, second.ord = TRUE,
       nobs = NULL)

#library(effects)
summary(Cand.models[[6]])
Cand.models6.eff <- allEffects(Cand.models[[6]], confidence.level=0.9)
effect("dewp.C", Cand.models[[6]])

trellis.par.get()
axis.text <- trellis.par.get("axis.text")
par.ylab.text <- trellis.par.get("par.ylab.text")
par.xlab.text <- trellis.par.get("par.xlab.text")

axis.text$cex <- 2
par.ylab.text$cex <- 2
par.xlab.text$cex <- 2



####Fig 3
library(grid)

png(filename="//.../Figure3a.png", 
    units="in", 
         width=5, 
         height=5, 
         pointsize=12, 
         res=900)

plot(Cand.models6.eff, ylim=c(0, 70), main="", cex.axis=2, cex=2,
     xlab= expression("Mean Dew Point ("* degree *"C)"), ylab="Density (millions/ha)")
trellis.focus("toplevel") ## has coordinate system [0,1] x [0,1]
panel.text(0.55, .85, labels="Density = 64.5 - 25.3 Dew Point", cex=1.3)
panel.text(0.55, .80, labels= "(adj. R-squ = 0.51, p = 0.07)", cex=1.3)
trellis.unfocus()
dev.off()

Cand.models7.eff <- allEffects(Cand.models[[7]])
effect("dewp", Cand.models[[7]])
plot(Cand.models7.eff, ylim=c(0, 90), xlab= expression("Mean Temperature ("* degree *"F)"), ylab="Density (millions/ha)")


fitdewdate <- lm(days~dewp.C)
summary(fitdewdate)
fitdewdate.eff <- allEffects(fitdewdate)
effect("dewp.C", fitdewdate)
plot(fitdewdate.eff)

fitdatedew <- lm(dewp.C~days)
summary(fitdatedew)
fitdatedew.eff <- allEffects(fitdatedew)
effect("days", fitdatedew)
plot(fitdatedew.eff)


dewpoints <- data.frame(dewp.C = seq(from=1.6, to=1.8, length.out=6))
pred.w.plim <- predict(fitg.4, dewpoints, interval = "prediction")
pred.w.clim <- predict(fitg.4, dewpoints, interval = "confidence")
matplot(dewpoints$dewp.C, cbind(pred.w.clim, pred.w.plim[,-1]), ylim=c(0,100),
        lty = c(1,2,2,3,3), col=c(1,1,1,2,2), lwd=c(2,1,1,1,1), type = "l", 
        xlab = expression("Mean Dew Point ("* degree *"C)"), ylab="Density (millions/ha)")



png(filename="//.../Figure3b.png", 
    units="in", 
    width=7, 
    height=4, 
    pointsize=12, 
    res=900)

ggplot(subset(OWtemp, (Month=="11"|Month=="12"|Month=="1"|Month=="2"|Month=="3")), 
                       aes(x=Day, y=DEW.C, color=factor(month(yearmonthday)))) + 
  theme_bw() +
  geom_boxplot(fill = "grey80", colour = "black", outlier.shape = NA) + 
  geom_point(size = 2, alpha = 3/10) +
  facet_grid(. ~ fall2fall) +
  ylab("Mean Dew Point") + theme(legend.position="none")

dev.off()





################Density effect on population size estimates
#library(matrixStats)
#View(N.est)    ####From Semmens et al. (2016)

N.estsummary <- colQuantiles(N.est[,2:23], probs=c(0.025, 0.5, 0.975))
exp(N.estsummary)
N.estPop <- exp(N.est[,2:23])*quantile(rand.samples, probs=c(0.5))
summary(N.estPop)
quantPop <- colQuantiles(N.estPop, probs=c(0.025, 0.5, 0.975))

# NOW PLOT THE TIME SERIES ################################################################
y<-log(c(6.23,7.81,12.61,18.19,5.77,5.56,8.97,2.83,9.36,7.54,11.12,2.19,5.91,
         6.87,4.61,5.06,1.92,4.02,2.89,1.19,0.67,1.13))
n.years<-length(y)
y.y<-seq(1993,2014)
N.yr<-n.years

fitted<-lower<-upper<-numeric() 

# turn the following 3 lines on to plot in non-log space
#N.est<-exp(N.est)
#y<-exp(y)

for (i in 1:n.years){
  fitted[i]<-quantile(N.estPop[,i],0.50)
  lower[i]<-quantile(N.estPop[,i],0.025)
  upper[i]<-quantile(N.estPop[,i],0.975)
  
}
m1<-min(c(fitted,lower))
m2<-max(c(fitted,upper)) 



#Figure 4
png(filename="//.../Figure4.png", 
    units="in", 
    width=7, 
    height=4, 
    pointsize=12, 
    res=900)



par(mar=c(4.5,4,1,1),cex=1.2)

plot(0,0,ylim=c(m1,m2),xlim=c(min(y.y),max(y.y)+1),ylab="Monarch Population (millions)",
     xlab="Year",las=1,col="black",
     type="l",lwd=2,frame=FALSE,axes=FALSE)
axis(2,las=1)
axis(1,at=seq(min(y.y),max(y.y)+1,2), labels=seq(min(y.y),max(y.y)+1,2))
axis(1,at=0:N.yr,labels=rep("",N.yr+1),tcl=-0.25)
polygon(x=c(min(y.y):max(y.y),max(y.y):min(y.y)),y=c(lower,upper[N.yr:1]), col= rgb(0, 0, 1, 0.25), border=NA)

points(y.y,fitted,type="l",col= "black" , lwd=2)

abline(h=0.25*quantile(rand.samples, probs=c(0.5)), lwd=2, lty=4, col="dark gray" )  
#Upper-most quasi-extinction threshold as described by Semmens et al. (2016)



dev.off()

quantile(rand.samples, probs=c(0.5, 0.025, 0.975))*6

##########

meanlist <- c(quantile(N.estPop[,1],0.50),
              quantile(N.estPop[,2],0.50),
              quantile(N.estPop[,3],0.50),
              quantile(N.estPop[,4],0.50),
              quantile(N.estPop[,5],0.50),
              quantile(N.estPop[,6],0.50),
              quantile(N.estPop[,7],0.50),
              quantile(N.estPop[,8],0.50),
              quantile(N.estPop[,9],0.50),
              quantile(N.estPop[,10],0.50),
              quantile(N.estPop[,11],0.50),
              quantile(N.estPop[,12],0.50),
              quantile(N.estPop[,13],0.50),
              quantile(N.estPop[,14],0.50),
              quantile(N.estPop[,15],0.50),
              quantile(N.estPop[,16],0.50),
              quantile(N.estPop[,17],0.50),
              quantile(N.estPop[,18],0.50),
              quantile(N.estPop[,19],0.50),
              quantile(N.estPop[,20],0.50),
              quantile(N.estPop[,21],0.50),
              quantile(N.estPop[,22],0.50))


mean(quantile(N.estPop[,1],0.50),
     quantile(N.estPop[,2],0.50),
     quantile(N.estPop[,3],0.50),
     quantile(N.estPop[,4],0.50),
     quantile(N.estPop[,5],0.50),
     quantile(N.estPop[,6],0.50),
     quantile(N.estPop[,7],0.50),
     quantile(N.estPop[,8],0.50),
     quantile(N.estPop[,9],0.50),
     quantile(N.estPop[,10],0.50),
     quantile(N.estPop[,11],0.50),
     quantile(N.estPop[,12],0.50),
     quantile(N.estPop[,13],0.50),
     quantile(N.estPop[,14],0.50),
     quantile(N.estPop[,15],0.50),
     quantile(N.estPop[,16],0.50),
     quantile(N.estPop[,17],0.50),
     quantile(N.estPop[,18],0.50),
     quantile(N.estPop[,19],0.50),
     quantile(N.estPop[,20],0.50),
     quantile(N.estPop[,21],0.50),
     quantile(N.estPop[,22],0.50))

#mean(N.estPop) #deprecated
j <- unlist(N.estPop)
mean(j)
summary(N.estPop)
summary(j)

quantile(j,c(0.025,0.5,0.975))


lowerlist <- c(quantile(N.estPop[,1],0.0250),
               quantile(N.estPop[,2],0.0250),
               quantile(N.estPop[,3],0.0250),
               quantile(N.estPop[,4],0.0250),
               quantile(N.estPop[,5],0.0250),
               quantile(N.estPop[,6],0.0250),
               quantile(N.estPop[,7],0.0250),
               quantile(N.estPop[,8],0.0250),
               quantile(N.estPop[,9],0.0250),
               quantile(N.estPop[,10],0.0250),
               quantile(N.estPop[,11],0.0250),
               quantile(N.estPop[,12],0.0250),
               quantile(N.estPop[,13],0.0250),
               quantile(N.estPop[,14],0.0250),
               quantile(N.estPop[,15],0.0250),
               quantile(N.estPop[,16],0.0250),
               quantile(N.estPop[,17],0.0250),
               quantile(N.estPop[,18],0.0250),
               quantile(N.estPop[,19],0.0250),
               quantile(N.estPop[,20],0.0250),
               quantile(N.estPop[,21],0.0250),
               quantile(N.estPop[,22],0.0250))

upperlist <- c(quantile(N.estPop[,1],0.9750),
               quantile(N.estPop[,2],0.9750),
               quantile(N.estPop[,3],0.9750),
               quantile(N.estPop[,4],0.9750),
               quantile(N.estPop[,5],0.9750),
               quantile(N.estPop[,6],0.9750),
               quantile(N.estPop[,7],0.9750),
               quantile(N.estPop[,8],0.9750),
               quantile(N.estPop[,9],0.9750),
               quantile(N.estPop[,10],0.9750),
               quantile(N.estPop[,11],0.9750),
               quantile(N.estPop[,12],0.9750),
               quantile(N.estPop[,13],0.9750),
               quantile(N.estPop[,14],0.9750),
               quantile(N.estPop[,15],0.9750),
               quantile(N.estPop[,16],0.9750),
               quantile(N.estPop[,17],0.9750),
               quantile(N.estPop[,18],0.9750),
               quantile(N.estPop[,19],0.9750),
               quantile(N.estPop[,20],0.9750),
               quantile(N.estPop[,21],0.9750),
               quantile(N.estPop[,22],0.9750))



idlist <- 1994:2015
poplists <- cbind(meanlist, lowerlist, upperlist)
row.names(poplists) <- idlist
colnames(poplists) <- c("Mean", "Lower", "Upper")

l = signif(poplists,4)
dfl <- data.frame(matrix(unlist(l), nrow=22, byrow=T),stringsAsFactors=FALSE)
dfl

write.table(poplists, "clipboard", sep="\t", row.names=FALSE, col.names=FALSE)


mean(l[,1])

mean(l[,2])

mean(l[,3])


quantile(rand.samples, probs=c(0.2, 0.5, 0.9))

mean(rand.samples)*6 
quantile(rand.samples, probs=c(0.025, 0.5, 0.975))*6 

mean(rand.samples)*6 *28.5
quantile(rand.samples, probs=c(0.025, 0.5, 0.975))*6 *28.5

(quantile(N.estPop[,6],probs=c(0.025, 0.5, 0.975)) *28.5) - (quantile(N.estPop[,19], probs=c(0.025, 0.5, 0.975)) *28.5)















##################5 estimate approach using published Brower mean estimate of 50 million monarchs per ha
#Brower, (1/mortality rate)*no. of monarchs killed per m2*m2/ha=Total no. of monarchs/ha before storm
# 90.55/2 #Conejos colony, (1/0.801)*7253*10000
# 35.37/2 #Zapatero colony, (1/0.743)*2626*10000
# (90.55+35.37)/2    #Brower et al originally reported 60 million per ha, but revised this estimate down in subsequent publications


a <- 35.37
c <- 50   #Brower mean
b <- 90.55




Browermin <- a
Browerlikely <- c
Browermax <- b

#The number of samples from the mixture distribution
N = 1000000      

#Sample N random uniforms U
U = runif(N)

#Variable to store the samples from the mixture distribution                                                                   
rand.samples = rep(NA,N)

#Sampling from the mixture
for(i in 1:N){
   if(U[i]<0.25){
      rand.samples[i] = rlnorm(1, lmeanP, SDP)
   }else if(U[i]<0.50){
      rand.samples[i] = rlnorm(1, lmeanJS, SDJS)
   }else if(U[i]<0.75){
      rand.samples[i] = rlnorm(1, lmeanB, SDP)
   }else {
      rand.samples[i] = rlnorm(1, lmeanBrower, SDBrower)
   }}