## Covariate SSM used in:
## "Taking animal tracking to new depths: synthesizing horizontal-vertical movement 
#   relationships for four marine predators"
## Sophie Bestley, Ian D Jonsen, Mark A Hindell, Robert G Harcourt, Nicholas J Gales 
##  Contact sophie.bestley@aad.gov.au or ian.jonsen@mq.edu.au 

model;
{
pi <- 3.1416

## PRIORS ##
## Process uncertainty (iSigma)
Omega[1,1] <- 1
Omega[1,2] <- 0
Omega[2,1] <- 0
Omega[2,2] <- 1
iSigma[1:2,1:2] ~ dwish(Omega[,], 2)	
Sigma[1:2,1:2] <- inverse(iSigma[,])

## Turn angle (theta)
tmp[1] ~ dbeta(20, 20) 
tmp[2] ~ dbeta(10, 10) 
theta[1] <- (2 * tmp[1] -1) * pi 	## theta for directed movement (bmode[1])
theta[2] <- tmp[2] * pi * 2 		## theta for resident movement (bmode[2])

## Autocorrelation parameter/Persistence term (gamma)
gamma[2] ~ dbeta(1.5, 2)		## gamma for resident movement  
dev ~ dbeta(1,1)			## a random deviate to ensure that gamma[1] > gamma[2]
gamma[1] <- gamma[2] + dev 		## gamma for directed movement 

## Logistic function parameters relating switching probablities to covariate(s) 
beta[1] ~ dnorm(0, 0.1) 	# intercept for directed movement 
beta[2] ~ dnorm(0, 0.1)		# intercept for resident movement 

m[1,1] ~ dnorm(0, 0.1)		## coefficient to relate covariate to switching probability from directed to resident 
m[2,1] <- 0		## set to 0 if assume no r'ship b/w covariate and switching from resident to directed 

lambda[1] ~ dbeta(1, 1)		## prob of directed movement at t=1
lambda[2] <- 1 - lambda[1]	## prob of resident movement at t=1

## MODEL ##
## We use a hierarchical formulation of the Jonsen et al. (2005) switching model to fit to multiple individuals
##	simultaneously. This approach assumes that individuals share a common behavioural response 
##	but allows variability in the nature of individual responses.

for(k in 1:N){ # each individual animal
	logpsi[k] ~ dunif(-10, 10)  #prior on scaling parameter for t-dist. error variances
	psi[k] <- exp(logpsi[k])
	bmode[Xidx[k]] ~ dcat(lambda[])  # assign behav. mode for first obs

	for(i in 1:2){
		x[Xidx[k],i] ~ dt(y[Yidx[k],i],itau2[Yidx[k],i],nu[Yidx[k],i]) # 'true' position x at 1st regular timestep t=1		
	}

	# assume simple random walk to estimate 2nd regular position
	x[(Xidx[k]+1),1:2] ~ dmnorm(x[Xidx[k],], iSigma[,])

	## Transition equation ##
	## The transition equation models the 'true' location and unobserved behavioural states as a function of their previous
	##	 states (ie. a first order Markov process)
	## The transition equation has been modified from Jonsen et al. 2005 Ecology 86:2874-2880 by including 2 logistic
	## 	functions to model the switching probabilities as a function of covariate(s). Whereas the Jonsen
	##	et al. (2005) model estimated 2 constant switching probabilities, this model estimates switching probabilities at
	##	time step (phi[t,]).  
  ## This implements the one-way switch influence single covariate model from Bestley et al 2013 PRSB 280: 20122262 
  #   (where covariates are not included in the determination of the reverse switch probabilities)
	
  for(t in (Xidx[k]+1):(Xidx[k+1]-2)){ # regular timesteps t 
	  # covariate influences Pr(directed) at time t directly
	  mu[t,1] <- exp(beta[1] + m[1,1] * covar[t,1] )    # (mu1) staying directed
	  mu[t,2] <- exp(beta[2] + m[2,1] * covar[t,1] )    # (mu2) switching to directed (in one-way model m[2,1] is set to zero)

	  phi[t,1] <- mu[t,1] / (1 + mu[t,1]) # pr(directed|directed)
	  phi[t,2] <- mu[t,2] / (1 + mu[t,2]) # Pr(directed|resident)
	
	  eta[t,1] <- phi[t,bmode[t-1]] # pr(directed|directed) or pr(directed|resident)
	  eta[t,2] <- 1 - phi[t,bmode[t-1]] # Pr(resident|directed) or Pr(resident|resident)
	
	  bmode[t] ~ dcat(eta[t,])

	  Tdx[t,1] <- cos(theta[bmode[t]]) * (x[t,1] - x[t-1,1]) + sin(theta[bmode[t]]) * (x[t,2] - x[t-1,2])
 	  Tdx[t,2] <- -sin(theta[bmode[t]]) * (x[t,1] - x[t-1,1]) + cos(theta[bmode[t]]) * (x[t,2] - x[t-1,2]) 

	  x.mn[t,1] <- x[t,1] +  Tdx[t,1] * gamma[bmode[t]]
	  x.mn[t,2] <- x[t,2] +  Tdx[t,2] * gamma[bmode[t]]
  
	  # predict next location (with error)
	  x[t+1,1:2] ~ dmnorm(x.mn[t,], iSigma[,])	
	  }

	## estimate final bmode  
	mu[Xidx[k+1]-1,1] <- exp(beta[1] + m[1,1] * covar[Xidx[k+1]-1,1] )
	mu[Xidx[k+1]-1,2] <- exp(beta[2] + m[2,1] * covar[Xidx[k+1]-1,1] )

	phi[Xidx[k+1]-1,1] <- mu[Xidx[k+1]-1,1] / (1 + mu[Xidx[k+1]-1,1])
	phi[Xidx[k+1]-1,2] <- mu[Xidx[k+1]-1,2] / (1 + mu[Xidx[k+1]-1,2])
	
	eta[Xidx[k+1]-1,1] <- phi[Xidx[k+1]-1,bmode[Xidx[k+1]-2]]
	eta[Xidx[k+1]-1,2] <- 1 - phi[Xidx[k+1]-1,bmode[Xidx[k+1]-2]]
	
	bmode[Xidx[k+1]-1] ~ dcat(eta[Xidx[k+1]-1,])

	## Robust measurement equation ##
	## The measurement equation relates the estimated 'true' location states to the observed locations
	## As in Jonsen et al. (2005), we use a robust likelihood (t-distribution) to model the observation because they are
	##	prone to occasionally large errors. We use the independently estimated Argos error variances (1/itau2[i,q]) and
	##	degrees of freedom (nu[i,q]) from Jonsen et al. (2005).

	for(t in (Xidx[k]+1):(Xidx[k+1]-1)){  # loops over regular time intervals (t)
	  for(i in idx[t]:(idx[t+1]-1)){  # loops over observed locations within interval t
	    for(q in 1:2){       # loops over longitude and latitude
	      zhat[i,q] <- (1 - j[i]) * x[t-1,q] + j[i] * x[t,q]
	      itau2.psi[i,q] <- psi[k] * itau2[i,q]  # re-scale error variances to account for variation in accuracy among tags
	      y[i,q] ~ dt(zhat[i,q], itau2.psi[i,q], nu[i,q])
	      } #q
	    } #i
	  }	#t
	} #k
} #end