# Winbugs code for the model published in
# Antoniazza M, Korner-Nievergelt F & Keller V (2012) Les mouvements des Grands Cormorans Phalacrocorax carbo
# bagueés dans la colonie du Fanel, lac de Neuchâtel. Nos Oiseaux, 59, 11-22.


#---------------------------------------------------------------------------------
model{
####################################################
# Parameters:
# phi: survival probability
# mAB: movement probability from site A to site B
# mAC: movement probability from site A to site C
# mBA: movement probability from site B to site A
# mBC: movement probability from site B to site C
# mCA: movement probability from site C to site A
# mCB: movement probability from site C to site B
# pA: recapture probability at site A
# pB: recapture probability at site B
# pC: recapture probability at site C
###################################################
# States and observations:
# STATES (S)
# 1 alive at A
# 2 alive at B
# 3 alive at C
# 4 dead
# OBSERVATIONS (O)
# 1 seen at A
# 2 seen at B
# 3 seen at C
# 4 not seen
###################################################

# priors for movement probabilities
# for first year of life
mAB[1,1] ~ dbeta(1,1)
invmAB[1,1] <- 1-mAB[1,1]
mAC[1,1] ~ dunif(0, invmAB[1,1])

mBA[1, 1] <-0
invmBC[1, 1] <- 0
mBC[1, 1] <-0

mCB[1, 1] <-0
invmCB[1, 1] <- 0
mCA[1, 1]<-0

for(sais in 2:3){
  mAB[1,sais] ~ dbeta(1,1)
  invmAB[1,sais] <- 1-mAB[1,sais]
  mAC[1,sais] ~ dunif(0, invmAB[1,sais])

  mBA[1, sais] ~ dunif(0, invmBC[1,sais])
  invmBC[1, sais] <- 1-mBC[1,sais]
  mBC[1, sais] ~ dbeta(1,1)

  mCB[1, sais]  ~ dbeta(1,1)
  invmCB[1, sais] <- 1-mCB[1,sais]
  mCA[1, sais]~dunif(0, invmCB[1,sais])
}

mAB[1,4] <-0
invmAB[1,4] <-0
mAC[1,4] <-0

mBA[1, 4] ~ dunif(0, invmBC[1,4])
invmBC[1, 4] <- 1
mBC[1, 4] <- 0

mCB[1, 4]  ~ dbeta(1,1)
invmCB[1, 4] <- 1-mCB[1,4]
mCA[1, 4]~dunif(0, invmCB[1,4])


# for older birds
    mAB[2, 1] ~ dbeta(1, 1)
    invmAB[2, 1] <- 1-mAB[2, 1]
    mAC[2, 1]~ dunif(0, invmAB[2, 1])

    mBA[2, 1] <- 0
    invmBC[2, 1] <-  1
    mBC[2, 1] ~ dunif(0, invmBC[2, 1])

    mCB[2, 1] <- 0
    invmCB[2, 1] <- 0
    mCA[2, 1] <- 0


  for(sais in 2:3){
    mAB[2, sais] ~ dbeta(1, 1)
    invmAB[2, sais] <- 1-mAB[2, sais]
    mAC[2, sais]~ dunif(0, invmAB[2, sais])

    mBA[2, sais] ~ dbeta(1, 1)
    invmBC[2, sais] <-  1-mBA[2, sais]
    mBC[2, sais] ~ dunif(0, invmBC[2, sais] ) # constrain the sum of transitions from state 3 to one

    mCB[2, sais] ~ dbeta(1, 1)
    invmCB[2, sais] <- 1-mCB[2, sais]
    mCA[2, sais]~dunif(0, invmCB[2, sais])
    }
mAB[2, 4] <-0
invmAB[2, 4] <- 0
mAC[2, 4]<-0

mBA[2, 4] ~ dbeta(1, 1)
invmBC[2, 4] <-  1
mBC[2, 4] <- 0 # constrain the sum of transitions from state 3 to one

mCB[2, 4] ~ dbeta(1, 1)
invmCB[2, 4] <- 1-mCB[2, 4]
mCA[2, 4]~dunif(0, invmCB[2, 4])


a0~dunif(-10, 10)
bA~dunif(-10, 10) # intercept of resighting probability constraint
bB~dunif(-10, 10) # intercept of resighting probability constraint
bC~dunif(-10, 10) # intercept of resighting probability constraint
age.phi~dunif(-5, 5)  # age effect on survival

ageseas2~dunif(-5, 5)  # interaction age x saison
ageseas3~dunif(-5, 5)
ageseas4~dunif(-5, 5)

for(k in 1:3){
  sais.phi[k]~dunif(-5, 5)
  sais.p[k]~dunif(-5, 5)
  }

#  constraints   (assumption: movement at the beginning of the time interval, survival conditional on movement)
for(i in 1:nind){
  for(j in first[i]:(nocc-1)){
    logit(phi[i,j]) <- a0 + age.phi*agemat23[i,j+1] + sais.phi[1]*season2[j+1] +
        sais.phi[2]*season3[j+1]+ sais.phi[3]*season4[j+1]

    logit(pA[i,j]) <- bA  + sais.p[1]*season2[j+1] + sais.p[2]*season3[j+1]+ sais.p[3]*season4[j+1]

    logit(pB[i,j]) <- bB  + sais.p[1]*season2[j+1] + sais.p[2]*season3[j+1]+ sais.p[3]*season4[j+1]

    logit(pC[i,j]) <- bC  + sais.p[1]*season2[j+1] + sais.p[2]*season3[j+1]+ sais.p[3]*season4[j+1]

    }
  }


# Define parameters
for (i in 1:nind){
   # define probabilities of state S(t+1) given S(t)
   for (j in first[i]:(nocc-1)){  # loop over time
      # first index = states at time t-1, last index = states at time t
      ps[1,i,j,1] <- phi[i,j] * (1-mAB[agemat[i,j], season[j]]-mAC[agemat[i,j], season[j]])
      ps[1,i,j,2] <- phi[i,j] * mAB[agemat[i,j], season[j]]
      ps[1,i,j,3] <- phi[i,j] * mAC[agemat[i,j], season[j]]
      ps[1,i,j,4] <- 1-phi[i,j]
      ps[2,i,j,1] <- phi[i,j] * mBA[agemat[i,j], season[j]]
      ps[2,i,j,2] <- phi[i,j] * (1-mBA[agemat[i,j], season[j]]-mBC[agemat[i,j], season[j]])
      ps[2,i,j,3] <- phi[i,j] * mBC[agemat[i,j], season[j]]
      ps[2,i,j,4] <- 1-phi[i,j]
      ps[3,i,j,1] <- phi[i,j] * mCA[agemat[i,j], season[j]]
      ps[3,i,j,2] <- phi[i,j] * mCB[agemat[i,j], season[j]]
      ps[3,i,j,3] <- phi[i,j] * (1-mCA[agemat[i,j], season[j]]-mCB[agemat[i,j], season[j]])
      ps[3,i,j,4] <- 1-phi[i,j]
      ps[4,i,j,1] <- 0
      ps[4,i,j,2] <- 0
      ps[4,i,j,3] <- 0
      ps[4,i,j,4] <- 1

      # define probabilities of O(t) given S(t)
      # first index = states at time t, last index = observations at time t
      po[1,i,j,1] <- pA[i,j]
      po[1,i,j,2] <- 0
      po[1,i,j,3] <- 0
      po[1,i,j,4] <- 1-pA[i,j]
      po[2,i,j,1] <- 0
      po[2,i,j,2] <- pB[i,j]
      po[2,i,j,3] <- 0
      po[2,i,j,4] <- 1-pB[i,j]
      po[3,i,j,1] <- 0
      po[3,i,j,2] <- 0
      po[3,i,j,3] <- pC[i,j]
      po[3,i,j,4] <- 1-pC[i,j]
      po[4,i,j,1] <- 0
      po[4,i,j,2] <- 0
      po[4,i,j,3] <- 0
      po[4,i,j,4] <- 1
      } #close j
   } #i

# State-space model likelihood
for (i in 1:nind){
   for (t in (first[i]+1):nocc){  # loop over time
      # State equation: draw S(t) given S(t-1)
      alive[i,t] ~ dcat(ps[alive[i,t-1],i,t-1,])

      # Observation equation: draw O(t) given S(t)
      indhist[i,t] ~ dcat(po[alive[i,t],i,t-1,])
      } #t
   } # i
}