`SCRmultinomial.fn` <-
function(){
# This R script generates spatial capture recapture data under a
# multinomial observation model, defines the WinBUGS model, and then
# fits the model to the simulated data in WinBUGS. The model is similar
# to that described in Chapter 7 but this is from Gardner et al. (2009; Ecology).
#
#
# The model assumes the individual probability of capture is:
# pi[i]<-po*exp(beta/-E[i])
# where E[i] is a function of trap locations and individual activity center
# and a parameter "sigma" which is estimated along with beta and p0 and the
# number of individuals (population size)
# The total number of captures from J samples is binomial
# Nsum[i]~dbin(mu[i], J)
# Conditional on capture, the trap-of-capture is a multinomial random
# variable which is described in WinBUGS as a categorical random variable
# using the distribution dcat()
# Specifically:
#       Hvec[n,1] ~ dcat(PI[Hvec[n,2], ])
# The data are stored in a compact form where Hvec[n,1] is the trap-of-capture
# and Hvec[n,2] is the individual that was captured.
# The WinBUGS implementation uses data augmentation (Royle et al. 2007; J. Comp. and Graph. Stats)

library("spdep")
library("R2WinBUGS")

#set the number of sample occassions, M, and N (the number of animals in the bounding box)
J=10  # number reps
M=180  # number excess zeros used for data augmentation
N=80   # population size
#Create a 10x10 grid
x=4:13
y=4:13
locs=expand.grid(x,y)
ntraps=100
#Set sigma and po parameters for the simulated data
sigma=1.5
po=.4
#Set larger area boundaries (0,0) to (18,18) for "buffer" area
llx=0
upx=18
#Simulate N home range centers
bearx=runif(N, 0,upx)
beary=runif(N, 0,upx)
#Make the Dij matrix of distances from home center to trap locations
D=matrix(0, nrow=length(locs[,1]), ncol=N)
for(i in 1:N){
      D[ ,i] = spDistsN1(as.matrix(locs), c(bearx[i], beary[i]))
}
#Make the K matrix
k=matrix(0, nrow=length(locs[,1]), ncol=N)
k=exp(-D^2/sigma)
#Now, calculate PIij matrix (conditional [on n] capture probability)
PI=matrix(0, nrow=length(locs[,1]), ncol=N)
for(j in 1:N) {
for(i in 1:length(locs[,1])){
     PI[i,j] = k[i,j]/sum(k[,j])
}
}
#to make the encounter history for each simulated bear
#first get at probability of encounter
prob=NULL
for(j in 1:N){
prob[j] = po*exp(1/(-1*sum(k[,j])))
}
#Now, make enounter history for each bear over J sample periods
nij=matrix(0, nrow=J, ncol=N)
for(j in 1:N) {
   nij[ ,j]=rbinom(J, 1, prob[j])
}
nij=t(nij)
hij=matrix(0, nrow=N, ncol=J)
for(i in 1:N){
a=sum(nij[i,])
b=sample(1:100, a, replace = TRUE, prob = PI[,i])
c=which(nij[i,] > 0)
if(length(c) >0) {hij[i,c]=b}
}
if(TRUE){   #set to TRUE, if you want to make some plots of the simulated data
plot(c(llx,upx), c(llx, upx), typ='n')
points(locs)
text(bearx, beary, col='red')
hold=which(apply(nij, 1, sum) > 1)
for(mm in 1:10) {
bno=hold[mm]
for(i in 1:J){
a=hij[bno,]
lines(c(bearx[bno],locs[a[i],1]), c(beary[bno],locs[a[i],2]))
}
}
rowSums(nij)
}

##########Set up the variables for WinBUGS
Nsumm=as.vector(apply(nij, 1, sum))
aug=rep(0, M-N)
##### Ncap = number of captuers for each individual
Ncap=as.vector(c(Nsumm, aug))
trapmat=as.matrix(locs)
##### Totcaps = total number of captures, over all individuals
Totcaps=sum(Ncap)
indx=matrix(0, nrow=N, ncol=J)
for(i in 1:N) {
a=which(nij[i,] > 0)
indx[i,a]=i
}
hvec=as.matrix(cbind(c(t(hij)), c(t(indx))))
a=which(hvec[,1] >0)
hvec=hvec[a,]
hvec=as.matrix(hvec)
sink("multinomialSCR.txt")
cat("
model {
for (i in 1:M){
 
   for(j in 1:ntraps) {
  D2[i,j] <- pow(SX[i]-trapmat[j,1], 2) + pow(SY[i]-trapmat[j,2],2)
  K[i,j] <- exp(-D2[i,j]/sigma)
  PI[i,j] <- K[i,j]/E[i]
  }
 
 z[i]~dbern(psi)
 E[i] <- sum(K[i, ])
 Ncap[i]~dbin(mu[i], J)
 mu[i]<-pi[i]*z[i]
 pi[i]<-po*exp(beta/-E[i])
 SX[i]~dunif(0, 18)
 SY[i]~dunif(0, 18)
 }
for(n in 1:Totcaps){
Hvec[n,1] ~ dcat(PI[Hvec[n,2], ])
}
beta~dunif(0,10)
psi~dunif(0,1)
po~dunif(0, 2)
sigma~dunif(0,10)
N<-sum(z[1:M])
}
",fill=TRUE)
sink()

data1<-list(M=M,Hvec=hvec, Ncap=Ncap, Totcaps=Totcaps, ntraps=ntraps, trapmat=trapmat, J=J)
params = list('po', 'sigma', 'psi', 'N','beta')
po1=runif(1,.1,1)
sigma1=runif(1,0,5)
psi=.5
z=as.vector(rbinom(M, 1, .5))
SX=as.vector(runif(M,3,15))
SY=as.vector(runif(M,3,15))
inits =  function() {list(po=po1, z=z, SX=SX, SY=SY, psi=psi, sigma=sigma1) }
fit = bugs(data1, inits, params, "multinomialSCR.txt", debug=TRUE, n.chains=3, n.iter=2000, n.burnin=1000, n.thin=5)
 
}