#------------#
# JAGS model #
#------------#

# this code runs the WT_Climate model

cat("
model{
  # Prior distributions
  # potential Regression parameters
  alpha0 ~ dunif(-10,10)  # intercept detection prob on sigma (shape parameter)
  rugd ~ dunif(-10,10)
  asp ~ dunif(-10,10)
  slp ~ dunif(-10,10)
  
  for(j in 1:Nlev){           #3 levels fixed effect
    beta0[j] ~ dunif(-10,10)  # intercept initial density/Abundance
    betat0[j] ~ dunif(-10,10) # intercept density dynamic model
  }
  
  taubtDD <- pow(2,-2)
  btDD ~ dnorm(0,taubtDD)I(-2,2) # DD parameter
  btTREND  ~ dnorm(0,100)  # Excess trend in growth
  btPrect ~ dunif(-10,10)   # max precipitation effect 
  btROSt ~ dunif(-10,10)    # ROS days effect
  btTempt ~ dunif(-10,10)   # temperature effect
  btWiOnt ~ dunif(-10,10)    # Winter Onset effect
  btWiTemp ~ dunif(-10,10)    # Winter Temperature effect
  btCarct ~ dunif(-10,10)    # temporal carcass effect
  btHarvt ~ dunif(-10,10)    # Harvest effect
  #btSeaIcet ~ dunif(-10,10)    # Sea ice effect

  
  ## Alternative specification of precision via inv gamma distr
  PrOc ~ dgamma(alphaProc, betaProc) # approximates inv.gamma with vague priors, alpha and beta = 0.01
  PrEc ~ dgamma(alphaPrec, betaPrec) # alpha and beta = 0.01
  sdproctau <- 1/sqrt(PrOc)
  sdprectau <- 1/sqrt(PrEc)

  
  ## Definition of random transect cluster effect - to take into account that points are clustered and non-independent
  for (j in 1:NClust) # NClust = number of clusters (areas)
  {rCl1[j]~ dnorm(0,rtau) 
    rCl[j]~ dnorm(0,rtau2)}
  rtau ~ dgamma(alphaTau, betaTau)
  rtau2 ~ dgamma(alphaTau2, betaTau2)
  sdrtau <- 1/sqrt(rtau)
  sdrtau2 <- 1/sqrt(rtau2)
  
  # 'Likelihood'
  for (s in 1:nsites){
    # Linear model for detection function scale
    log(sigma[s]) <- alpha0 + rugd*vrm_d[s] + asp*aspect_d[s] + slp*slope_d[s]
    # Compute detection probability
    for(k in 1:nD){
      log(p[s,k]) <- -midpt[k]*midpt[k]/(2*sigma[s]*sigma[s]) # Half-normal detection function
      f[s,k] <- p[s,k]*pi[s,k]
      fc[s,k] <- f[s,k]/pcap[s]
      fct[s,k] <- fc[s,k]/sum(fc[s,1:nD])
      pi[s,k] <- (2*midpt[k]*delta )/(B*B)
    }
    pcap[s]<-sum(f[s,1:nD])  # Overall detection probability, i.e. sum over all bins
    
    # Process model
    # Abundance/density model for Yr1 as in Sillett et al 2012
    y[s,1] ~ dbin(pcap[s], (N[s,1]*Nrep[s,1]))  # measurement error
    N[s,1] ~ dpois( lambda[s,1] )  # poisson variation     # N is poisson with expected value lambda
    lambda[s,1] <-  D[s,1] * areadet[s]
    logD[s,1] ~ dnorm( mu[s,1], PrEc )  
    mu[s,1] <- beta0[ThreeLev[s]] + rCl1[Clust[s]] 
    D[s,1] <- exp(logD[s,1])
    
    # Population dynamics model for subsequent years
    for (t in 2:T){
      y[s,t] ~ dbin(pcap[s], (N[s,t]*Nrep[s,t]))
      N[s,t] ~ dpois(lambda[s,t])  ## poisson variation
      lambda[s,t] <- D[s,t] * areadet[s]
      logD[s,t] ~ dnorm( mu[s,t] , PrOc ) # precision Process variance/error. High value means high precision.
      # Autoregressive model: mu is the latent state to be estimated
      mu[s,t] <- betat0[ThreeLev[s]] + rCl[Clust[s]] + btDD * mu[s,t-1] +
      btHarvt*Harvt[t-1] +
      btTempt*Temp1Julyt[t-1] +
      btPrect*Prec1Julyt[t-1] +
      btWiOnt*WiOnt[t-1] +
      btROSt*ROS_days[t] +
      btWiTemp*WiTemp[t] +
      btCarct*Carcass[t-1] +
      btTREND * (t-1) #+
      # btSeaIcet*SeaIce[t]
      # model on growth rates because of delayed effect (btDD), i.e. effect estimates are on the growth rate
      D[s,t] <- exp(logD[s,t])
    }
  }
  
  # Distance sampling observation model for observed (binned) distance data
  for(i in 1:nobs){
    dclass[i] ~ dcat(fct[site[i],]+0.01)  
  }
  
  for(s in 1:nsites){
    for(t in 1:T){
      PredY[s,t] <- ((exp(mu[s,t]+ 0.5*(sdproctau*sdproctau))*areadet[s])*Nrep[s,t]) * pcap[s]     
    }    } # prediction
  
  # Derived quantities:
  for(t in 1:T){
    Ntot[t] <- sum(N[,t])
    Dest[t] <- Ntot[t] / sum(areadet) 
  }
}
", file="DynDensSvalbardPt.txt")