model;
{ #  2-d random-walk-with-drift hierarchical model
#    Default model from Ian D Jonsen, Ransom A Myers, & Michael C James
#    "Robust hierarchical state-space models reveal diel variation in travel 
#       rates of migrating leatherback turtles"
#    jonsen@mathstat.dal.ca, IDJ 02/10/2021
# Syntax Key:
#  - model parameters are "stochastic nodes" and are assigned with the symbol ~
#      which can be read as "is distributed as"
#  - derived parameters are "logical nodes" and are assigned with the symbol <-
#  - square brackets [ ] denote subsets, e.g., mu.alpha[1,1] denotes the 
#      1st row & first column of the matrix mu.alpha
#  - prior & sampling distributions are preceded by a "d"
#  - all variances are specified as precisions (1/variance), for convenience we
#      place priors on standard deviations and subsequently convert these 
#      parameters to precisions
#  - I(0, ) denotes censoring of values <= 0
#  for further details see http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/manual14.pdf
pi <- 3.1415     # declare pi
# Specify vague hyper-priors for all model parameters
# Hyper-priors on mean & sd of travel rates (alpha)
mu.alpha[1,1] ~ dnorm(0, 0.1)                   # day, lon
sd.alpha[1,1] ~ dunif(0, 0.1)
isd2.alpha[1,1] <- pow(sd.alpha[1,1], -2)       # convert sd to precision (1/var)
mu.alpha[1,2] ~ dnorm(0, 0.1)                   # day, lat
sd.alpha[1,2] ~ dunif(0, 0.1)
isd2.alpha[1,2] <- pow(sd.alpha[1,2], -2)       # convert sd to precision (1/var)
mu.alpha[2,1] ~ dnorm(0, 0.1)                   # night, lon
sd.alpha[2,1] ~ dunif(0, 0.1)
isd2.alpha[2,1] <- pow(sd.alpha[2,1], -2)       # convert sd to precision (1/var)
mu.alpha[2,2] ~ dnorm(0, 0.1)                   # night, lat
sd.alpha[2,2] ~ dunif(0, 0.1)
isd2.alpha[2,2] <- pow(sd.alpha[2,2], -2)       # convert sd to precision (1/var)
# Hyper-priors on mean & sd of process variability (constant b/w day & night)
mu.sigma[1] ~ dunif(0, 0.3)                     # lon
sd.sigma[1] ~ dunif(0, 0.1)
isd2.sigma[1] <- pow(sd.sigma[1], -2)           # convert sd to precision (1/var)
mu.sigma[2] ~ dunif(0, 0.3)                     # lat
sd.sigma[2] ~ dunif(0, 0.1)
isd2.sigma[2] <- pow(sd.sigma[2], -2)           # convert sd to precision (1/var)
# Derive travel rate on the plane (hierarchical means)
mu.r[1] <- pow(pow(mu.alpha[1,1], 2) + pow(mu.alpha[1,2], 2), 0.5)      # day
mu.r[2] <- pow(pow(mu.alpha[2,1], 2) + pow(mu.alpha[2,2], 2), 0.5)      # night
# Derive mu.delta (hierarchical mean of log ratio of day-night travel rates)
mu.delta <- log(mu.r[1]) - log(mu.r[2])
# Estimate Bayes predictive distributions by sampling from hyper-priors
Bpd.alpha11 ~ dnorm(mu.alpha[1,1], isd2.alpha[1,1])             # day, lon
Bpd.alpha12 ~ dnorm(mu.alpha[1,2], isd2.alpha[1,2])             # day, lat
Bpd.alpha21 ~ dnorm(mu.alpha[2,1], isd2.alpha[2,1])             # night, lon
Bpd.alpha22 ~ dnorm(mu.alpha[2,2], isd2.alpha[2,2])             # night, lat
Bpd.sigma1 ~ dnorm(mu.sigma[1], isd2.sigma[1])I(0,)             # lon
Bpd.sigma2 ~ dnorm(mu.sigma[2], isd2.sigma[2])I(0,)             # lat
Bpd.r1 <- pow(pow(Bpd.alpha11, 2) + pow(Bpd.alpha12, 2), 0.5)   # day
Bpd.r2 <- pow(pow(Bpd.alpha21, 2) + pow(Bpd.alpha22, 2), 0.5)   # night
Bpd.delta <- log(Bpd.r1) - log(Bpd.r2)	                   # day-night travel rate ratio
# Loop over J pathways
for(j in 1:J){
  # Priors for travel rate (alpha)
  alpha[1,1,j] ~ dnorm(mu.alpha[1,1], isd2.alpha[1,1])  # day, lon
  alpha[2,1,j] ~ dnorm(mu.alpha[2,1], isd2.alpha[2,1])  # day, lat
  alpha[1,2,j] ~ dnorm(mu.alpha[1,2], isd2.alpha[1,2])  # night, lon
  alpha[2,2,j] ~ dnorm(mu.alpha[2,2], isd2.alpha[2,2])  # night, lat
  r[1,j] <- pow(pow(alpha[1,1,j], 2) + pow(alpha[1,2,j], 2), 0.5)  # day
  r[2,j] <- pow(pow(alpha[2,1,j], 2) + pow(alpha[2,2,j], 2), 0.5)  # night
  delta[j] <- log(r[1,j]) - log(r[2,j])	# day-night travel rate ratio
  # Priors for process variability
  sigma[1,j] ~ dnorm(mu.sigma[1], isd2.sigma[1])I(0, )  # lon (half-Normal prior)
  isigma2[1,j] <- pow(sigma[1,j], -2)                   # convert sd to precision
  sigma[2,j] ~ dnorm(mu.sigma[2], isd2.sigma[2])I(0, )  # lat (half-Normal prior)
  isigma2[2,j] <- pow(sigma[2,j], -2)                   # convert sd to precision
  # Prior on first state = first observed location with uncertainty in dataset j
  # Ridx indexes the first and last regular time steps for dataset j
  # Oidx indexes the first observation in dataset j
  x[Ridx[j],1] ~ dt(y[Oidx[j],1], itau[Oidx[j],1], nu[Oidx[j],1])
  x[Ridx[j],2] ~ dt(y[Oidx[j],2], itau[Oidx[j],2], nu[Oidx[j],2])
  for(t in Ridx[j]:(Ridx[j+1]-2)){                     # loop over regular time intervals
    # account for earth's curvature & different day/night lengths (h)
    tmp[t,1] <-  x[t,1] + alpha[day[t],1,j] / cos(x[t,2] / 180 * pi) * h[t]
    tmp[t,2] <-  x[t,2] + alpha[day[t],2,j] * h[t]
    for(d in 1:2){                                     # loop over lon & lat
      # account for different day/night lengths (h)
      sig.tmp[t,d,j] <- isigma2[d,j] / h[t]            
      x[t+1,d] ~ dnorm(tmp[t,d], sig.tmp[t,d,j])       # process equation
      }
    }
  # Prior for scaling factor (accounts for among tag performance variability)
  logpsi[j] ~ dunif(-10, 10)
  psi[j] <- exp(logpsi[j])
  # idx indexes the observed locations in regular time step t
  for(t in (Ridx[j]):(Ridx[j+1]-1)){                   # loop over regular time intervals
    for(i in idx[t]:(idx[t+1]-1)){                     # loop over observations in each t
      for(d in 1:2){                                   # loop over lon & lat
        itau.psi[i,d] <- itau[i,d] * psi[j]            # re-scale observation error
        y[i,d] ~ dt(x[t,d], itau.psi[i,d], nu[i,d])    # Observation equation
        }
      }
    }
  }
}