# Meta-analysis of band-reporting rates in North American waterfowl
# T. Arnold, Univ of Minnesota, arnol065@umn.edu
# 31 Dec 2018, added data from Anderson & Henny 1972
# 24 Jun 2019, added variable to exclude pseudo-replicates of pooled year studies

library(jagsUI)
library(ggplot2)
reporting <- read.csv("C:/Users/arnol065/Documents/projects/submitted/band_reporting_meta-analysis/data/reporting_rate_summary.csv", header=TRUE)
summary(reporting) # 421 observations of 38 variables
# subset to remove pseudoreplicates from pooled year studies
reporting <- subset(reporting,reporting$xtra_yr!=1) # 337 observations

reporting$survey3 <- reporting$survey + 2*reporting$survey2
reporting$survey3 <- as.factor(reporting$survey3)

# survey3: 0 = reward band, 1 = harvest survey, 2 = Padding & Royle methodology
table(reporting$survey3)
# 337 reporting rate estimates
table(reporting$spp)
table(reporting$r_yr)

sum(reporting$survey) # Martinson methodology 141
sum(reporting$survey2) # Padding & Royle 50
sum(reporting$reward) # reward band methodology 142 (4 more from Dubovsky, unpubl.)

plot(reporting$rho~reporting$r_yr, ylim=c(0,1), ylab="Reporting probability (p)",xlab="", pch=19, cex=0.6)
segments(reporting$r_yr,reporting$rho-reporting$se_rho,reporting$r_yr,reporting$rho+reporting$se_rho)

# Fig 1 w/ all symbols identical
rho <- ggplot(reporting, aes(x=r_yr, y=rho)) + geom_point() +
  geom_smooth(mapping=aes(x=r_yr, y=rho), method='loess', span=0.3, inherit.aes = FALSE)
rho + theme_classic() + coord_cartesian(ylim=c(0,1)) + 
  labs(x = "Year", y = "Reporting probability (p)") + 
  theme(axis.title.y = element_text(size=14)) +
  theme(axis.title.x = element_text(size=14)) +
  theme(axis.text.y = element_text(size=12)) +
  theme(axis.text.x = element_text(size=12)) +
  scale_x_continuous(breaks=c(1950,1960,1970,1980,1990,2000,2010))

# Fig 1b w/ color coded symbols (duplicates from multi-year estimates removed)
rho <- ggplot(reporting, aes(x=r_yr, y=rho, colour=survey3)) + geom_point() +
  geom_smooth(mapping=aes(x=r_yr, y=rho), method='loess', span=0.3, inherit.aes = FALSE)
rho + theme_classic() + coord_cartesian(ylim=c(0,1)) + 
  labs(x = "", y = "Reporting probability (p)") + 
  theme(axis.title.y = element_text(size=14)) +
  theme(axis.text.y = element_text(size=12)) +
  theme(axis.text.x = element_text(size=12)) +
  scale_x_continuous(breaks=c(1950,1960,1970,1980,1990,2000,2010))

# Fig 1c w/ shape coded symbols (duplicates from multi-year estimates removed)
rho <- ggplot(reporting, aes(x=r_yr, y=rho, shape=survey3)) + geom_point() +
  geom_smooth(mapping=aes(x=r_yr, y=rho), method='loess', span=0.3, inherit.aes = FALSE)
rho + theme_classic() + coord_cartesian(ylim=c(0,1)) + 
  scale_shape_manual(values=c(19,3,2)) +
  labs(x = "Year", y = "Reporting probability (p)") + 
  theme(axis.title.y = element_text(size=14)) +
  theme(axis.title.x = element_text(size=14)) +
  theme(axis.text.y = element_text(size=12)) +
  theme(axis.text.x = element_text(size=12)) +
  scale_x_continuous(breaks=c(1950,1960,1970,1980,1990,2000,2010))

# scale early reward band studies to $100 based on Nichols et al. 1991 relationship
reporting$rho <- reporting$rho*reporting$correction
# subset to remove early Martinson papers
reporting2 <- subset(reporting,reporting$survey!=1)

# plot raw data
plot(reporting$r_yr,reporting$rho,ylim=c(0,1),ylab="Reporting rate",xlab="")
legend(1970,1,legend=c("Inflation adjusted rates"),bty="n")
points(reporting2$r_yr,reporting2$rho,pch=19,col="red")
# note extreme outlier effect of Geis & Martinson papers

reporting$year <- reporting$r_yr-1950
reporting$rev.yr <- abs(reporting$r_yr-2011)
max.yr <- max(reporting$rev.yr)

# exploratory models
summary(m0<-lm(rho~1,data=reporting))
# all spp, all flyways
summary(m1<-lm(rho~year+I(year*year)+reward+I(reward*value)+tollfree+spp+flyway+I(year*tollfree),data=reporting))
# trophy spp only, E Canada only
summary(m2<-lm(rho~year+I(year*year)+reward+I(reward*value)+tollfree+duck+trophy+E_Can+I(year*tollfree),data=reporting))
AIC(m0,m1,m2) # MLE favors simpler model with no flyway or species effects (E_Can & trophy only).

# this covariate based on proportion 1-800 or web bands in recovery year
new.bands.recov <- c(0.975,0.974,0.969,0.960,0.955,0.943,0.919,0.889,0.843,0.782,0.712,0.616,0.527,0.384,0.03,rep(0,48))
new.bands.deploy <- c(0.995,0.995,0.985,0.987,0.989,0.983,0.981,0.967,0.917,0.912,0.857,0.832,0.756,0.583,0.125,rep(0,48))
year <- (1948:2010)
plot(rev(year),new.bands.recov,xlim=c(1995,2010), type="b", ylab="Proportion of new bands", xlab="Year")
points(rev(year),new.bands.deploy,col="red", type="b")
legend(2002,0.75,legend=c("Deployed","Recovered"),text.col = c("red","black"),bty="n")

old.bands.recov <- 1-new.bands.recov
old.bands.recov2 <- c(old.bands.recov[-1],1)
change1 <- old.bands.recov2-old.bands.recov
old.bands.deploy <- 1-new.bands.deploy
old.bands.deploy2 <- c(old.bands.deploy[-1],1)
change2 <- old.bands.deploy2-old.bands.deploy
plot(change2)
sum(change2)

reporting.data <- list(change1 = change1, change2 = change2, trophy = reporting$trophy,
                       max.yr = max.yr, rev.yr = reporting$rev.yr, reward = reporting$reward,
                       pub = reporting$n.pub, AdM = reporting$AdM, survey = reporting$survey,  
                       spp = reporting$n.spp, nspp = max(reporting$n.spp), 
                       rho.obs = reporting$rho, rho.tau = 1/(reporting$se_rho^2), 
                       survey2 = reporting$survey2, solicit = reporting$solicit, 
                       tollfree = reporting$tollfree, duck = reporting$duck, 
                       flyway = reporting$n.fly, nfly = max(reporting$n.fly),
                       E_Can = reporting$E_Can, n.obs = dim(reporting)[1])

sink("reporting.jags")  
cat("
model {
  # Priors
    rho.init ~ dunif(0.6,0.9) # somewhat informative prior for 2011 reporting rate
    l.rho.init <- log(rho.init/(1-rho.init))
    b.band ~ dunif(-4,4) # change in rho due to change in band type
    b.duck ~ dunif(-2,2) # contrast ducks vs. geese
    b.trophy ~ dunif(-2,2) # fixed effect for CANV, NOPI, REDH
    b.E_Can ~ dunif(-2,2) # Francophone effect
    b.solicit ~ dunif(-2,2) # solicitation effect
    b.survey ~ dunif(-2,2) # early PCS surveys (e.g. Atwood, Martinson)
    b.survey2 ~ dunif(-2,2) # differentiate Padding * Royle
    sd.yr ~ dunif(0,2)
    tau.yr <- pow(sd.yr,-2)

  # likelihood
    # assign initial values to year 1
      eps.yr[1] ~ dnorm(0,tau.yr)
      l.rho[1] <- l.rho.init + b.band*change2[1] + eps.yr[1]
    # trend is 1st-order Markovian state-space process with band type effect
    for (i in 2:max.yr){
      eps.yr[i] ~ dnorm(0,tau.yr)
      l.rho[i] <- l.rho[i-1] + b.band*change2[i-1] + eps.yr[i]      
    } # close i
  # obs model
    for (n in 1:n.obs){
      logit(rho[n]) <- l.rho[rev.yr[n]] + b.trophy*trophy[n] + b.survey2*survey2[n]
        + b.E_Can*E_Can[n] + b.survey*survey[n] + b.solicit*solicit[n] + b.duck*duck[n]
      rho.obs[n] ~ dnorm(rho[n],rho.tau[n])
      }
    }    # end jags model
    ",fill=TRUE)
sink()

parms <- c("rho","l.rho","b.band","b.trophy","b.E_Can","b.survey","b.survey2",
           "b.solicit","b.duck","sd.yr")
na = 1000; nc = 3; nt = 1; ni = 60000; nb = 10000
n.post <- nc*(ni-nb)/nt

# change name to rho3 for reporting3 data
rho.out <- jagsUI(reporting.data, inits=NULL, parms, "reporting.jags", 
                  n.adapt=na, n.chains=nc, n.thin=nt, n.iter=ni, n.burnin=nb,
                  parallel=TRUE, codaOnly = c("rho"))

print(rho.out)

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

save(rho.out, file="rho_out.Rdata")

means <- c(rho.out$mean$b.band,rho.out$mean$b.trophy,rho.out$mean$b.E_Can,
           rho.out$mean$b.survey,rho.out$mean$b.survey2,rho.out$mean$b.solicit,
           rho.out$mean$b.duck)
sds <- c(rho.out$sd$b.band,rho.out$sd$b.trophy,rho.out$sd$b.E_Can,
           rho.out$sd$b.survey,rho.out$sd$b.survey2,rho.out$sd$b.solicit,
           rho.out$sd$b.duck)

plot(means,ylim=c(-2.5,2.5),ylab="Parameter values",pch=19)
segments(1:7,means-1.96*sds,1:7,means+1.96*sds)
abline(h=0,col="black")

# calculate effect sizes for mean 0.5 or 0.7 on real scale
int.50 <- rep(0,7)
int.70 <- rep(0.85,7)
plogis(int.50)-plogis(int.50+means)
plogis(int.70)-plogis(int.70+means)

plot(rev(1948:2010),rho.out$sd$l.rho,pch=19,ylim=c(0,0.3),
     type='b',ylab="SD logit p",xlab="")

redh.rho <- E_Can.rho <- matrix(NA,nrow=n.post,ncol=max.yr)
for (i in 1:n.post){
  redh.rho[i,1:63] <- plogis(rho.out$sims.list$l.rho[i,1:63]+rho.out$sims.list$b.trophy[i])
  E_Can.rho[i,1:63] <- plogis(rho.out$sims.list$l.rho[i,1:63]+rho.out$sims.list$b.E_Can[i])
}
redh.mean <- redh.sd <- E_Can.mean <- E_Can.sd <- numeric()
redh.mean <- colMeans(redh.rho)
E_Can.mean <- colMeans(E_Can.rho)
for (i in 1:63){
  redh.sd[i] <- sd(redh.rho[,i])
  E_Can.sd[i] <- sd(E_Can.rho[,i])
}
# estimates and SD for trophy (REDH), Eastern Canada, and "mallards"
rev(redh.mean)
rev(redh.sd)
rev(E_Can.mean)
rev(E_Can.sd)
rev(plogis(rho.out$mean$l.rho))
rev(plogis(rho.out$mean$l.rho+rho.out$sd$l.rho)-plogis(rho.out$mean$l.rho-rho.out$sd$l.rho))/2

# plot annual reporting rates for reward band studies (typical spp)
plot(rev(1948:2010), plogis(rho.out$mean$l.rho), ylim=c(0,1), pch=19, type="b",
     ylab="Reporting probability", xlab="Year", cex.lab=1.2, cex.axis=1.2, bty="l")
segments(rev(1948:2010), plogis(rho.out$mean$l.rho-1.96*rho.out$sd$l.rho), 
         rev(1948:2010), plogis(rho.out$mean$l.rho+1.96*rho.out$sd$l.rho))
# include only years when trophy ducks were sampled
points(rev(1954:1960)-0.2, redh.mean[51:57], ylim=c(0,1), pch=21, type="b")
segments(rev(1954:1960)-0.2, redh.mean[51:57]-1.96*redh.sd[51:57], rev(1954:1960)-0.2, redh.mean[51:57]+1.96*redh.sd[51:57], lty=2)
# include only years when Eastern Canada was segregated
points(rev(2002:2010)+0.2, E_Can.mean[1:9], ylim=c(0,1), pch=22, type="b")
segments(rev(2002:2010)+0.2, E_Can.mean[1:9]-1.96*E_Can.sd[1:9], rev(2002:2010)+0.2, E_Can.mean[1:9]+1.96*E_Can.sd[1:9], lty=2)
legend(1970, 0.8, legend=c("Trophy","Mallard","E Canada"), bty="n", lty=c(2,1,2), 
       pch=c(21,19,22), cex=1.2)