# Stage-2 state-space model from:
#	Eckert, S., J. Moore, D. Dunn, R. Sagarminaga, K. Eckert, P. Halpin
#	"Hierarchical state-space models of loggerhead sea turtle (Caretta caretta)
#	movement in relation to turtle size and oceanographic features
#	in the western Mediterranean Sea"
#	jemoore@duke.edu, 20 July 2007

# Based on the hierarchical switch model from:
#	Morales, J., D. Haydon, J. Frair, K. Holsinger, and J. Fryxell.	
#	"Extracting more out of reloction data: building movement models
#	as mixtures of random walks"
#	Ecology 85:2436-2445 



model{

# Mode 1 is slow, Mode 2 is fast
# For hyper-parameters, mean tortousity (rho.mu) is constrained to be greater for slow mode (rho.mu[1]<rho.mu[2])
# mean turn angle is assumed to be zero for all turtles and for both modes

# Linear model (logit function) for switch probabilities is:
#	B0 + Bsize(size) + Bdepth(depth) + Bmsla(msla) + BsizeMsla(sizeMsla)
# There are two switch functions estimated:
# 1. Staying in state 1 at t, given state 1 at t-1 (1 minus this value = switching from 1 to 2),
# 2. Switching to state 1 at t, given state 2 at t-1 (1 minus this value = staying in state 2, given 2).

# Covariate data were standardized to z-scores


Pi <- 3.14159265359		## define pi
nPi <- Pi * -1



## PRIORS FOR FIXED TURTLE-SIZE EFFECTS

bSize[1] ~ dnorm(0, 0.001)		# effect of size on switching to state 1 at t, given 1 or 2 at t-1
bSize[2] ~ dnorm(0, 0.001)	

bSizeMsla[1] ~ dnorm(0, 0.001)		# interaction of size and msla, on switching
bSizeMsla[2] ~ dnorm(0, 0.001)



## PRIORS FOR HYPER-PARAMETERS


# hyperpriors for environmental effects (slope coefficients)

B0.mu[1] ~ dnorm(0,0.001)		# Intercept; mean logit for switching from mode 1 to 1
B0.sig[1] ~ dnorm(0,0.001) I(0,)
B0.tau[1] <- 1/(B0.sig[1]*B0.sig[1])
B0.mu[2] ~ dnorm(0,0.001)		# mean logit for switching from mode 2 to 1
B0.sig[2] ~ dnorm(0,0.001) I(0,)
B0.tau[2] <- 1/(B0.sig[2]*B0.sig[2])

bMsla.mu[1] ~ dnorm(0,0.001)
bMsla.sig[1] ~ dnorm(0,0.001) I(0,)
bMsla.tau[1] <- 1/(bMsla.sig[1]*bMsla.sig[1])
bMsla.mu[2] ~ dnorm(0,0.001)
bMsla.sig[2] ~ dnorm(0,0.001) I(0,)
bMsla.tau[2] <- 1/(bMsla.sig[2]*bMsla.sig[2])

bDepth.mu[1] ~ dnorm(0,0.001) I(-3,3)
bDepth.sig[1] ~ dnorm(0,0.001) I(0,)
bDepth.tau[1] <- 1/(bDepth.sig[1]*bDepth.sig[1])
bDepth.mu[2] ~ dnorm(0,0.001) I(-3,3)
bDepth.sig[2] ~ dnorm(0,0.001) I(0,)
bDepth.tau[2] <- 1/(bDepth.sig[2]*bDepth.sig[2])


# hyperpriors for movemement parameters for each behavioral mode

# overall mean rates for each mode; mode 1 is constrained to be slower 
rateMean.mu[1] ~ dnorm(0,0.001) I(0,)	
rateMean.muDiff ~ dnorm(0,0.001) I(0,)
rateMean.mu[2] <- rateMean.mu[1]+rateMean.muDiff

# individual variation in mode-specific movement rates					
rateMean.tau1 ~ dnorm(0,0.001) I(0,)		
rateMean.var[1] <- 1/rateMean.tau1			# variance across turtles of mode 1 mean-rate
rateMean.sig[1] <- pow( rateMean.var[1],0.5)		# sd across turtles of mode 1 mean-rate

rateMean.tauDiff ~ dnorm(0,0.001) I(0,)			# diff in variance between the two modes
rateMean.varDiff <- 1/rateMean.tauDiff
rateMean.var[2] <- rateMean.var[1] + rateMean.varDiff	# variance across turtles of mode 2 mean-rate; assume greater than for mode 1
rateMean.sig[2] <- pow( rateMean.var[2] , 0.5)		# sd across turtles of mode 2 mean-rate	

# mean variance in movement rates for each mode (hyper-parm for indiv-level variance)
# expressed here as precision (1/var)
rateTau.mu[2] ~ dnorm(0,0.001) I(0,)	
rateTau.muDiff ~ dnorm(0,0.001) I(0,)
rateTau.mu[1] <- rateTau.mu[2] + rateTau.muDiff		# constrain mode 1 to have higher precision (lower variance) in movement rates	

# individual variation in mode-specific variance for rate (hyper-parm for indiv-level varianc)
rateTau.tau2 ~ dnorm(0,0.001) I(0,)	
rateTau.tauDiff ~ dnorm(0,0.001) I(0,)	# i.e., some turtles are more variable in their mode-specific rates

# mean tortuosity of movement; assume greater for mode 1

rho.mu[1] ~ dbeta(1,1) 		# mean tortuosity for mode 1
r[1] ~ dnorm(0,0.001) I(0,)	# beta parameter	
s[1] <- r[1]/rho.mu[1] - r[1]	# second beta parameter, calucated from rho.mu and r
rho.var[1] <- r[1]*s[1]/( (r[1]+s[1])*(r[1]+s[1])*(r[1]+s[1]+1) )	# variance in tortuosity for mode 1, calculated from r and s
rho.sig[1] <- pow(rho.var[1], 0.5)	# sd in mode 1 tortuosity

rhoDiff.mu ~ dnorm(0,0.01) I(0,)	# difference in mean tortusosity between modes

rho.mu[2] <- rho.mu[1]+rhoDiff.mu	# mean tortuosity for model 2; is less tortuous than mode 1
r[2] ~ dnorm(0,0.001) I(0,)	# beta parameter	
s[2] <- r[2]/rho.mu[2] - r[2]	# second beta parameter
rho.var[2] <- r[2]*s[2]/( (r[2]+s[2])*(r[2]+s[2])*(r[2]+s[2]+1) )	# variance in tortuosity for mode 2
rho.sig[2] <- pow( rho.var[2], 0.5)	# sd
rhoDiff.sig <- pow( (rho.var[1]+rho.var[2]) , 0.5)	# difference in SD of tortuosity between the two modes



## LOOP THROUGH EACH TURTLE

for (k in 1:K){		# k is number of turtles in analysis
	
	# TURTLE-SPECIFIC PRIORS FOR BEHAVIORAL MODE 1 AND 2
	
		# priors for environmental effects; normally distributed with mean and variance described by hyper-parameters
		
		B0[k,1] ~ dnorm(B0.mu[1], B0.tau[1])			# Intercept; mean logit for switching to mode 1 from 1
		B0[k,2] ~ dnorm(B0.mu[2], B0.tau[2])			# Intercept; mean logit for switching to mode 1 from 2
		
		bMsla[k,1] ~ dnorm(bMsla.mu[1], bMsla.tau[1])		# coefficient for MSLA 
		bMsla[k,2] ~ dnorm(bMsla.mu[2], bMsla.tau[2]) 
		
		bDepth[k,1] ~ dnorm(bDepth.mu[1], bDepth.tau[1]) 	# coefficient for depth
		bDepth[k,2] ~ dnorm(bDepth.mu[2], bDepth.tau[2]) 
	
	
		# priors characterizing movement parms for the two modes 

		# each turtle's mean rate (mode-specific)
		km.mu[k,1] ~ dnorm(rateMean.mu[1], rateMean.tau1) I(0,)	
		km.muDiff[k] ~ dnorm(rateMean.muDiff, rateMean.tauDiff) I(0,)
		km.mu[k,2] <- km.mu[k,1]+km.muDiff[k]
		
		# each turtle's variance in rate (mode-specific)
		km.tau[k,2] ~ dnorm(rateTau.mu[2], rateTau.tau2) I(0,)	
		km.tauDiff[k] ~ dnorm(rateTau.muDiff, rateTau.tauDiff) I(0,)
		km.tau[k,1] <- km.tau[k,2] + km.tauDiff[k]
		km.sig[k,2] <- pow(1/km.tau[k,2], 0.5)
		km.sig[k,1] <- pow(1/km.tau[k,1], 0.5)
	
		# each turtle's tortuosity
		rho[k,1] ~ dbeta(r[1],s[1])
		rho[k,2] ~ dbeta(r[2],s[2]) I(0.000005, ) # ad-hoc lower bound added (arbitrarily small number) to avoid bracket-slice failures
		rhoDiff[k] <- rho[k,2]-rho[k,1]
	
		# each turtle's avg turn angle
		mu[k,1] <- 0	# assume that mean turn angle for each movement type (for each turtle) is zero
		mu[k,2] <- 0
		
	# INITIAL STATE FOR EACH TURTLE
	
		phi[k,1] ~ dunif(0,1)	# probability of being in bmode 1 at time 1
		phi[k,2] <- 1-phi[k,1]
		bmode[ Ridx[ Kidx[k] ] ] ~ dcat(phi[k,])	
	
	# LOOP THROUGH LOCATIONS FOR EACH TURTLE
	
	# First location for each turtle is indexed by Ridx
	# start with 2nd loc, and ending w 2nd to last
   for ( t in (Ridx[k]+1):(Ridx[k+1]-2) ) {	
		
		msla[t] ~ dnorm(0,1)				# priors for missing covariate data
		depth[t] ~ dnorm(0,1)

		
		# MODEL Q(switch prob) as a function of covariates

		logit.Q[t,1] <- exp( B0[k,1] + bDepth[k,1]*depth[t] + bMsla[k,1]*msla[t] + bSize[1]*size[k] + bSizeMsla[1]*size[k]*msla[t] )
		logit.Q[t,2] <- exp( B0[k,2] + bDepth[k,2]*depth[t] + bMsla[k,2]*msla[t]  + bSize[2]*size[k] + bSizeMsla[2]*size[k]*msla[t])	
		
		Q[t,1] <- logit.Q[t,1]/(1+logit.Q[t,1])		# tranform logits into switch probabilities
		Q[t,2] <- logit.Q[t,2]/(1+logit.Q[t,2])
		
		alpha[t,1] <- Q[t,bmode[t-1]]	# probability of being in bmode 1 at times > 1
		alpha[t,2] <- 1-alpha[t,1]
   	 	bmode[t] ~ dcat(alpha[t,])  	# estimated bmode is a Bernoulli draw 

		## Likelihood for steps
		kmDay[t] ~ dnorm(km.mu[k,bmode[t-1]], km.tau[k,bmode[t-1]]) I(0,)	# Normal distribution for step length

		## Likelihood for turns.  
		## use the “ones” trick (see WinBUGS manual) to sample from the Wrapped Cauchy distribution

		ones[t] <- 1
		ones[t] ~ dbern(wC[t])
		## pdf for Wrapped Cauchy distribution, divided by 500 (arbitrary) to ensure that wC[t] will be less than one
		wC[t] <- ( 1/(2*Pi)*(1-rho[k,bmode[t-1]]*rho[k,bmode[t-1]])/(1+rho[k,bmode[t-1]]*rho[k,bmode[t-1]]-2*rho[k,bmode[t-1]]*cos(deltRad[t]-mu[k,bmode[t-1]])) )/500
	
# END TURTLE-SPECIFIC LOCATION LOOP
  }


# END LOOPING THROUGH TURTLES
}
	

# END FULL MODEL LOOP
}