#--------------------------------------------------------------------------------------------------------------
# R and JAGS code for estimating total woodcock harvest, 1964-2016 (53 years)
# Combined analysis of duck stamp survey (DSS: 1964-2001) and 
#   Harvest Information Program (HIP: 1999-2016)
#   reconciling different sampling frames of hunters
# Some states pooled (20 jurisdictions total)
#
# Todd Arnold, Proc. 11th American Woodcock Symposium, arnol065@umn.edu
#--------------------------------------------------------------------------------------------------------------

library(jagsUI)

# LOAD DATA #
# modify working directory to specify path on your own computer
setwd("C:/Users/arnol065/Documents/projects/submitted/AMWO_harvest/data") 
# these 2 files represent historical data
stamps.data <- read.csv("Stamps.csv",header=TRUE)
DSS.data <- read.csv("DSS.csv",header=TRUE)

# these 2 files are current through 2016 harvest
 # analysis can be updated by inserting additional columns in each file
HIP.data <- read.csv("HIP.csv",header=TRUE)
HIP_var.data <- read.csv("HIP_var.csv",header=TRUE)

n.states.orig <- dim(HIP.data)[1] # data from 35 states
n.states <- 20 # consolidate to 20 states
n.HIP <- dim(HIP.data)[2]-1
n.DSS <- dim(DSS.data)[2]-1
n.yrs <- n.HIP + n.DSS - 3 # subtract off 3 overlap years (1999-2001)
n.stamps <- dim(stamps.data)[2]-1

# replace state specific NA & 0 with small but non-zero values
# Paul Padding, pers. comm., these values are true point estimates of zero harvest
# first convert 0 to NA
DSS.data[DSS.data == 0] <- NA
HIP.data[HIP.data == 0] <- NA
HIP_var.data[HIP_var.data == 0] <- NA
# then convert NA to min non-zero value observed for that state
for(i in 1:n.states.orig){
  for (j in 2:(n.HIP+1)){
    DSS.data[i,j][is.na(DSS.data[i,j] ) ] = apply(HIP.data[i,2:19],1,min,na.rm=TRUE) 
    HIP.data[i,j][is.na(HIP.data[i,j] ) ] = apply(HIP.data[i,2:19],1,min,na.rm=TRUE) 
    HIP_var.data[i,j][is.na(HIP_var.data[i,j] ) ] = apply(HIP_var.data[i,2:19],1,min,na.rm=TRUE) 
  }}

## read data into 20x53 matrices for jags (combined state specific data)
# col 1 is 2016, col 53 is 1964 (reverse time)
HIP2 <- HIP_var2 <- DSS2 <- DSS_var2 <- stamps2 <- matrix(NA,nrow=n.states,ncol=n.yrs)

for (j in 1:n.HIP){
  HIP2[1,j] <- HIP.data[12,(j+1)] # LA
  HIP_var2[1,j] <- HIP_var.data[12,(j+1)] 
  HIP2[2,j] <- HIP.data[16,(j+1)] # MI
  HIP_var2[2,j] <- HIP_var.data[16,(j+1)] 
  HIP2[3,j] <- HIP.data[34,(j+1)] # WI
  HIP_var2[3,j] <- HIP_var.data[34,(j+1)] 
  HIP2[4,j] <- HIP.data[17,(j+1)] # MN
  HIP_var2[4,j] <- HIP_var.data[17,(j+1)] 
  HIP2[5,j] <- HIP.data[25,(j+1)] # OH
  HIP_var2[5,j] <- HIP_var.data[25,(j+1)] 
  HIP2[6,j] <- HIP.data[8,(j+1)] + HIP.data[9,(j+1)] # IL, IN
  HIP_var2[6,j] <- HIP_var.data[8,(j+1)] + HIP_var.data[9,(j+1)] # add variances
  HIP2[7,j] <- HIP.data[2,(j+1)] + HIP.data[7,(j+1)] + HIP.data[18,(j+1)] # AR, IA, MO
  HIP_var2[7,j] <- HIP_var.data[2,(j+1)] + HIP_var.data[7,(j+1)] + HIP_var.data[18,(j+1)] 
  HIP2[8,j] <- HIP.data[11,(j+1)] + HIP.data[30,(j+1)] # KY, TN
  HIP_var2[8,j] <- HIP_var.data[11,(j+1)] + HIP_var.data[30,(j+1)] 
  HIP2[9,j] <- HIP.data[1,(j+1)] + HIP.data[19,(j+1)] # AL, MS
  HIP_var2[9,j] <- HIP_var.data[1,(j+1)] + HIP_var.data[19,(j+1)] 
  HIP2[10,j] <- HIP.data[10,(j+1)] + HIP.data[21,(j+1)] + HIP.data[26,(j+1)] + HIP.data[31,(j+1)] # KS, NE, OK, TX
  HIP_var2[10,j] <- HIP_var.data[10,(j+1)] + HIP_var.data[21,(j+1)] + HIP_var.data[26,(j+1)] + HIP_var.data[31,(j+1)]
  HIP2[11,j] <- HIP.data[24,(j+1)] # NY
  HIP_var2[11,j] <- HIP_var.data[24,(j+1)] 
  HIP2[12,j] <- HIP.data[15,(j+1)] # ME
  HIP_var2[12,j] <- HIP_var.data[15,(j+1)] 
  HIP2[13,j] <- HIP.data[23,(j+1)]  # NJ
  HIP_var2[13,j] <- HIP_var.data[23,(j+1)]  
  HIP2[14,j] <- HIP.data[27,(j+1)] # PA
  HIP_var2[14,j] <- HIP_var.data[27,(j+1)] 
  HIP2[15,j] <- HIP.data[13,(j+1)]  # MA
  HIP_var2[15,j] <- HIP_var.data[13,(j+1)]  
  HIP2[16,j] <- HIP.data[22,(j+1)] + HIP.data[33,(j+1)] # NH, VT
  HIP_var2[16,j] <- HIP_var.data[22,(j+1)] + HIP_var.data[33,(j+1)] 
  HIP2[17,j] <- HIP.data[3,(j+1)] + HIP.data[28,(j+1)] # CT, RI
  HIP_var2[17,j] <- HIP_var.data[3,(j+1)] + HIP_var.data[28,(j+1)] 
  HIP2[18,j] <- HIP.data[32,(j+1)] + HIP.data[35,(j+1)] + HIP.data[4,(j+1)] + HIP.data[14,(j+1)] # DE, MD, VA, WV
  HIP_var2[18,j] <- HIP_var.data[32,(j+1)] + HIP_var.data[35,(j+1)] + HIP_var.data[4,(j+1)] + HIP_var.data[14,(j+1)]
  HIP2[19,j] <- HIP.data[20,(j+1)] + HIP.data[29,(j+1)] # NC, SC
  HIP_var2[19,j] <- HIP_var.data[20,(j+1)] + HIP_var.data[29,(j+1)] 
  HIP2[20,j] <- HIP.data[5,(j+1)] + HIP.data[6,(j+1)] # FL, GA
  HIP_var2[20,j] <- HIP_var.data[5,(j+1)] + HIP_var.data[6,(j+1)] 
}  
  for (j in 1:n.DSS){
    DSS2[1,j+15] <- DSS.data[12,(j+1)] # LA
    DSS_var2[1,j+15] <- exp(1.384 +1.359*log(DSS.data[12,(j+1)]+1))
    DSS2[2,j+15] <- DSS.data[16,(j+1)] # MI
    DSS_var2[2,j+15] <- exp(1.384 +1.359*log(DSS.data[16,(j+1)]+1))
    DSS2[3,j+15] <- DSS.data[34,(j+1)] # WI
    DSS_var2[3,j+15] <- exp(1.384 +1.359*log(DSS.data[34,(j+1)]+1))
    DSS2[4,j+15] <- DSS.data[17,(j+1)] # MN
    DSS_var2[4,j+15] <- exp(1.384 +1.359*log(DSS.data[17,(j+1)]+1))
    DSS2[5,j+15] <- DSS.data[25,(j+1)] # OH
    DSS_var2[5,j+15] <- exp(1.384 +1.359*log(DSS.data[25,(j+1)]+1))
    DSS2[6,j+15] <- DSS.data[8,(j+1)] + DSS.data[9,(j+1)] # IL, IN
    DSS_var2[6,j+15] <- exp(1.384 +1.359*log(DSS.data[8,(j+1)] + DSS.data[9,(j+1)]+1)) # add variances
    DSS2[7,j+15] <- DSS.data[2,(j+1)] + DSS.data[7,(j+1)] + DSS.data[18,(j+1)] # AR, IA, MO
    DSS_var2[7,j+15] <- exp(1.384 +1.359*log(DSS.data[2,(j+1)] + DSS.data[7,(j+1)] + DSS.data[18,(j+1)] + 1))
    DSS2[8,j+15] <- DSS.data[11,(j+1)] + DSS.data[30,(j+1)] # KY, TN
    DSS_var2[8,j+15] <- exp(1.384 +1.359*log(DSS.data[11,(j+1)] + DSS.data[30,(j+1)]+1))
    DSS2[9,j+15] <- DSS.data[1,(j+1)] + DSS.data[19,(j+1)] # AL, MS
    DSS_var2[9,j+15] <- exp(1.384 +1.359*log(DSS.data[1,(j+1)] + DSS.data[19,(j+1)]+1))
    DSS2[10,j+15] <- DSS.data[10,(j+1)] + DSS.data[21,(j+1)] + DSS.data[26,(j+1)] + DSS.data[31,(j+1)] # KS, NE, OK, TX
    DSS_var2[10,j+15] <- exp(1.384 +1.359*log(DSS.data[10,(j+1)] + DSS.data[21,(j+1)] + DSS.data[26,(j+1)] + DSS.data[31,(j+1)]+1)) 
    DSS2[11,j+15] <- DSS.data[24,(j+1)] # NY
    DSS_var2[11,j+15] <- exp(1.384 +1.359*log(DSS.data[24,(j+1)]+1))
    DSS2[12,j+15] <- DSS.data[15,(j+1)] # ME
    DSS_var2[12,j+15] <- exp(1.384 +1.359*log(DSS.data[15,(j+1)]+1))
    DSS2[13,j+15] <- DSS.data[23,(j+1)] # NJ
    DSS_var2[13,j+15] <- exp(1.384 +1.359*log(DSS.data[27,(j+1)]+1))
    DSS2[14,j+15] <- DSS.data[27,(j+1)] # PA
    DSS_var2[14,j+15] <- exp(1.384 +1.359*log(DSS.data[27,(j+1)]+1))
    DSS2[15,j+15] <- DSS.data[13,(j+1)] # MA
    DSS_var2[15,j+15] <- exp(1.384 +1.359*log(DSS.data[13,(j+1)]+1)) 
    DSS2[16,j+15] <- DSS.data[22,(j+1)] + DSS.data[33,(j+1)] # NH, VT
    DSS_var2[16,j+15] <- exp(1.384 +1.359*log(DSS.data[22,(j+1)] + DSS.data[33,(j+1)]+1))
    DSS2[17,j+15] <- DSS.data[3,(j+1)] + DSS.data[28,(j+1)] # CT, RI
    DSS_var2[17,j+15] <- exp(1.384 +1.359*log(DSS.data[3,(j+1)] + DSS.data[28,(j+1)]+1)) 
    DSS2[18,j+15] <- DSS.data[4,(j+1)] + DSS.data[14,(j+1)] + DSS.data[32,(j+1)] + DSS.data[35,(j+1)] # DE, MD, VA, WV
    DSS_var2[18,j+15] <- exp(1.384 +1.359*log(DSS.data[4,(j+1)] + DSS.data[14,(j+1)] + DSS.data[32,(j+1)] + DSS.data[35,(j+1)]+1))
    DSS2[19,j+15] <- DSS.data[20,(j+1)] + DSS.data[29,(j+1)] # NC, SC
    DSS_var2[19,j+15] <- exp(1.384 +1.359*log(DSS.data[20,(j+1)] + DSS.data[29,(j+1)]+1))
    DSS2[20,j+15] <- DSS.data[5,(j+1)] + DSS.data[6,(j+1)] # FL, GA
    DSS_var2[20,j+15] <- exp(1.384 +1.359*log(DSS.data[5,(j+1)] + DSS.data[6,(j+1)]+1))
  }    
  for (j in 6:(n.stamps+5)){
    stamps2[1,j] <- stamps.data[12,(j-4)] # LA
    stamps2[2,j] <- stamps.data[16,(j-4)] # MI
    stamps2[3,j] <- stamps.data[34,(j-4)] # WI
    stamps2[4,j] <- stamps.data[17,(j-4)] # MN
    stamps2[5,j] <- stamps.data[25,(j-4)] # OH
    stamps2[6,j] <- stamps.data[8,(j-4)] + stamps.data[9,(j-4)] # IL, IN
    stamps2[7,j] <- stamps.data[2,(j-4)] + stamps.data[7,(j-4)] + stamps.data[18,(j-4)] # AR, IA, MO
    stamps2[8,j] <- stamps.data[11,(j-4)] + stamps.data[30,(j-4)] # KY, TN
    stamps2[9,j] <- stamps.data[1,(j-4)] + stamps.data[19,(j-4)] # AL, MS
    stamps2[10,j] <- stamps.data[10,(j-4)] + stamps.data[21,(j-4)] + stamps.data[26,(j-4)] + stamps.data[31,(j-4)] # KS, NE, OK, TX
    stamps2[11,j] <- stamps.data[24,(j-4)] # NY
    stamps2[12,j] <- stamps.data[15,(j-4)] # ME
    stamps2[13,j] <- stamps.data[23,(j-4)] # NJ
    stamps2[14,j] <- stamps.data[27,(j-4)] # PA
    stamps2[15,j] <- stamps.data[13,(j-4)] # MA
    stamps2[16,j] <- stamps.data[22,(j-4)] + stamps.data[33,(j-4)] # NH, VT
    stamps2[17,j] <- stamps.data[3,(j-4)] + stamps.data[28,(j-4)] # CT, RI
    stamps2[18,j] <- stamps.data[4,(j-4)] + stamps.data[14,(j-4)] + stamps.data[32,(j-4)] + stamps.data[35,(j-4)] # DE, MD, VA, WV
    stamps2[19,j] <- stamps.data[20,(j-4)] + stamps.data[29,(j-4)] # NC, SC
    stamps2[20,j] <- stamps.data[5,(j-4)] + stamps.data[6,(j-4)] # FL, GA
  }

DSS_sd2 <- sqrt(DSS_var2)
HIP_sd2 <- sqrt(HIP_var2)

# plot HIP and DSS estimates with 95% CI derived from variance estimates 
 # Wisconsin = row 3, modify matrix notation "[3," to other rows to plot other states (e.g. "[1,]" for Louisiana)
plot(1964:2016, rev(DSS2[3,]/1000), ylab="Harvest (1000's)", xlab="", main="Wisconsin", col="red")
segments(1964:2016, rev((DSS2[3,]-1.28*DSS_sd2[3,])/1000), 1964:2016, rev((DSS2[3,]+1.28*DSS_sd2[3,])/1000), col="red")
points(1964:2016, rev(HIP2[3,]/1000), col="blue")
segments(1964:2016, rev((HIP2[3,]-1.28*HIP_sd2[3,])/1000), 1964:2016, rev((HIP2[3,]+1.28*HIP_sd2[3,])/1000),col="blue")

# correlate DSS and stamps
DSS_stamps.cor <- numeric()
for (i in 1:n.states){
  DSS_stamps.cor[i] <- cor(DSS2[i,16:53],stamps2[i,16:53])
}
print(round(DSS_stamps.cor,3))
# [1]  0.565  0.536  0.732 -0.062  0.356  0.336  0.493  0.356  0.551  0.424  0.928  0.826
# [13]  0.911  0.781  0.777  0.810  0.654  0.652  0.144  0.224
mean(DSS_stamps.cor[1:10]) # 0.429 arithmatic means, not weighted by harvest
mean(DSS_stamps.cor[11:20]) # 0.671
mean(DSS_stamps.cor[1:20]) # 0.545

# calculate difference between HIP and DSS in common years
HIP_DSS <- matrix(NA,nrow=n.states,ncol=4)
for (i in 1:n.states){
  for (j in 1:3){
    HIP_DSS[i,j] <- HIP2[i,(j+15)]/DSS2[i,(j+15)]}
  HIP_DSS[i,4] <- mean(HIP_DSS[i,1:3])
  }
plot(HIP_DSS[,4], ylab="HIP/DSS ratio", xlab="State", pch=19)
abline(h=1)
# HIP > DSS in most states, but not ubiquitous

# difference between HIP and DSS in overlap years (1999-2001)
mean(HIP2[,16:18])/mean(DSS2[,16:18]) # overall 1.86 higher under HIP
mean(HIP2[1:10,16:18])/mean(DSS2[1:10,16:18]) # central 2.03 higher in HIP
mean(HIP2[11:20,16:18])/mean(DSS2[11:20,16:18]) # eastern 1.54 higher under HIP
mean(HIP2[1,16:18])/mean(DSS2[1,16:18]) # LA, 1.75
mean(HIP2[2,16:18])/mean(DSS2[2,16:18]) # MI, 2.40
mean(HIP2[8,16:18])/mean(DSS2[8,16:18]) # IL/IN, 0.48

#-------------------------------------------------------------------------------------
# Fit splines to annual stamp sales to measure annual departures from long-term trends
# I'm assuming the spline indexes total hunters per state (waterfowl and woodcock)
#  whereas residuals indicate short term participation in waterfowl hunting per se 
# Define knots at a specified number of sample quantiles of the covariate
#-------------------------------------------------------------------------------------

years <- seq(2005,1964,-1) # start in 2005, poor stamp data in late 2000's
X <- years
Y1 <- Y2 <- Y3 <- Y4 <- Y5 <- Y6 <- Y7 <- Y8 <- Y9 <- Y10 <- numeric()
Y11 <- Y12 <- Y13 <- Y14 <- Y15 <- Y16 <- Y17 <- Y18 <- Y19 <- Y20 <- numeric()

# create vectors of response variables
Y1 <- log(stamps2[1,12:53]); Y11 <- log(stamps2[11,12:53])
Y2 <- log(stamps2[2,12:53]); Y12 <- log(stamps2[12,12:53])
Y3 <- log(stamps2[3,12:53]); Y13 <- log(stamps2[13,12:53])
Y4 <- log(stamps2[4,12:53]); Y14 <- log(stamps2[14,12:53])
Y5 <- log(stamps2[5,12:53]); Y15 <- log(stamps2[15,12:53])
Y6 <- log(stamps2[6,12:53]); Y16 <- log(stamps2[16,12:53])
Y7 <- log(stamps2[7,12:53]); Y17 <- log(stamps2[17,12:53])
Y8 <- log(stamps2[8,12:53]); Y18 <- log(stamps2[18,12:53])
Y9 <- log(stamps2[9,12:53]); Y19 <- log(stamps2[19,12:53])
Y10 <- log(stamps2[10,12:53]); Y20 <- log(stamps2[20,12:53])

# assess appropriate model complexity by altering df
Mich3 <- smooth.spline(X,Y2,df=3)
Mich4 <- smooth.spline(X,Y2,df=4)
Mich6 <- smooth.spline(X,Y2,df=6)
Mich10 <- smooth.spline(X,Y2,df=10)

# I'm using df=4 because it seems to fit appropriately, but similar results w/ 3 or 6
df=4

# fit splines to each state or state group
s1 <- smooth.spline(X,Y1,df=df); s11 <- smooth.spline(X,Y11,df=df)
s2 <- smooth.spline(X,Y2,df=df); s12 <- smooth.spline(X,Y12,df=df)
s3 <- smooth.spline(X,Y3,df=df); s13 <- smooth.spline(X,Y13,df=df)
s4 <- smooth.spline(X,Y4,df=df); s14 <- smooth.spline(X,Y14,df=df)
s5 <- smooth.spline(X,Y5,df=df); s15 <- smooth.spline(X,Y15,df=df)
s6 <- smooth.spline(X,Y6,df=df); s16 <- smooth.spline(X,Y16,df=df)
s7 <- smooth.spline(X,Y7,df=df); s17 <- smooth.spline(X,Y17,df=df)
s8 <- smooth.spline(X,Y8,df=df); s18 <- smooth.spline(X,Y18,df=df)
s9 <- smooth.spline(X,Y9,df=df); s19 <- smooth.spline(X,Y19,df=df)
s10 <- smooth.spline(X,Y10,df=df); s20 <- smooth.spline(X,Y20,df=df)

# plot for Fig. 2 in manuscript
par(mfrow=c(3,2))
plot(X,Y2,main="",xlab="",ylab="log Stamps",ylim=c(10.5,12.5)); lines(Mich3,col="black");lines(Mich4,col="red")
lines(Mich6,col="blue");lines(Mich10,col="green")
legend(1977,12.7,legend="Michigan",bty="n")
plot(X,Y4,main="",xlab="",ylab="",ylim=c(10.5,12.5)); lines(s4,col="red")
legend(1977,12.7,legend="Minnesota",bty="n")
plot(X,Y1,main="",xlab="",ylab="log Stamps",ylim=c(10,12)); lines(s1,col="red")
legend(1977,12.2,legend="Louisiana",bty="n")
plot(X,Y14,main="",xlab="",ylab="",ylim=c(10,12)); lines(s14,col="red")
legend(1977,12.2,legend="New York",bty="n")
plot(X,Y12,main="",xlab="",ylab="log Stamps",ylim=c(9,11)); lines(s12,col="red")
legend(1977,11.2,legend="Maine",bty="n")
plot(X,Y13,main="",xlab="",ylab="",ylim=c(9,11)); lines(s13,col="red")
legend(1977,11.2,legend="New Jersey",bty="n")

#plot stamps and fitted splines for all state groups
par(mfrow=c(2,1))
plot(X,Y1,main="LA"); lines(s1,col="red")
plot(X,Y2,main="MI"); lines(s2,col="red")
plot(X,Y3,main="WI"); lines(s3,col="red")
plot(X,Y4,main="MN"); lines(s4,col="red")
plot(X,Y5,main="OH"); lines(s5,col="red")
plot(X,Y6,main="IL/IN"); lines(s6,col="red")
plot(X,Y7,main="IA/MO/AR"); lines(s7,col="red")
plot(X,Y8,main="KY/TN"); lines(s8,col="red")
plot(X,Y9,main="MS/AL"); lines(s9,col="red")
plot(X,Y10,main="NE/KS/OK/TX"); lines(s10,col="red")
plot(X,Y11,main="NY"); lines(s11,col="red")
plot(X,Y12,main="ME"); lines(s12,col="red")
plot(X,Y13,main="NJ"); lines(s13,col="red")
plot(X,Y14,main="PA"); lines(s14,col="red")
plot(X,Y15,main="MA"); lines(s15,col="red")
plot(X,Y16,main="VT/NH"); lines(s16,col="red")
plot(X,Y17,main="CT/RI"); lines(s17,col="red")
plot(X,Y18,main="DE/MD/VA/WV"); lines(s18,col="red")
plot(X,Y19,main="NC/SC"); lines(s19,col="red")
plot(X,Y20,main="GA/FL"); lines(s20,col="red")

stamp.residuals <- matrix(NA,nrow=n.states,ncol=n.yrs)
  stamp.residuals[1,12:n.yrs]<-residuals(s1)
  stamp.residuals[2,12:n.yrs]<-residuals(s2)
  stamp.residuals[3,12:n.yrs]<-residuals(s3)
  stamp.residuals[4,12:n.yrs]<-residuals(s4)
  stamp.residuals[5,12:n.yrs]<-residuals(s5)
  stamp.residuals[6,12:n.yrs]<-residuals(s6)
  stamp.residuals[7,12:n.yrs]<-residuals(s7)
  stamp.residuals[8,12:n.yrs]<-residuals(s8)
  stamp.residuals[9,12:n.yrs]<-residuals(s9)
  stamp.residuals[10,12:n.yrs]<-residuals(s10)
  stamp.residuals[11,12:n.yrs]<-residuals(s11)
  stamp.residuals[12,12:n.yrs]<-residuals(s12)
  stamp.residuals[13,12:n.yrs]<-residuals(s13)
  stamp.residuals[14,12:n.yrs]<-residuals(s14)
  stamp.residuals[15,12:n.yrs]<-residuals(s15)
  stamp.residuals[16,12:n.yrs]<-residuals(s16)
  stamp.residuals[17,12:n.yrs]<-residuals(s17)
  stamp.residuals[18,12:n.yrs]<-residuals(s18)
  stamp.residuals[19,12:n.yrs]<-residuals(s19)
  stamp.residuals[20,12:n.yrs]<-residuals(s20)
  max(stamp.residuals,na.rm=TRUE) # 1.58

# calculate correlation between residuals and DSS
  resid_DSS <- numeric()
  for (i in 1:20){
    resid_DSS[i] <- cor(stamp.residuals[i,16:53],DSS2[i,16:53])
  }
print(round(resid_DSS,3))
# [1] 0.324 0.407 0.478 0.090 0.125 0.217 0.320 0.388 0.610 0.356 0.425 0.420
# [13] 0.529 0.499 0.220 0.554 0.366 0.481 0.348 0.127
  mean(resid_DSS[1:10]) # 0.331 Central
  mean(resid_DSS[11:20]) # 0.397 Eastern
  mean(resid_DSS[1:20]) # 0.364 All
  
# calculate mean residuals 1999-2001 (index of waterfowl hunting during survey overlap)
mean.resid <- max.resid <- min.resid <- fraction <-numeric()
for (i in 1:n.states){
  mean.resid[i] <- mean(stamp.residuals[i,16:18])
  min.resid[i] <- min(stamp.residuals[i,],na.rm=TRUE)
  max.resid[i] <- max(stamp.residuals[i,],na.rm=TRUE)
  fraction[i] <- (mean.resid[i]-min.resid[i])/(max.resid[i]-min.resid[i])  
}
print(round(fraction,3))
# [1] 0.947 0.466 0.582 0.653 0.699 0.637 0.860 0.816 0.239 0.898 0.588 0.887 0.390 0.528 0.177 0.642
# [17] 0.357 0.540 0.660 0.430
mean(mean.resid[1:20]) # 0.045 (slightly more duck stamps sold in overlap years)
mean(fraction[1:10]) # 0.680 fraction of maximum residuals
mean(fraction[11:20]) # 0.520
mean(fraction[1:20]) # 0.600

# use mean and sd from last 5 years HIP (2012-2016) for initial priors (reverse time)
# for each state, calculate residuals as fraction of max - min
mean5.HIP2 <- sd5.HIP2 <- max.stamps <- min.stamps <- numeric()
rel.stamps <- matrix(NA,nrow=n.states,ncol=n.yrs)
for (i in 1:n.states){
  mean5.HIP2[i] <- mean(HIP2[i,1:5],na.rm=TRUE)
  sd5.HIP2[i] <- sd(HIP2[i,1:5],na.rm=TRUE)
  max.stamps[i] <- max(stamp.residuals[i,],na.rm=TRUE)
  min.stamps[i] <- min(stamp.residuals[i,],na.rm=TRUE)
  for (j in 6:n.yrs){
    rel.stamps[i,j] <- (stamp.residuals[i,j]-min.stamps[i])/(max.stamps[i]-min.stamps[i])
  }}

######################################################
# specify jags model to predict total annual harvest
# Three components to harvest, 
# 1) only hunt woodcock, never ducks (HIP only, never DSS)
# 2) always hunt ducks too, always captured in DSS
# 3) sometimes hunt ducks, occasionally captured in DSS when sales high relative to long-term trend
######################################################

sink("AMWO.state.harvest.bug")
cat("
  model {
  # specify means and variances for state-level effects
    mean.r ~ dunif(-1,1) # variables to describe mean trend in harvest
    sd.r ~ dunif(0.01,1) # state-specific variation in mean trend
    tau.r <- pow(sd.r,-2)
    mu.sigmaH ~ dunif(0.01,10) # prior for mean variance
    sd.sigmaH ~ dunif(0.01,10)
    tau.sigmaH <- pow(sd.sigmaH,-2)
  for (i in 1:n.states){
    # priors for state-specific 2016 harvest (mean HIP estimate from last 5 years, converted to log scale)
    Harvest[i,1] ~ dnorm(log.mean5.HIP[i],tau.mean5.HIP[i])
    # priors for annual trend in harvest growth rate (reverse time)
    trend[i] ~ dnorm(mean.r,tau.r) 
    sigma.H[i] ~ dlnorm(mu.sigmaH,tau.sigmaH)
    tau.H[i] <- pow(sigma.H[i],-2)
    
    # prior for state-specific proportion of AMWO harvest by hunters who never, always, or sometimes hunt waterfowl
    # note these are treated as fixe effects, but could be coded as random effects using logit and mlogit links
    pi.1[i] ~ dunif(0,0.7)
    pi.2[i] ~ dunif(0,1-pi.1[i])
    pi.3[i] <- max(0,(1 - pi.1[i] - pi.2[i]))  # don't allow pi.3 to be negative
  } # end i (state) loop
    
  # State process (estimate true Harvest in years 2-54)
  for (i in 1:n.states){
    for (j in 2:n.yrs){
      eps[i,j] ~ dnorm(trend[i],tau.H[i]) 
      Harvest[i,j] <- Harvest[i,j-1] + eps[i,j] # autoregressive model, don't allow negative harvest
    } # end j (year) loop
    } # end i (state) loop
    
    # observation process, recoveries and harvest data
    for (i in 1:n.states){
    for (j in 1:n.HIP){
      log.HIP[i,j] ~ dnorm(Harvest[i,j],tau.HIP[i,j])} # end HIP years, measures all harvest
    for (j in 16:n.yrs){
    # consider adding state-specific multiplier to H3
      log.DSS[i,j] ~ dnorm(log(pi.2[i]*exp(Harvest[i,j])+frac[i,j]*pi.3[i]*exp(Harvest[i,j])),tau.DSS[i,j])} # end DSS yrs
    } # end i loop
    
    } # end jags model
    ",fill = TRUE)
sink()

# Bundle data, set inits, parameters monitored
jags.data <- list(n.yrs=n.yrs, n.states=n.states, n.HIP=n.HIP, 
                  log.mean5.HIP=log(mean5.HIP2+1), 
                  tau.mean5.HIP=1/(sd5.HIP2/mean5.HIP2^2)^2,
                  log.HIP=log(HIP2+1), tau.HIP=(HIP_var2/HIP2^2)^-1, 
                  log.DSS=log(DSS2+1), tau.DSS=(DSS_var2/DSS2^2)^-1, 
                  frac=rel.stamps)
inits <- NULL
parameters <- c("pi.1","pi.2","pi.3","trend","sigma.H","Harvest")

# MCMC settings (110000 ni, 10000 nb leads to full convergence, R-hat < 1.03, ~106 min)
# 11000, 1000 took 12 min
na <- 1000; ni <- 60000; nt <- 10
nb <- 10000; nc <- 3
n.iter <- ((ni-nb)/nt)*nc # calculate number of retained iterations

# call jags
AMWO.Harvest <- jagsUI(jags.data, inits, parameters, "AMWO.state.harvest.bug", n.adapt = na, n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, 
                  parallel=TRUE)

print(AMWO.Harvest, digits=2)

max(max.rhat<-c(sapply(AMWO.Harvest$Rhat,max,na.rm=TRUE))) # report max rhat of all variables

save(AMWO.Harvest, file="AMWO_Harvest.Rdata")

# plot posteriors for state-specific harvest trends and annual variation in harvest
par(mfrow=c(2,1))
plot(AMWO.Harvest$mean$trend, ylim=c(0,0.08), pch=19, ylab="Annual rate of decline in harvest", xlab="State #")
segments(1:20,AMWO.Harvest$mean$trend-1.96*AMWO.Harvest$sd$trend,
         1:20,AMWO.Harvest$mean$trend+1.96*AMWO.Harvest$sd$trend,lty=2)
plot(AMWO.Harvest$mean$sigma.H, pch=19, ylab="Annual variation in harvest", xlab="State #", ylim=c(0,1))
  segments(1:20,AMWO.Harvest$mean$sigma.H-1.96*AMWO.Harvest$sd$sigma.H,
           1:20,AMWO.Harvest$mean$sigma.H+1.96*AMWO.Harvest$sd$sigma.H,lty=2)

# plot posteriors for estimated proportions of hunters by state who never, always, or sometimes hunt waterfowl
par(mfrow=c(3,1))
plot(AMWO.Harvest$mean$pi.1, pch=19, ylab="Proportion never hunting ducks", xlab="State #", ylim=c(0,1),
     main="Never hunt waterfowl")
  segments(1:20,AMWO.Harvest$mean$pi.1-1.96*AMWO.Harvest$sd$pi.1,
           1:20,AMWO.Harvest$mean$pi.1+1.96*AMWO.Harvest$sd$pi.1,lty=2)
plot(AMWO.Harvest$mean$pi.2, pch=19, ylab="Proportion never hunting ducks", xlab="State #", ylim=c(0,1),
     main="Always hunt waterfowl")
  segments(1:20,AMWO.Harvest$mean$pi.2-1.96*AMWO.Harvest$sd$pi.2,
           1:20,AMWO.Harvest$mean$pi.2+1.96*AMWO.Harvest$sd$pi.2,lty=2)
plot(AMWO.Harvest$mean$pi.3, pch=19, ylab="Proportion never hunting ducks", xlab="State #", ylim=c(0,1),
     main="Sometimes hunt waterfowl")
  segments(1:20,AMWO.Harvest$mean$pi.3-1.96*AMWO.Harvest$sd$pi.3,
           1:20,AMWO.Harvest$mean$pi.3+1.96*AMWO.Harvest$sd$pi.3,lty=2)
  
round(AMWO.Harvest$mean$pi.1,3) # never hunt ducks (labeled group 3 in paper)

# calculate means of state averages (not weighted by harvest)
  mean(AMWO.Harvest$mean$pi.1[1:10]) 
  mean(AMWO.Harvest$mean$pi.1[11:20])
  mean(AMWO.Harvest$mean$pi.1[1:20])
  
round(AMWO.Harvest$mean$pi.2,3) # always hunt waterfowl too
  mean(AMWO.Harvest$mean$pi.2[1:10]) 
  mean(AMWO.Harvest$mean$pi.2[11:20]) 
  mean(AMWO.Harvest$mean$pi.2[1:20])

round(AMWO.Harvest$mean$pi.3,3) # sometimes hunt waterfowl
  mean(AMWO.Harvest$mean$pi.3[1:10]) 
  mean(AMWO.Harvest$mean$pi.3[11:20])
  mean(AMWO.Harvest$mean$pi.3[1:20])
  
state.id <- c("LA","MI","WI","MN","OH","IL","MO","KY","MS","TX","NY","ME","NJ","PA","MA","VT","CT","MD","NC","GA")

# calculate mean harvest during all years (mean.H) and HIP:DSS overlap years (mean.overlap.H)

mean.H <- mean.overlap.H <- numeric()
for (i in 1:20){
mean.H[i]<-mean(exp(AMWO.Harvest$mean$Harvest[i,1:53]))
mean.overlap.H[i]<-mean(exp(AMWO.Harvest$mean$Harvest[i,16:18]))
}
print(round(mean.H,0))
print(round(mean.overlap.H,0))

sum(mean.H[1:10]) 
sum(mean.H[11:20])
sum(mean.H[1:20])

sum(mean.overlap.H[1:10]) 
sum(mean.overlap.H[11:20])
sum(mean.overlap.H[1:20])


# create plot parameters (here I'm plotting 90% credible intervals)

lower5 <- upper5 <- matrix(NA, nrow=n.states+3, ncol=n.yrs)
E.mean5 <- C.mean5 <- C.mean5.sansLA <- numeric()
for (s in 1:n.states){
  for (t in 1:n.yrs){
    lower5[s,t] <- exp(quantile(AMWO.Harvest$sims.list$Harvest[,s,t], 0.05))
    upper5[s,t] <- exp(quantile(AMWO.Harvest$sims.list$Harvest[,s,t], 0.95))}}

min(exp(AMWO.Harvest$sims.list$Harvest[,,])) # 97 avoids 0 harvest

# summarize eastern and central regions
eastern5 <- central5.sansLA <- central5 <- matrix(NA, nrow=n.iter, ncol=n.yrs)
E.stamps <- C.stamps <- E.DSS <- C.DSS <- E.HIP <- C.HIP <- numeric() 
  E.HIP.sd <- C.HIP.sd <- E.DSS.sd <- C.DSS.sd <- numeric()
for (i in 1:n.iter){
  for (t in 1:n.yrs){
    central5[i,t] <- sum(exp(AMWO.Harvest$sims.list$Harvest[i,1:10,t]))
    central5.sansLA[i,t] <- sum(exp(AMWO.Harvest$sims.list$Harvest[i,2:10,t]))
    eastern5[i,t] <- sum(exp(AMWO.Harvest$sims.list$Harvest[i,11:20,t]))
  }}
for (t in 1:n.yrs){
  E.mean5[t] <- mean(eastern5[,t])
  C.mean5[t] <- mean(central5[,t])
  C.mean5.sansLA[t] <- mean(central5.sansLA[,t])
  lower5[(n.states+1),t] <- quantile(eastern5[,t],0.05)  
  upper5[(n.states+1),t] <- quantile(eastern5[,t],0.95)  
  lower5[(n.states+2),t] <- quantile(central5[,t],0.05)  
  upper5[(n.states+2),t] <- quantile(central5[,t],0.95)
  lower5[(n.states+3),t] <- quantile(central5.sansLA[,t],0.05)  
  upper5[(n.states+3),t] <- quantile(central5.sansLA[,t],0.95)
  E.stamps[t] <- sum(stamps2[11:20,t])
  C.stamps[t] <- sum(stamps2[1:10,t])
  E.DSS[t] <- sum(DSS2[11:20,t])
  C.DSS[t] <- sum(DSS2[1:10,t])
  E.HIP[t] <- sum(HIP2[11:20,t])
  C.HIP[t] <- sum(HIP2[1:10,t])
  E.HIP.sd[t] <- sqrt(sum(HIP_var2[11:20,t]))
  C.HIP.sd[t] <- sqrt(sum(HIP_var2[1:10,t]))
  E.DSS.sd[t] <- sqrt(sum(DSS_var2[11:20,t]))
  C.DSS.sd[t] <- sqrt(sum(DSS_var2[1:10,t]))
}  

# correlations of DSS harvest estimates w stamp sales, mgmt unit scale and total US
cor(E.stamps[16:53],E.DSS[16:53]) 
cor(C.stamps[16:53],C.DSS[16:53])
cor(C.stamps[16:53]+E.stamps[16:53],C.DSS[16:53]+E.DSS[16:53])

# weighted harvest components in years 16-18 (1999-2001)
C.total.1 <- E.total.1 <- C.total.2 <- E.total.2 <- C.total.3 <- E.total.3 <- 0
for (i in 1:10){
  C.total.1 <- C.total.1 + AMWO.Harvest$mean$pi.1[i]*AMWO.Harvest$mean$Harvest[i,16:18] 
  E.total.1 <- E.total.1 + AMWO.Harvest$mean$pi.1[(i+10)]*AMWO.Harvest$mean$Harvest[(i+10),16:18] 
  C.total.2 <- C.total.2 + AMWO.Harvest$mean$pi.2[i]*AMWO.Harvest$mean$Harvest[i,16:18] 
  E.total.2 <- E.total.2 + AMWO.Harvest$mean$pi.2[(i+10)]*AMWO.Harvest$mean$Harvest[(i+10),16:18] 
  C.total.3 <- C.total.3 + AMWO.Harvest$mean$pi.3[i]*AMWO.Harvest$mean$Harvest[i,16:18] 
  E.total.3 <- E.total.3 + AMWO.Harvest$mean$pi.3[(i+10)]*AMWO.Harvest$mean$Harvest[(i+10),16:18] 
}
mean(C.total.1/(C.total.1+C.total.2+C.total.3)) # 0.217 proportion Central AMWO only
mean(C.total.2/(C.total.1+C.total.2+C.total.3)) # 0.492 proportion Central ducks always
mean(C.total.3/(C.total.1+C.total.2+C.total.3)) # 0.291 proportion Central ducks sometimes
mean(E.total.1/(E.total.1+E.total.2+E.total.3)) # 0.152 proportion Eastern AMWO only
mean(E.total.2/(E.total.1+E.total.2+E.total.3)) # 0.612 proportion Eastern ducks always
mean(E.total.3/(E.total.1+E.total.2+E.total.3)) # 0.236 proportion Eastern ducks sometimes
# combined all states
mean((C.total.1+E.total.1)/(C.total.1+C.total.2+C.total.3+E.total.1+E.total.2+E.total.3)) # 0.185
mean((C.total.2+E.total.2)/(C.total.1+C.total.2+C.total.3+E.total.1+E.total.2+E.total.3)) # 0.551
mean((C.total.3+E.total.3)/(C.total.1+C.total.2+C.total.3+E.total.1+E.total.2+E.total.3)) # 0.264

# plot harvest estimates

year <- 1964:2016
par(mfrow=c(1,1))
# plot Central harvest
plot(year+0.1,rev(C.mean5/1000),pch=19,main="",ylab="Harvest (1000's)",
     xlab="",type="b",ylim=c(0,1150),cex.lab=1.5,cex.axis=1.2)
segments(year+0.1, rev(lower5[22,]/1000), year+0.1, rev(upper5[22,]/1000))
points(year,rev(C.DSS/1000),pch=21,col="red")
segments(year,rev((C.DSS-1.28*C.DSS.sd)/1000),
         year,rev((C.DSS+1.28*C.DSS.sd)/1000),col="red")
points(year+0.2,rev(C.HIP/1000),pch=21,col="blue")
segments(year+0.2,rev((C.HIP-1.28*C.HIP.sd)/1000),
         year+0.2,rev((C.HIP+1.28*C.HIP.sd)/1000),col="blue")
legend(1980,1160,legend=c("Central AMWO Harvest"),cex=1.3,bty="n")
legend(1983,1070,legend=c("State-space model (1964-2016)",
                          "Duck stamp survey (1964-2001)", 
                          "Harvest Information Program (1999-2016)"),
       bty="n",text.col=c("black","red","blue"))

# plot Central harvest (with and w/out Louisiana)
 # is peak in late 1970's caused by high harvest estimates in LA? Somewhat, but not exclusively
plot(year+0.1,rev(C.mean5/1000),pch=19,ylim=c(0,1150),main="",ylab="Harvest (thousands)",xlab="",type="b")
segments(year+0.1, rev(lower5[22,]/1000), year+0.1, rev(upper5[22,]/1000))
points(year+0.1,rev(C.mean5.sansLA/1000),pch=19,type="b",col="red")
segments(year+0.1, rev(lower5[23,]/1000), year+0.1, rev(upper5[23,]/1000),col="red")
legend(1990,900,legend=c("All Central Population states","Louisiana excluded"),bty="n",text.col=c("black","red"))

# plot Eastern harvest
par(mfrow=c(1,1)) # one plot per panel
plot(year+0.1,rev(E.mean5/1000),pch=19,ylim=c(0,550),main="",
     ylab="Harvest (1000's)",xlab="",type="b",cex.lab=1.5,cex.axis=1.2)
segments(year+0.1, rev(lower5[21,]/1000), year+0.1, rev(upper5[21,]/1000))
points(year,rev(E.DSS/1000),pch=21,col="red")
segments(year,rev((E.DSS-1.28*E.DSS.sd)/1000),
         year,rev((E.DSS+1.28*E.DSS.sd)/1000),col="red")
points(year+0.2,rev(E.HIP/1000),pch=21,col="blue")
segments(year+0.2,rev((E.HIP-1.28*E.HIP.sd)/1000),
         year+0.2,rev((E.HIP+1.28*E.HIP.sd)/1000),col="blue")
legend(1980,570,legend=c("Eastern AMWO Harvest"),cex=1.3,bty="n")

# Plot MI harvest (state 2)
par(mfrow=c(1,1))
# data from MI independent survey of AMWO harvest
MI.AMWO <- c(173361,154918,177049,177049,165984,140164,154918,195492,188115,
             250820,280328,387295,387295,354098,309836,324590,302459,250820,261885,
             206557,NA,239754,228689,250820,254508,295082,239754,276639,165984,
             195492,202869,236066,228689,158607,191803,177049,143852,154918,129098,
             140164,118033,121721,136475,110656,121721,92213,99590,95902)
plot(1964:2016,rev(HIP2[2,1:53]/1000),ylim=c(0,400),
     ylab="Harvest (1000's)",xlab="",pch=21,col="blue",main="",cex.lab=1.5,cex.axis=1.2)
segments(1964:2016, rev((HIP2[2,]-1.28*HIP_sd2[2,])/1000), 1964:2016, rev((HIP2[2,]+1.28*HIP_sd2[2,])/1000),col="blue")
points(1964:2001, rev(DSS2[2,16:53]/1000), pch=21, col="red")
segments(1964:2001, rev((DSS2[2,]-1.28*DSS_sd2[2,])/1000), 1964:2001, rev((DSS2[2,]+1.28*DSS_sd2[2,])/1000),col="red")
# add MI state harvest estimate
points(1964:2011,MI.AMWO/1000,col="green",pch=19)
# add Bayesian estiamtes (state-specific slopes)
points(year-0.1,rev(exp(AMWO.Harvest$mean$Harvest[2,1:53])/1000),pch=19,col="black",type="b")
segments((1964:2016)-0.1, rev(lower5[2,]/1000), (1964:2016)-0.1, rev(upper5[2,]/1000),col="black")
legend(1982,420,legend=c("Michigan"),bty="n",cex=1.8)
# underestimates relative to MI state survey in late 1980s and 90s, but does well in 60s through early 80s

par(mfrow=c(1,1))
# Plot MN harvest (2)
# data for independent AMWO harvest in Minnesota (MN DNR: 1993-2015)
MN.AMWO2 <- c(68000, 74000, 82000, 58000, 58000, 63000, 54000, 45341, 26662, 28230, 30438, 41479, 27919, 39907, 27866, 29210, 35430, 29770, 24980, 30360, 31920, 25810, 37270) 
plot(1964:2016,rev(HIP2[4,1:53]/1000),ylim=c(0,250),
     ylab="Harvest (1000's)",xlab="",pch=21,col="blue",main="",cex.lab=1.5,cex.axis=1.2)
segments(1964:2016, rev((HIP2[4,]-1.28*HIP_sd2[4,])/1000), 1964:2016, rev((HIP2[4,]+1.28*HIP_sd2[4,])/1000),col="blue")
points(1964:2001, rev(DSS2[4,16:53]/1000), pch=21, col="red")
segments(1964:2001, rev((DSS2[4,]-1.28*DSS_sd2[4,])/1000), 1964:2001, rev((DSS2[4,]+1.28*DSS_sd2[4,])/1000),col="red")
# add MI state harvest estimate
points(1993:2015,MN.AMWO2/1000,col="green",pch=19) # MN DNR independent survey of AMWO harvest
# add Bayesian estiamtes (state-specific slopes)
points(year-0.1,rev(exp(AMWO.Harvest$mean$Harvest[4,1:53])/1000),pch=19,col="black",type="b")
segments((1964:2016)-0.1, rev(lower5[4,]/1000), (1964:2016)-0.1, rev(upper5[4,]/1000),col="black")
legend(1985,250,legend=c("Minnesota"),bty="n",cex=1.8)

length(MI.AMWO) #48
53-16+1 # 38 years overlap with DSS
cor.test(MI.AMWO[1:48],rev(AMWO.Harvest$mean$Harvest[2,6:53]),use="complete.obs") # r = 0.78, MI DNR vs our estimates
cor.test(MN.AMWO2[1:23],rev(AMWO.Harvest$mean$Harvest[4,2:24]),use="complete.obs") # r = 0.73, MN DNR vs our estimates

# Plot LA harvest (not very different from DSS)
par(mfrow=c(1,1))
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[1,1:53])/1000),ylim=c(0,550),pch=19,type="b",main="",ylab="Harvest (1000,s)",cex.lab=1.5,cex.axis=1.2,xlab="")
segments((1964:2016), rev(lower5[1,]/1000), (1964:2016), rev(upper5[1,]/1000),lty=1)
points(1964:2001,rev(DSS2[1,16:53]/1000),col="red")
segments(1964:2001, rev((DSS2[1,]-1.28*DSS_sd2[1,])/1000), 1964:2001, rev((DSS2[1,]+1.28*DSS_sd2[1,])/1000),col="red")
points(1964:2016,rev(HIP2[1,1:53]/1000),col="blue")
segments(1964:2016, rev((HIP2[1,]-1.28*HIP_sd2[1,])/1000), 1964:2016, rev((HIP2[1,]+1.28*HIP_sd2[1,])/1000),col="blue")
legend(1982,570,legend=c("Louisiana"),bty="n",cex=1.8)

# Maine harvest (DSS grossly underestimates due to large proportion of non-waterfowl hunters)
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[12,1:53])/1000),ylim=c(0,110),pch=19,type="b",main="",ylab="Harvest (1000's)",cex.lab=1.5,cex.axis=1.2,xlab="")
segments((1964:2016), rev(lower5[12,]/1000), (1964:2016), rev(upper5[12,]/1000),lty=1)
points(1964:2001,rev(DSS2[12,16:53]/1000),col="red")
segments(1964:2001, rev((DSS2[12,]-1.28*DSS_sd2[12,])/1000), 1964:2001, rev((DSS2[12,]+1.28*DSS_sd2[12,])/1000),col="red")
points(1964:2016,rev(HIP2[12,1:53]/1000),col="blue")
segments(1964:2016, rev((HIP2[12,]-1.28*HIP_sd2[12,])/1000), 1964:2016, rev((HIP2[12,]+1.28*HIP_sd2[12,])/1000),col="blue")
legend(1982,115,legend=c("Maine"),bty="n",cex=1.8)

# Plot NJ harvest (nearly identical to DSS estimates)
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[13,1:53])/1000),ylim=c(0,60),pch=19,type="b",main="",ylab="",cex.lab=1.5,cex.axis=1.2,xlab="")
segments((1964:2016), rev(lower5[13,]/1000), (1964:2016), rev(upper5[13,]/1000),lty=1)
points(1964:2001-0.1,rev(DSS2[13,16:53]/1000),col="red")
segments(1964:2001-0.1, rev((DSS2[13,]-1.28*DSS_sd2[13,])/1000), 1964:2001-0.1, rev((DSS2[13,]+1.28*DSS_sd2[13,])/1000),col="red")
points(1964:2016,rev(HIP2[13,1:53]/1000),col="blue")
segments(1964:2016, rev((HIP2[13,]-1.28*HIP_sd2[13,])/1000), 1964:2016, rev((HIP2[13,]+1.28*HIP_sd2[13,])/1000),col="blue")
legend(1980,60,legend=c("New Jersey"),bty="n",cex=1.8)

# Plot NY harvest
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[14,1:53])/1000),ylim=c(0,110),pch=19,type="b",main="",ylab="Harvest (1000's)",cex.lab=1.5,cex.axis=1.2,xlab="")
segments((1964:2016), rev(lower5[14,]/1000), (1964:2016), rev(upper5[14,]/1000),lty=1)
points(1964:2001-0.1,rev(DSS2[14,16:53]/1000),col="red")
segments(1964:2001-0.1, rev((DSS2[14,]-1.28*DSS_sd2[14,])/1000), 1964:2001-0.1, rev((DSS2[14,]+1.28*DSS_sd2[14,])/1000),col="red")
points(1964:2016,rev(HIP2[14,1:53]/1000),col="blue")
segments(1964:2016, rev((HIP2[14,]-1.28*HIP_sd2[14,])/1000), 1964:2016, rev((HIP2[14,]+1.28*HIP_sd2[14,])/1000),col="blue")
legend(1980,110,legend=c("New York"),bty="n",cex=1.8)

# Plots of additional states (including ones not included in manuscript)
# (WI, MN, IA/MO/AR, IL/IN)
par(mfrow=c(2,2)) # plot in 2x2 panels
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[3,1:53])),ylim=c(0,250000),pch=19,type="b",main="Wisconsin",ylab="Harvest")
segments((1964:2016), rev(lower5[3,]), (1964:2016), rev(upper5[3,]),lty=1)
points(1964:2001,rev(DSS2[3,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[3,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[4,1:53])),ylim=c(0,250000),pch=19,type="b",main="Minnesota",ylab="Harvest")
segments((1964:2016), rev(lower5[4,]), (1964:2016), rev(upper5[4,]),lty=1)
points(1964:2001,rev(DSS2[4,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[4,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[7,1:53])),ylim=c(0,50000),pch=19,type="b",main="IA/MO/AR",ylab="Harvest")
segments((1964:2016), rev(lower5[7,]), (1964:2016), rev(upper5[7,]),lty=1)
points(1964:2001,rev(DSS2[7,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[7,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[6,1:53])),ylim=c(0,50000),pch=19,type="b",main="Illinois/Indiana",ylab="Harvest")
segments((1964:2016), rev(lower5[6,]), (1964:2016), rev(upper5[6,]),lty=1)
points(1964:2001,rev(DSS2[6,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[6,1:53]),pch=19,col="blue")

# OH, KY/TN, MS/AL, NE/KS/OK/TX
par(mfrow=c(2,2))
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[5,1:53])),ylim=c(0,30000),pch=19,type="b",main="Ohio",ylab="Harvest")
segments((1964:2016), rev(lower5[5,]), (1964:2016), rev(upper5[5,]),lty=1)
points(1964:2001,rev(DSS2[5,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[5,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[8,1:53])),ylim=c(0,30000),pch=19,type="b",main="Kentucky/Tennessee",ylab="Harvest")
segments((1964:2016), rev(lower5[8,]), (1964:2016), rev(upper5[8,]),lty=1)
points(1964:2001,rev(DSS2[8,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[8,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[9,1:53])),ylim=c(0,50000),pch=19,type="b",main="MS/AL",ylab="Harvest")
segments((1964:2016), rev(lower5[9,]), (1964:2016), rev(upper5[9,]),lty=1)
points(1964:2001,rev(DSS2[9,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[9,1:53]),pch=19,col="blue")
# this is eastern edge of Central Flyway
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[10,1:53])),ylim=c(0,50000),pch=19,type="b",main="NE/KS/OK/TX",ylab="Harvest")
segments((1964:2016), rev(lower5[10,]), (1964:2016), rev(upper5[10,]),lty=1)
points(1964:2001,rev(DSS2[10,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[10,1:53]),pch=19,col="blue")

# plot 11-14 (NY, PA, MA, VT/NH)
par(mfrow=c(2,2))
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[11,1:53])),ylim=c(0,150000),pch=19,type="b",main="New York",ylab="Harvest")
segments((1964:2016), rev(lower5[11,]), (1964:2016), rev(upper5[11,]),lty=1)
points(1964:2001,rev(DSS2[11,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[11,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[14,1:53])),ylim=c(0,150000),pch=19,type="b",main="Pennsylvania",ylab="Harvest")
segments((1964:2016), rev(lower5[14,]), (1964:2016), rev(upper5[14,]),lty=1)
points(1964:2001,rev(DSS2[14,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[14,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[15,1:53])),ylim=c(0,50000),pch=19,type="b",main="Massachusetts",ylab="Harvest")
segments((1964:2016), rev(lower5[15,]), (1964:2016), rev(upper5[15,]),lty=1)
points(1964:2001,rev(DSS2[15,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[15,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[16,1:53])),ylim=c(0,50000),pch=19,type="b",main="VT/NH",ylab="Harvest")
segments((1964:2016), rev(lower5[16,]), (1964:2016), rev(upper5[16,]),lty=1)
points(1964:2001,rev(DSS2[16,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[16,1:53]),pch=19,col="blue")

# plot 15-18 (MA, VT/NH, CT/RI, DE/MD/VA/WV)
par(mfrow=c(2,2))
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[17,1:53])),ylim=c(0,50000),pch=19,type="b",main="CT/RI",ylab="Harvest")
segments((1964:2016), rev(lower5[17,]), (1964:2016), rev(upper5[17,]),lty=1)
points(1964:2001,rev(DSS2[17,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[17,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[18,1:53])),ylim=c(0,50000),pch=19,type="b",main="DE/MD/VA/WV",ylab="Harvest")
segments((1964:2016), rev(lower5[18,]), (1964:2016), rev(upper5[18,]),lty=1)
points(1964:2001,rev(DSS2[18,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[18,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[19,1:53])),ylim=c(0,40000),pch=19,type="b",main="N & S Carolina",ylab="Harvest")
segments((1964:2016), rev(lower5[19,]), (1964:2016), rev(upper5[19,]),lty=1)
points(1964:2001,rev(DSS2[19,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[19,1:53]),pch=19,col="blue")
plot(year,rev(exp(AMWO.Harvest$mean$Harvest[20,1:53])),ylim=c(0,40000),pch=19,type="b",main="Georgia/Florida",ylab="Harvest")
segments((1964:2016), rev(lower5[20,]), (1964:2016), rev(upper5[20,]),lty=1)
points(1964:2001,rev(DSS2[20,16:53]),pch=19,col="red")
points(1964:2016,rev(HIP2[20,1:53]),pch=19,col="blue")

# write model-based estimates and SD (equivalent to frequentist SE) to files
mean.harvest <- sd.harvest <- matrix(NA, nrow=n.yrs, ncol=24)

for (j in 1:n.yrs){
  for (i in 1:20){
    mean.harvest[j,i+1] <- round(exp(AMWO.Harvest$mean$Harvest[i,j]),1)
    # SD as one half upper plus one half lower (asymmetrical due to log-normal dist)
    sd.harvest[j,i+1] <- round(0.5*(exp(AMWO.Harvest$mean$Harvest[i,j]+AMWO.Harvest$sd$Harvest[i,j])-exp(AMWO.Harvest$mean$Harvest[i,j])) +
      0.5*(exp(AMWO.Harvest$mean$Harvest[i,j])-exp(AMWO.Harvest$mean$Harvest[i,j]-AMWO.Harvest$sd$Harvest[i,j])),1)
  }
  mean.harvest[j,22] <- sum(mean.harvest[j,2:11]) # central
  mean.harvest[j,23] <- sum(mean.harvest[j,12:21]) # eastern
  mean.harvest[j,24] <- sum(mean.harvest[j,22:23]) # total
  sd.harvest[j,22] <- sqrt(sd.harvest[j,2]^2 + sd.harvest[j,3]^2 + sd.harvest[j,4]^2 + sd.harvest[j,5]^2 + 
                               sd.harvest[j,6]^2 + sd.harvest[j,7]^2 + sd.harvest[j,8]^2 + sd.harvest[j,9]^2 + 
                               sd.harvest[j,10]^2 + sd.harvest[j,11]^2)
  sd.harvest[j,23] <- sqrt(sd.harvest[j,12]^2 + sd.harvest[j,13]^2 + sd.harvest[j,14]^2 + sd.harvest[j,15]^2 + 
                               sd.harvest[j,16]^2 + sd.harvest[j,17]^2 + sd.harvest[j,18]^2 + sd.harvest[j,19]^2 + 
                               sd.harvest[j,20]^2 + sd.harvest[j,21]^2)
  sd.harvest[j,24] <- sqrt(sd.harvest[j,22]^2 + sd.harvest[j,23]^2)
}

colnames(mean.harvest) <- c("Year","LA","MI","WI","MN","OH","IL-IN","IA-MO-AR","KY-TN","MS-AL","TX-OK-KS-NE",
                         "NY","ME","NJ","PA","MA","VT-NH", "CT-RI","MD-DE-VA-WV","NC-SC","GA-FL",
                         "Central","Eastern","Total")
colnames(sd.harvest) <- c("Year","LA","MI","WI","MN","OH","IL-IN","IA-MO-AR","KY-TN","MS-AL","TX-OK-KS-NE",
                      "NY","ME","NJ","PA","MA","VT-NH", "CT-RI","MD-DE-VA-WV","NC-SC","GA-FL",
                      "Central","Eastern","Total")
mean.harvest[1:n.yrs,1] <- seq(2016,1964,-1)
sd.harvest[1:n.yrs,1] <- seq(2016,1964,-1)

write.csv(mean.harvest,file="Arnold_harvest_estimates.csv")
write.csv(sd.harvest,file="Arnold_harvest_estimates_SD.csv")