---
title: "Analysis Code - Burt et al."
author: "Sudden collapse of a mesopredator reveals its complementary role in mediating rocky reef regime shifts"
date: "Mar 01, 2018"
output: 
  html_document: 
    toc: yes
  pdf_document: 
    highlight: tango
    toc: yes
    toc_depth: 6
---

### Load packages
```{r load_libraries, message=FALSE}
library(stats)
library(lattice)
library(coda)
library(boot)
library(readxl)
library(rjags)
library(runjags)
library(matrixStats)
library(ggplot2)
library(cowplot)
library(reshape2)
library(magrittr)
library(plyr)
library(dplyr)
library(tidyr)
library(parallel)
library(loo)
library(Rcpp)
library(MASS)
```

### Load the data
** see README file for details
```{r load_data}
# set appropriate file directory
DatQ <- read.csv("Burt_CC-DiveSurveys2013-2016_QuadCounts.csv") #Quadrat-level data, counts of urchins and kelp
DatS <- read.csv("Burt_CC-DiveSurveys2013-2016_Site_Sum_stats.csv") #Site-level data for otters, pycno, urchins, kelp
DatU = read.csv("Burt_CC-DiveSurveys2013-2016_UrchSize.csv")  #Every measured urchin in all quadrats from 2013-2016. 
attach(DatS)
attach(DatQ)
attach(DatU)
```

**Data legend for DatQ 
Y = survey year  
S = site  
Q = quadrat  
sm = number of small red urchins (<=3cm) in a given 1x1m quadrat  
med = number of medium red urchins (>=4cm & <=7cm) in a given 1x1m quadrat  
large = number of large red urchins (>=8cm) in a given 1x1m quadrat  
TotalCount = total count of all red urchins in a 1x1m quadrat  
KelpCount = total count of adult kelp (stipe length >15cm) for quadrats deeper than 5m depth  

**Data legend for DatS
Year = survey year  
SiteNum = site, numbered 1-11  
Otter = otter occupancy of a site (otter occupied = 1, otters absent = 0)  
OtterCat = otter status at a site is "continuous", "absent", or "newly arrived"  
Pycno = biomass (kg) of Pycnopodia per 10m2  
Predscen = "predator scenario" = "both" (otters present & Pycno present ie. >0.09kg/10m2), "PYConly" (otters absent, Pycno present), "OTTonly" (otters present & Pycno absent - ie. < 0.09kg/10m2), "none" (otters absent & Pycno absent)  
AvgOfSm = mean density of small red urchins at a site (inidivuals/m2)  
AvgOfmed = mean density of medium red urchins at a site  (inidivuals/m2)  
AvgOflg = mean density of large red urchins at a site  (inidivuals/m2)   
AvgDens = mean density of all red urchins at a site (individuals/m2)  
AvgDens5 = mean density of all red urchins at a site (individuals/m2) for quads deeper than 5m  (Fig. 6)  
Dens5se = standard error around mean density of all red urchins at a site  
BiomOfsm = mean biomass of small red urchins at a site (kg urchins/m2)  
BiomOfmed = mean biomass of medium red urchins at a site (kg urchins/m2)  
BiomOflg = mean biomass of large red urchins at a site (kg urchins/m2)  
AvgBiom5 = mean biomass of all red urchins at a site (kg urchins/m2) for quad deeper than 5m (Fig. 6)  
AvgKelp = mean density of adult kelp spp. (stipes/m2) for quad deeper than 5m (used in Figs. 1 & 5)  
KelpSE = standard error around mean density of adult kelp spp. (stipes/m2)  
Perennials = mean perennial kelp species at a site (stipes/m2) 

**Data legend for DatU
UrchSize = individual sizes (measured to nearest 1cm) of all the individual red urchins 
SizeCat = assigned red urchin size category (“sm” = 1-3cm, “med” = 4-7cm, “lg” = 8cm and up)



##Obtain urchin growth transition probabilities 
### Run simulation to estimate growth transition probabilities 
Simulate annual growth increment for a large sample of urchins whose starting size is given by the size distribution from urchin plots.
Simulated growth based on the stochastic growth model presented in: Zane Zhang, Alan Campbell, Dominique Bureau 2008. Growth and Natural Mortality Rates of Red Sea Urchin (Strongylocentrotus franciscanus) in British Columbia. Journal of Shellfish Research, 27(5):1291-1299.
``` {r growth_sim, eval=FALSE}
# From Zhang et al. (2008) - use Price Islands Islands data for growth curve, Tanaka function
# Paramaters for Tanaka function:

a = 0.00539
d = 118.2
f = 0.00342
sig = 1.21

reps = 5000
# Define size classes:
Szclases = c('sm', 'med', 'lg')
Grow = c(10, 40, 75)
G = matrix(nrow=reps,ncol=2)
for (r in 1:reps){
  Dt = DatU$UrchSize*10 + runif(length(DatU$UrchSize), -4.999, 4.999)
  E = exp(sqrt(f)*(Dt-d))
  Dt1 = 1/sqrt(f)*log(2*f*(E/(4*f) - a/E +1) + 2*sqrt((f^2)*(E/(4*f) - a/E +1)^2 + f*a )) + d + rnorm(length(Dt),0,sig) 
  Dt1 = pmax(Dt,Dt1)
  for(i in 1:2){
    ii = which(Dt>=Grow[i]-5 & Dt<=(Grow[i+1]-5))
    tmp1 = Dt[ii]
    tmp2 = Dt1[ii]
    j = which(tmp2>=(Grow[i+1]-5))
    G[r,i] = length(j)/length(ii)
  }
}
G_hat = colMeans(G);
print(paste0('G prob values for sm and med: ', format(G_hat[1],digits=3),
             ', ', format(G_hat[2],digits=3)))
#
#create dataframe for plotting
# urchin growth increments to produce a scatter plot of random growth generated from the tanaka model
datagraph = data.frame(Dt = Dt[ii], Dt1 = Dt1[ii]-Dt[ii] ) 
# here is the output of the simulated growth transition probabilities to make a boxplot
databoxplot = data.frame(Small = G[,1], Medium = G[,2], Run="NA")

# Make plots
#  Scatter plot of random growth 
ii =  sample(seq(1,length(DatU$UrchSize)),500)
smth = smooth.spline(Dt,Dt1-Dt) 
plot(Dt[ii],Dt1[ii]-Dt[ii],xlim=c(0,120),xlab="Uchin size (mm)",ylab="Growth increment (mm)",
     main = "Urchin Growth, Tanaka model")
lines(smth, col='blue')
# Boxplot of growth transition probabilities
boxplot(G, names = c('Small','Medium'), xlab = "Size class", ylab = "Proportion Growing") 

```


##Bayesian model description
### Prepare data for running the model
```{r prep_data}
# These are the outputs from the growth model, stage-specific transition probabilities
Growth_Sm_Md = 0.52
Growth_Md_Lg = 0.09
# Log normal params for Recruitment prior: gives mean ~5
Rpri1 = 1.4
Rpri2 = 2.2
#
SiteN = SiteNum; 
YearN = Year - min(Year) + 1
Yr = Y -  min(Year) + 1
St = S 
Nsites =  max(SiteN)
Nyears =  max(YearN)
Nobs = length(Yr) 
Ott = matrix(0,nrow = Nsites, ncol = Nyears)
Pyc = matrix(0,nrow = Nsites, ncol = Nyears)
counts = cbind(sm, med, lg)
N0 = cbind(AvgOfsm[1:Nsites], AvgOfmed[1:Nsites], AvgOflg[1:Nsites])

# Deternine Otter and Pycno status (relative abund, 0-1 scale) at each site/year
for (s in 1:Nsites){
   for (y in 1:Nyears){
     Ott[s,y] = mean(Otter[which(SiteN==s & YearN==y)])
     Pyc[s,y] = min(1,mean(Pycno[which(SiteN==s & YearN==y)]))
   }
} 
```

### Create data and initials objects for JAGS
```{r jags.data}
jags.data <- list(NS = Nsites, obs = counts, NY = Nyears, NEQ = 20, 
                  Y = Yr, S = St, NObs = Nobs, g1 = Growth_Sm_Md,  
                  g2 = Growth_Md_Lg, O = Ott, P = Pyc, N0 = N0, 
                  Rpri1=Rpri1, Rpri2=Rpri2)  
inits <- function() list(Pobs = runif(1, .1, .5), dbar = runif(1, .5, 2),
                         beta1 = runif(1, .1, .5), beta2 = runif(1, .05, .3),
                         sigS = runif(1, .5, 2), sigY = runif(1, .5, 2),
                         Dispers = runif(1, .2, 1), 
                         alphaP = runif(3, .1, 1),
                         alphaO = runif(2, .1, 1)) 

params <- c("dbar","Pobs","sigS","sigY","Dispers","R",
            "beta1","beta2","d","alphaP","alphaO","Rt",
            "Ns","Nm","Nl","NsP","NmP","NlP",
            "gamma_s", "gamma_m", "gamma_l") 
nsamples <- 1000
nt <- 5
nb <- 2000

# For parallel (comment out for serial)
library(parallel)
cores = detectCores()
ncore = min(20,cores-1)
cl <- makeCluster(ncore)
nc <- ncore
```


### Create the jags model
```{r jags_file, eval=FALSE}

jags_model <- "model {
  # Model: Otter and Pycno size-specific pred, R variable (fitted params)
	# Set Priors: 
    Dispers ~ dgamma(1.1,0.1)   # over-dispersion parameter (for Neg Binom counts)
    # Mean per-capita death rate for small urchins - fitted 
    dbar ~ dunif(.0001,10) # 
    Pobs ~ dbeta(1.5,1.5) # vague prior, prob of counting small szclass urchins
    # Mortality rate params for med and large: computed as a fraction
    # of the rate for small (so that mortality at each site correlated accross size classes)
    beta1 ~ dbeta(1, 2)
    beta2 ~ dbeta(1, 2)
    # Predator-specific mortality rates for pycno and otters (flat priors)
    # ***NOTE*** IN THIS CODE & ANALYSIS size class 1 = rate for MED urchins, 2 = rate for LG urchins, and 3 = rate for SM urchins
    # ********** IN THE PAPER, for clarity, size class 1 = SM urchins, size class 2 = MED urchins, and size class 3 = LG urchins
	# NOTE: use vague normal priors to allow positive OR negative estimates
    alphaP[1] ~ dunif(0.0001,5)
    alphaP[2] ~ dunif(0.0001,5)
    alphaP[3] ~ dunif(0.0001,5)
    alphaO[1] ~ dunif(0.0001,5)
    alphaO[2] ~ dunif(0.0001,5)
    # Variance params: variance among sites/years in mortality rates
    sigS ~ dunif(0.001,10)
    tauS <- pow(sigS,-2)
    sigY ~ dunif(0.001,10)
    tauY <- pow(sigY,-2)
	# R is allowed to vary over years, calculate R[t]...
	for (t in 1:NY){
		# Annual Apparent Annual Recruitment to smallest size class 
		Rt[t] ~ dlnorm(Rpri1,Rpri2) # slightly informed prior centered around 5
	}
	# Mean Apparent Annual Recruitment to smallest size class 
	R <- mean(Rt)
    # Compute hidden state dynamics at each site  
    for (s in 1:NS){
      # Site-specific mortality (leads to variation in equilibrial abundance)
      d[s] ~ dnorm((dbar),tauS) T(0,)
      b1[s] <- d[s]*beta1
      b2[s] <- d[s]*beta2
      # Initialize site abundances at year 1 
      NsEq[s,1] <- N0[s,1]*(1/Pobs)
      NmEq[s,1] <- N0[s,2]
      NlEq[s,1] <- N0[s,3]
      # Cycle through initialization years to reach equilibrium abundance for each size class
      for (t in 2:(NEQ)){
        NsEqP[s,t] <- NsEq[s,t-1] + R - g1*NsEq[s,t-1] 
        NmEqP[s,t] <- NmEq[s,t-1] + g1*NsEq[s,t-1] - g2*NmEq[s,t-1] 
        NlEqP[s,t] <- NlEq[s,t-1] + g2*NmEq[s,t-1]        
        ### A) INTERMEDIATE STEP TO GET GAMMAS
        gamtmp_s[s,t-1] <- d[s] + alphaP[3]*P[s,1]
        gamtmp_m[s,t-1] <- b1[s] + alphaP[1]*P[s,1] + alphaO[1]*O[s,1]
        gamtmp_l[s,t-1] <- b2[s] + alphaP[2]*P[s,1] + alphaO[2]*O[s,1]        
        ### B) COMPUTE DETERMINISTIC SURVIAL
        NsEq[s,t] <- NsEqP[s,t]*exp(-1*(gamtmp_s[s,t-1]))
        NmEq[s,t] <- NmEqP[s,t]*exp(-1*(gamtmp_m[s,t-1]))
        NlEq[s,t] <- NlEqP[s,t]*exp(-1*(gamtmp_l[s,t-1]))
      }
	  # Year 1 abundances (equilibrium with stochasticity)	 
	  # First, the pre-survival abundance  based on growth and recruitment
	  NsP[s,1] <- NsEq[s,NEQ] + Rt[1] - g1*NsEq[s,NEQ]  	
      NmP[s,1] <- NmEq[s,NEQ] + g1*NsEq[s,NEQ] - g2*NmEq[s,NEQ] 
      NlP[s,1] <- NlEq[s,NEQ] + g2*NmEq[s,NEQ]    	       
      # GAMMAS: cumulative log hazards, which determine mean expected survival
      gamma_s[s,1] <- d[s] + alphaP[3]*P[s,1]
      gamma_m[s,1] <- b1[s] + alphaP[1]*P[s,1] + alphaO[1]*O[s,1]
      gamma_l[s,1] <- b2[s] + alphaP[2]*P[s,1] + alphaO[2]*O[s,1]
	  # Stochasticity in log hazards, Year 1
	  eps0[s,1] ~ dnorm(0,tauY)
      eps1[s,1] ~ dnorm(0,tauY)
      eps2[s,1] ~ dnorm(0,tauY)	  
      # Calculate N at year 1, factoring in survival with stochasticity
      Ns[s,1] <- NsP[s,1]*exp(-1*(gamma_s[s,1])*exp(eps0[s,1]))
      Nm[s,1] <- NmP[s,1]*exp(-1*(gamma_m[s,1])*exp(eps1[s,1]))
      Nl[s,1] <- NlP[s,1]*exp(-1*(gamma_l[s,1])*exp(eps2[s,1]))
      # Now cycle through remaining real Years (allow for year to year stochasticity)
      for (t in 2:NY){
		# Stochasticity in log hazards, Year t
        eps0[s,t] ~ dnorm(0,tauY)
		eps1[s,t] ~ dnorm(0,tauY)
		eps2[s,t] ~ dnorm(0,tauY)
		# pre-survival abundance Year t based on growth and recruitment
        NsP[s,t] <- Ns[s,t-1] + Rt[t] - g1*Ns[s,t-1] 
        NmP[s,t] <- Nm[s,t-1] + g1*Ns[s,t-1] - g2*Nm[s,t-1] 
        NlP[s,t] <- Nl[s,t-1] + g2*Nm[s,t-1]        
        #GAMMAS: cumulative log hazards, which determine mean expected survival
        gamma_s[s,t] <- d[s] + alphaP[3]*P[s,t]
        gamma_m[s,t] <- b1[s] + alphaP[1]*P[s,t] + alphaO[1]*O[s,t]
        gamma_l[s,t] <- b2[s] + alphaP[2]*P[s,t] + alphaO[2]*O[s,t]        
        #Calculate N at year t, factoring in survival with stochasticity
        Ns[s,t] <- NsP[s,t]*exp(-1*(gamma_s[s,t])*exp(eps0[s,t]))
        Nm[s,t] <- NmP[s,t]*exp(-1*(gamma_m[s,t])*exp(eps1[s,t]))
        Nl[s,t] <- NlP[s,t]*exp(-1*(gamma_l[s,t])*exp(eps2[s,t]))
      }
    }
    # Now determine stochastic nodes = measurements at each Quadrat, by size class
    for (i in 1:NObs){
	  # Mean expected abundance 
      mu[i,1] <- Ns[S[i],Y[i]]*Pobs
      mu[i,2] <- Nm[S[i],Y[i]]
      mu[i,3] <- Nl[S[i],Y[i]]
      # Scale Parameter for negative binomial
      p[i,1] <- Dispers/(Dispers+mu[i,1]) 
      p[i,2] <- Dispers/(Dispers+mu[i,2]) 
      p[i,3] <- Dispers/(Dispers+mu[i,3]) 
      # Observed Nodes: counts of small, medium and large urchins (with over dispersion):
      obs[i,1] ~ dnegbin(p[i,1],Dispers)
      obs[i,2] ~ dnegbin(p[i,2],Dispers)
      obs[i,3] ~ dnegbin(p[i,3],Dispers)
    }
  }"
```

### Call JAGS from R and run jags

```{r run_jags, eval=FALSE}
out <- run.jags(data = jags.data,
                inits = inits,
                monitor = params,
                model = jags_model,
                n.chains = nc,
                thin = nt,
                sample = nsamples,
                burnin = nb,
                method="rjparallel", cl=cl)

post = rbind(out$mcmc[[1]], out$mcmc[[2]])
for (i in 3:nc){
  post = rbind(post, out$mcmc[[i]])
}

sumstats = summary(out)
vn = row.names(sumstats)
# create matrix of all posterior chains (coda)
```


### Run posterior predictive checks
```{r post-pred-checks, eval=FALSE}
mcmc <- out$mcmc
require(dplyr)
data <- DatQ%>%
  dplyr::select(Y,S,Q,sm,med,lg)%>%
  dplyr::rename(s= sm,m=med,l= lg,time=Y,site=S)%>%
  melt(id.vars= c("time","site","Q"),value.name= "obs",variable.name="class")%>%
  mutate(time= time-min(time)+1)%>%
  mutate(time= as.character(time),
         site= as.character(site),
         class= as.character(class))

### function to generate pp checks
post_pred <- function(mcmc,data){

  ### need to get appropriate posteriors
  char_list_N<- c("Ns\\[","Nm\\[","Nl\\[")
  char_list_P <- c("Pobs")
  char_list_D <- c("Dispers")
  index_N <- grep(paste(char_list_N,collapse= "|"),dimnames(mcmc[[1]])[[2]])
  index_P <- grep(paste(char_list_P,collapse= "|"),dimnames(mcmc[[1]])[[2]])
  index_D <- grep(paste(char_list_D,collapse= "|"),dimnames(mcmc[[1]])[[2]])

  ### add iteration numbers
  post <- ldply(mcmc)%>%
    add_rownames(var= "iter")%>%
    mutate(iter= as.numeric(as.character(iter)))

  ### abundance estimate posteriors
  post_N <- post[,c(1,index_N+1)]%>%
    melt(id.vars= "iter",value.name= "N")%>%
    separate(variable,into= c("class","site","time"),sep= c("\\[|\\,|\\]"),
             extra = "drop", fill = "right")%>%
    mutate(class = gsub("N|P","",class))

  ### posterior expectation
  post_P <- post[,c(1,index_P+1)]%>%
    melt(id.vars= "iter",value.name= "P")%>%
    mutate(class= "s")%>%
    dplyr::select(-one_of("variable"))

  ### negative binomial count posterior
  post_D <- post[,c(1,index_D+1)]%>%
    melt(id.vars= "iter",value.name= "D")%>%
    dplyr::select(-one_of("variable"))

  ### join posteriors
  join_post <- join_all(list(post_N,post_D,post_P))
  join_post <- join_post%>%
    mutate(P= ifelse(is.na(P),1,P))%>%
    mutate(mu= N*P)

  join_post2 <- full_join(join_post,data)

  ### generate random draws for comparisons with posterior observations
  join_post2<- join_post2%>%
    filter(iter%in%base::sample(unique(iter),500))%>%
    mutate(rand = rnegbin(nrow(.),mu,D),
           rand2 = rnegbin(nrow(.),mu,D))%>%
    mutate(rand_pr = (rand-mu)/sqrt(mu+mu^2/D),
           rand_pr2 = (rand2-mu)/sqrt(mu+mu^2/D),
           obs_pr = (obs-mu)/sqrt(mu+mu^2/D))

  ### generate test statistics
  ppc <- join_post2%>%
    group_by(iter)%>%
    mutate(b1= rank(rand_pr),
           c1= rank(rand_pr2),
           a1= rank(obs_pr))%>%
    summarize(null.ks= ks.test(rand_pr,rand_pr2)$statistic[[1]],
              obs.ks= ks.test(obs_pr,rand_pr2)$statistic[[1]],
              null.chisq = sum((rand-mu)^2/mu),
              obs.chisq = sum((obs-mu)^2/mu),
              null.prcor = sum((a1-mean(a1))*(c1-mean(c1)))/sqrt(sum((a1-mean(a1))^2)*sum((c1-mean(c1))^2)),
              obs.prcor = sum((b1-mean(b1))*(c1-mean(c1)))/sqrt(sum((b1-mean(b1))^2)*sum((c1-mean(c1))^2)))

  ### reshape for plotting
  ppc2 <- reshape(ppc, varying = 2:7, sep = ".", direction = 'long',timevar= "metric")

  ppc2 <- ppc2%>%
    mutate(metric= factor(metric, labels= c("a) Goodness-of-fit","b) Kolmogorov Smirnov","c) Pearson Rank Correlation")))

  ### generate p-values
  ppc_stat = ppc2%>%
    group_by(metric)%>%
    summarize(p= as.character(paste("P =",sprintf("%.2f",round(mean(obs<null,na.rm=T),digits= 2)))))

  ### plot post-pred checks
  plot <- list()
  for(i in 1:length(unique(ppc2$metric))){
    plot_dat <-  ppc2%>%filter(metric==levels(ppc2$metric)[i])
    ranges <- range(c(plot_dat$obs,plot_dat$null),na.rm=T)
    plot[[i]] <- ggplot(aes(obs,null),data=plot_dat)+
      geom_point(alpha= 0.25)+
      facet_wrap(~metric)+
      geom_abline(slope=1,yintercept= 0)+
      geom_text(aes(label=p),x= -Inf,y= Inf,data= ppc_stat%>%filter(metric==levels(ppc2$metric)[i]),hjust=-1, vjust=2)+
      theme_bw()+
      coord_equal(ratio= 1)+
      xlim(ranges) +
      ylim(ranges) +
      ylab(ifelse(i==1,"expected",""))+
      xlab("observed")+
      theme(strip.background=element_blank(),
            panel.grid=element_blank(),
            strip.text=element_text(hjust=0,size= 9),
            axis.text= element_text(size= 9),
            axis.title= element_text(size= 9))
  }

  plots <- plot_grid(plot[[1]],plot[[2]],plot[[3]],ncol= 3,align= "hv")

  return(plots)
}

### run posterior predictive checks
mcmc <- out$mcmc
post_pred(mcmc,data)

```


### Diagnostic plots
``` {r diagnostic_plots, eval=FALSE}
#show trace plots and posterior distributions 
plot(out, c('trace','histogram'),
     vars=c("dbar", "alphaP","alphaO", "Pobs","Dispers","R","sigS","sigY","beta1","beta2","R"))

#check paired correlations of urchin mortality parameters to check there isn't identifiability issues (specifically betas vs alphas)
post2 <- as.data.frame(post)   #turn posterior matrix into a dataframe
rand.rows = sample(1:nrow(post2),500)
plot(post2[rand.rows,c("beta2","beta1","dbar","alphaP[1]","alphaP[2]","alphaP[3]","alphaO[1]","alphaO[2]")])
cor(post2[rand.rows,c("beta2","beta1","dbar","alphaP[1]","alphaP[2]","alphaP[3]","alphaO[1]","alphaO[2]")])

```


## Working with outputs
### Bias correction on posterior estimates for alphas and d_bar
We did a bias correction for our parameter estimates within gamma because the mean of the lognormal increases with the variance. In order to have unbiased parameter estimates, we need to account for the fact that higher variance in epsilon inflates gamma. To do the bias correction, we simply multiply by exp((sigma^2)/2) for each iteration in the posterior.  
```{r bias_correct, eval=FALSE}

# loop through posterior iterations (i) and create new bias-corrected parameters
for (i in 1:nrow(post2)) {
    post2$dbar_c[i]    <- post2$dbar[i]        * exp((post2$sigY[i]^2)/2)
    post2$alphaP1_c[i] <- post2$`alphaP[1]`[i] * exp((post2$sigY[i]^2)/2)
    post2$alphaP2_c[i] <- post2$`alphaP[2]`[i] * exp((post2$sigY[i]^2)/2)
    post2$alphaP3_c[i] <- post2$`alphaP[3]`[i] * exp((post2$sigY[i]^2)/2)
    post2$alphaO1_c[i] <- post2$`alphaO[1]`[i] * exp((post2$sigY[i]^2)/2)
    post2$alphaO2_c[i] <- post2$`alphaO[2]`[i] * exp((post2$sigY[i]^2)/2)
}

post3 <- as.matrix(post2)  # the summary commands below prefer to work on a matrix. 

# create a new summary dataframe that now has the bias-corrected values in it 
# the NEW CORRECTED parameters are at the bottom of this file

sumstatsCORR=data.frame(lower95=colQuantiles(post3, probs = c(0.025)),  
                        upper95=colQuantiles(post3, probs = c(0.975)),
                        mean=colMeans(post3),             
                        sd=colSds(post3))                           

# HERE are the new corrected values (get rid of all the previous rows)
sumstatsCORRshort <- sumstatsCORR[-c(1:424),]
# use these to make Fig. 2g,h,i 

```


### Calculating Total Annual Urchin Mortality in 4 Predation Scenarios
```{r totalmortality, eval=FALSE}
#get 5000 random samples from the posterior distributions of all params
simreps = 1000     
reps = dim(post2) 
reps = reps[1]  
rs = sample(reps,simreps) 
postsamp = post2[rs,]  

# select just the parameters needed to calculate total annual mortality
sampclean <- postsamp[,c("dbar_c", "beta1", "beta2", "alphaP1_c", "alphaP2_c", "alphaP3_c", "alphaO1_c", "alphaO2_c")]
sampdf <- as.data.frame(sampclean)  #convert to a dataframe

#calculate the beta terms (baseline mort for med, lg)
sampdf$B1<- sampdf$beta1*sampdf$dbar_c
sampdf$B2<- sampdf$beta2*sampdf$dbar_c

#use pop equation + relevant params for all combos of pred scenario & urchin size class
sampdf$Lg_none <- 1 - exp(-(sampdf$B2))
sampdf$Lg_PYConly <- 1 - exp(-(sampdf$B2 + sampdf$alphaP2_c))
sampdf$Lg_OTTonly <- 1 - exp(-(sampdf$B2 + sampdf$alphaO2_c))
sampdf$Lg_both <- 1 - exp(-(sampdf$B2 + sampdf$alphaO2_c + sampdf$alphaP2_c))

sampdf$Med_none <- 1 - exp(-(sampdf$B1))
sampdf$Med_PYConly <- 1 - exp(-(sampdf$B1 + sampdf$alphaP1_c))
sampdf$Med_OTTonly <- 1 - exp(-(sampdf$B1 + sampdf$alphaO1_c))
sampdf$Med_both <- 1 - exp(-(sampdf$B1 + sampdf$alphaO1_c + sampdf$alphaP1_c))

sampdf$Sm_none <- 1 - exp(-(sampdf$dbar_c))
sampdf$Sm_PYConly <- 1 - exp(-(sampdf$dbar_c + sampdf$alphaP3_c))
sampdf$Sm_OTTonly <- 1 - exp(-(sampdf$dbar_c + 0))
sampdf$Sm_both <- 1 - exp(-(sampdf$dbar_c + 0 + sampdf$alphaP3_c))

sampdfclean<-sampdf[,11:22] #clean df so just columns of interest

sampdflong <- melt(sampdfclean) #put data in long format
#create new columns for 'size' and 'pred scenario'
sampdflong$size <- ifelse(sampdflong$variable %in% c("Lg_none","Lg_PYConly","Lg_OTTonly","Lg_both"),"Large",
                   ifelse(sampdflong$variable %in% c("Med_none","Med_PYConly","Med_OTTonly","Med_both"),"Medium",
                   ifelse(sampdflong$variable %in% c("Sm_none","Sm_PYConly","Sm_OTTonly","Sm_both"),"Small",NA)))
                            
sampdflong$predscen <- ifelse(sampdflong$variable %in% c("Lg_none","Med_none","Sm_none"),"none",
                       ifelse(sampdflong$variable %in% c("Lg_PYConly","Med_PYConly","Sm_PYConly"),"PYConly",
                       ifelse(sampdflong$variable %in% c("Lg_OTTonly","Med_OTTonly","Sm_OTTonly"),"OTTonly",
                       ifelse(sampdflong$variable %in% c("Lg_both","Med_both","Sm_both"),"both",NA))))

#take the mean of the of the random samples, & upper/lower quantiles
summarydf <- sampdflong %>% group_by(size, predscen) %>% summarize(mean = mean(value),
                                                                   upper95 = quantile(value, .95),
                                                                   lower95 = quantile(value, .05))

summarydf$mean <- summarydf$mean*100
summarydf$upper95 <- summarydf$upper95*100
summarydf$lower95 <- summarydf$lower95*100

names(summarydf)[names(summarydf)=="size"]  <- "Urchin Size Class"

#now have all the values to create figure 4a.

```


### Deterministic simulation of urchin abundance (in Supplementary)
```{r simulation, eval=FALSE}

simreps = 1000
reps = dim(post)  
reps = reps[1]    
rs = sample(reps,simreps)  
postsamp = post[rs,]   

PobsMn = sumstats[which(vn=="Pobs"),4]  
#this is the estimated MEAN value for Pobs (probability of being observed). 

NYEQ = 4   # 4 scenerios: No preds, Both preds, Otters, Pycno
P = c(0, 1, 0, 1); O = c(0, 1, 1, 0);
# "Critical density" for large/medium urchins (above which mortality would increase)
Kdens = 50

g1 = Growth_Sm_Md    #proportion of 'smalls' that move to 'meds' 
g2 = Growth_Md_Lg    #proportion of 'meds' that move to 'large'  
R = postsamp[,which(vn=="R")]
alphaP = matrix(nrow=simreps,ncol=3); alphaO = matrix(nrow=simreps,ncol=2);   
dbar = postsamp[,which(vn=="dbar")]
beta1 = postsamp[,which(vn=="beta1")]   # propotion param for death in med urch
beta2 = postsamp[,which(vn=="beta2")]  # propotion param for death in lg urch
Pobs = postsamp[,which(vn=="Pobs")]   # proportion param: proportion of smalls those we see
b1 = dbar*beta1      
b2 = dbar*beta2      

alphaP[,1] = postsamp[,which(vn=="alphaP[1]")]    
alphaP[,2] = postsamp[,which(vn=="alphaP[2]")]     
alphaP[,3] = postsamp[,which(vn=="alphaP[3]")]     
alphaO[,1] = postsamp[,which(vn=="alphaO[1]")]     
alphaO[,2] = postsamp[,which(vn=="alphaO[2]")]     

sigY =  postsamp[,which(vn=="sigY")]  #posterior distribution for annual stochasticity.
sigYmn = mean(sigY)                      

#make a place to store outputs for each run (n=5000 simreps), 4 columns = 4 scenarios
NeqS = matrix(nrow=simreps,ncol=4)
NeqM = matrix(nrow=simreps,ncol=4)
NeqL = matrix(nrow=simreps,ncol=4)

# make the starting value the mean number of urchins in 2013 when otters == 0 
#(account for the fact we only see a proportion of smalls, ie. divide by PobsMn). 
Ns = numeric(length=NYEQ); Ns[1] = mean(DatS$AvgOfsm[DatS$Year==2013 & DatS$Otter==0])/PobsMn 
Nm = numeric(length=NYEQ); Nm[1] = mean(DatS$AvgOfmed[DatS$Year==2013 & DatS$Otter==0])
Nl = numeric(length=NYEQ); Nl[1] = mean(DatS$AvgOflg[DatS$Year==2013 & DatS$Otter==0])

### loop through posterior iterations (r)
for (r in 1:simreps){
  ### loop through 4 predator scenarios (s)
  for (s in 1:4){
    ### loop through each "t"ime (ie. year) from year 2 to year NYEQ=4
    for (t in 2:NYEQ){
      # first start with what generates urchins in each size category
      NsP = Ns[t-1] + R[r] - g1*Ns[t-1] 
      NmP = Nm[t-1] + g1*Ns[t-1] - g2*Nm[t-1] 
      NlP = Nl[t-1] + g2*Nm[t-1]
      # calculate home many urchins will be removed based on predation & baseline death
      gamma_sm = dbar[r] + alphaP[r,3]*P[s]
      gamma_md = alphaP[r,1]*P[s] + alphaO[r,1]*O[s] + b1[r]
      gamma_lg = alphaP[r,2]*P[s] + alphaO[r,2]*O[s] + b2[r]
      
      # NOTE- multiply gamma by exp(sigYmn^2/2) to correct for bias due to lognormal error.
      Ns[t] = NsP*(exp(-1*(gamma_sm*exp((sigY[r]^2)/2))))
      Nm[t] = NmP*(exp(-1*(gamma_md*exp((sigY[r]^2)/2)))) 
      Nl[t] = NlP*(exp(-1*(gamma_lg*exp((sigY[r]^2)/2))))  
    }
    NeqS[r,s] = Ns[NYEQ]*Pobs[r]
    NeqM[r,s] = Nm[NYEQ]
    NeqL[r,s] = Nl[NYEQ]
  }
}
# 

#create dataframes of the size-specific outputs
d4s=data.frame(names=c("No Pred","Both Pred","Otters only","Pycno only"),
              groupN=c("Small"),
              mean=colMeans(NeqS),  
              lower=colQuantiles(NeqS, probs = c(0.05)),   
              upper=colQuantiles(NeqS, probs = c(0.95)))

d4m=data.frame(names=c("No Pred","Both Pred","Otters only","Pycno only"),
               groupN=c("Medium"),
               mean=colMeans(NeqM),
               lower=colQuantiles(NeqM, probs = c(0.05)), 
               upper=colQuantiles(NeqM, probs = c(0.95)))

d4l=data.frame(names=c("No Pred","Both Pred","Otters only","Pycno only"),
               groupN=c("Large"),
               mean=colMeans(NeqL),
               lower=colQuantiles(NeqL, probs = c(0.05)), 
               upper=colQuantiles(NeqL, probs = c(0.95)))

#create single dataframe
d4 = rbind(d4s,d4m,d4l)  # stick the 3 size groups all together

# graphical specifics
filcol <- c(Small="black", Medium="grey50", Large="white")
pd <- position_dodge(0.4) # move them .05 to the left and right
names(d4)[names(d4)=="groupN"]  <- "Urchin Size Class"

#make it so the order is Large then Medium then Small 
d4$`Urchin Size Class` <- factor(d4$`Urchin Size Class`, levels = rev(levels(d4$`Urchin Size Class`)))  #this reverses order

# plot simulated data
#error bars here at the 0.05% lower quantile and 0.95% upper quantiles.
ggplot(d4, aes(x=names, y=mean, group=`Urchin Size Class`)) + 
  geom_hline(yintercept=0, colour="lightgrey", linetype="dashed")+
  geom_errorbar(aes(ymin=lower, ymax=upper), colour="grey55", width=.1, position=pd) +
  geom_point(aes(size=`Urchin Size Class`, fill=`Urchin Size Class`), shape=21, position=pd) +
  xlab("") + ylab("Equilibrium urchin abundance \n (urchins/m2)") +
  scale_fill_manual(values=filcol) +
  scale_size_manual(values=c(Small=1.5,Medium=4,Large=7))+
  scale_x_discrete(limits=c("No Pred","Pycno only", "Otters only", "Both Pred"), 
                   labels=c("No Pred"="No preds", "Pycno only", "Otters only","Both Pred"="Both preds")) +
  #scale_y_continuous(limits = c(0, 14), expand = c(0, 1), breaks= c(0, 3, 6, 9, 12)) +
  coord_cartesian(ylim=c(-1, 13), xlim = c(.5, 4.5), expand = c(0, 1)) +
  theme_cowplot() +
  theme(legend.justification=c(1,0),
        legend.position=c(.95,.6))
```
