# Code for the Columbia river example discussed in the paper

#    MUTHUKUMARANA S. SCHWARZ CJ, Swartz, TB (2007)
#    Bayesian analysis of mark-recapture data with
#    travel-time-dependent survival probabilities.
#    Canadian Journal of Statistics, **, ****-****.


model
{     

#  N is the no.of fish
#  m is the no.of listening lines

#  LT[i,j]  is the log of T[i,j]  i.e  log of cumulative travel time of ith fish at jth listening line 
#  l_LT[i,j] and u_LT[i,j] are lower and upper bound for LT[i,j]  
#  mu[j] is the mean of logT[,j]
#  l_mu[j] and u_mu[j] are lower and upper bound for mu[j] 
#  varcov is the sigma, variance covariance matrix of log T's
#  mean[j] is the mean of t[,j], i.e. the mean travel time between line j and line j+1
#  sd[j] is the standard deviation of t[,j], the standard deviation of the travel times between line j and line j+1
#  corr is the correlations of t's across the listening lines

#  phi[i,j] is the survival probability of ith fish between  j-1th and jth  listening line
#  q[j] is the daily survival probability between  (j-1)th and jth  listening line

#  p[j] is the detection probability at  jth  listening line

#  c[i,j] is the capture history  of ith fish at jth  listening line
#  s[i,j] is the survival state of ith fish at jth  listening line

#  P is the no.of parameters.
#     This is equal to the number of p's (capture rates), q's (daily survival rates),
#     mu's (mean time between lines), the sigma matrix (variance and covariances on travel time),
#     and the number of missing survival and travel time values (which must be estimated
#     by WinBugs
#  nstar is the "sample size". Following the usual capture-recapture conventions,
#     this should equal the number of "releases" of fish.
#     For example, capture history 10110 (with the initial release listed in the history)
#     counts as a sample size of 3, the 101 + 11 + 10 parts of the history
#     In this program, the initial releases are not recorded in the history, but this is available
#     from N.

#  loglik is the log likelihood of the data 
#  bic is the BIC of the model

# Calculating nstar - the sample size
   for(i in 1:m-1)
   {
      n[i]<-sum(c[,i])
   }
   nstar<-N+sum(n[])

# Calculating P - the number of parameters

   for(i in 1:N)
   {
      ss[i,m]<-c[i,m]
      for(j in 2:m)
      {
        ss[i,m+1-j]<-step(c[i,m+2-j]+c[i,m+1-j]-.5)
      }
   }
#  In the formula below, 
#     N*m-sum(ss[,]) is the no. of missing survival states.
#     3*m+m*(m+1)/2  is the no. of parameters from mu's, p's,q's and sigma matrix.
#     N*m-sum(c[,])  is the no. of missing T[i,j]'s.
   P<-N*m-sum(ss[,])+3*m+m*(m+1)/2 +N*m-sum(c[,])


# Setting up boundary  values

    for(i in 1:N)
    {  
        u_LT[i,m]<- 50        # upper and lower bounds on log-scale
        l_LT[i,1]<--10
      
        for(j in 1:m-1)
        {     
           u_LT[i,j]<-LT[i,j+1]          
           l_LT[i,j+1]<-LT[i,j]
         }
     }

#Calculating survival probabilities
     for(i in 1:N)
     {       
        phi[i,1]<-pow(q[1],exp(LT[i,1]))
        for(j in 2:m)
        {     
           phi[i,j]<-s[i,j-1]*pow(q[j],exp(LT[i,j])-exp(LT[i,j-1]))
        }
     }

# Evaluating likelihood 
     for(i in 1:N)  
     {
     #log of T is constrained MVN as defined in (5)
         LT[i,1:m] ~ dmnorm(mu[] , G[,])I(l_LT[i,],u_LT[i,])
          
         for(j in 1:m)     
           {
            s[i,j]~dbern(phi[i,j])     
            term[i,j] <- s[i,j]*p[j]
            c[i,j]~dbern(term[i,j]) 
                                      
            term1[i,j] <- s[i,j]*log(pow(p[j],c[i,j])*pow(1-p[j],1-c[i,j]))     
           }
     }

# Prior distributions as described in model development

# Priors for p's and q's
     
     for(j in 1:m)     
     {     
       q[j]~dbeta(1,1)  # prior for daily survival rate
       p[j]~dbeta(1,1)  # prior for catchability at each listening line
     }
          
# setting up boundary values for mu
     
     u_mu[m]<-50     
     l_mu[1]<- -10
          
     for(j in 1:m-1)
     {          
        u_mu[j]<-mu[j+1]          
        l_mu[j+1]<-mu[j]
     }
          
# constrained prior for mu as described in model development
     mu[1:m] ~ dmnorm(g0[] , gv[,])I(l_mu[],u_mu[])

    for(j in 1:m)  # set up mean and covariance matrix 
    {     
       g0[j] <- 0
       gv[j,j] <- .0001     
    }

    for(i in 1:m-1)
    {     
        for(j in i+1:m)
        {     
           gv[i,j]<-0
           gv[j,i]<-0
        }
     }

# Prior for sigma as described in model development
     G[1:m,1:m] ~ dwish(R[,],m)

     for(j in 1:m)
     {     
        R[j,j] <- 1/m     
     }

    for(i in 1:m-1)
    {     
       for(j in i+1:m)
       {     
           R[i,j]<-0
           R[j,i]<-0
        }
    }

    varcov[1:3,1:3] <- inverse(G[,])  

#  calculating BIC

     for(i in 1:N)     
     {          
        term2[i,1] <- log(pow(phi[i,1],s[i,1])*pow(1-phi[i,1],1-s[i,1]))     
        like[i,1]<-term1[i,1]+term2[i,1]
        for(j in 2:m)
        { 
          term2[i,j] <- s[i,j-1]*log(pow(phi[i,j],s[i,j])*pow(1-phi[i,j],1-s[i,j]))
          like[i,j]<-term1[i,j]+term2[i,j]
        }
     }
     
     sigma1[1]<-G[1,1]
     sigma1[2]<-G[2,1]
     sigma1[3]<-G[3,1]     
     
     sigma2[1]<-G[1,2]
     sigma2[2]<-G[2,2]
     sigma2[3]<-G[3,2]     
     
     sigma3[1]<-G[1,3]
     sigma3[2]<-G[2,3]
     sigma3[3]<-G[3,3]
     

     for(i in 1:N)
     {                         
        xm[i,1]<- LT[i,1]-mu[1]
        xm[i,2]<- LT[i,2]-mu[2]
        xm[i,3]<- LT[i,3]-mu[3]
          
        pro[i,1]<- inprod(xm[i,],sigma1[])
        pro[i,2]<- inprod(xm[i,],sigma2[])
        pro[i,3]<- inprod(xm[i,],sigma3[])
          
        final[i]<- -.5*inprod(pro[i,],xm[i,]) + .5*logdet(G[,])
     }

     loglik<-sum(like[,])+sum(final[])
     bic<- -2*(sum(like[,])+sum(final[]))+P*log(nstar)

 # calculating T[i,j]
     for(i in 1:N)
     {     
       for(j in 1:m)
       {      
         T[i,j]<-exp(LT[i,j])
       }
     }

# calculating mean and SD of  t's
     for(i in 1:N)
     {     
        t[i,1]<-T[i,1]     
        for(j in 2:m)
        {      
           t[i,j]<-T[i,j]-T[i,j-1]
        }
     }
     for(j in 1:m)
     {      
        mean[j]<-mean(t[,j])
        sd[j]<-sd(t[,j])
        var[j]<-sd[j]*sd[j]
        sdn[j]<-sqrt((N-1)*var[j]/N)
     }

# compute the correlation of the tavel times

     for(i in 1:N)
     {     
        crossprod1[i]<-(t[i,1]-mean[1])*(t[i,2]-mean[2])
        crossprod2[i]<-(t[i,1]-mean[1])*(t[i,3]-mean[3])
        crossprod3[i]<-(t[i,2]-mean[2])*(t[i,3]-mean[3])
     }
     crossprod11<-sum(crossprod1[])
     crossprod22<-sum(crossprod2[])
     crossprod33<-sum(crossprod3[])

     corr[1]<-crossprod11/((N-1)*sdn[1]*sdn[2])
     corr[2]<-crossprod22/((N-1)*sdn[1]*sdn[3])
     corr[3]<-crossprod33/((N-1)*sdn[2]*sdn[3])

} # end of model