`SCR.fn` <- function(){ ### This function simulates spatial capture-recapture data and fits the model # described in Royle and Gardner (Camera trapping book chapter) # which is a binomial observation model (not multinomial) # as described in Royle et al. (J. Appl. Ecology) library(spdep) library(R2WinBUGS) #Set sigma, lam0, and beta parameters sigma=3 beta=1 lam0=.6 J=10 # number of capture periods #Create a 10x10 grid of traps x=4:13 y=4:13 locs=expand.grid(x,y) ntraps=dim(locs)[1] #Set larger area boundaries (0,0) to (18,18) for "buffer" area llx=0 upx=18 plot(c(llx,upx), c(llx, upx), typ='n') points(locs) #Simulate N home range centers N=28 bearx=runif(N, 0,upx) beary=runif(N, 0,upx) points(bearx, beary, col='red') #Make the Dij matrix of distances from home center to trap locations D=matrix(0, nrow=ntraps, 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=ntraps, ncol=N) k=exp(-D^2/sigma) #to make the encounter history for each simulated bear #first get at probability of encounter Eo=matrix(NA, nrow=ntraps, ncol=N) pmean=matrix(NA, nrow=ntraps, ncol=N) for(i in 1:ntraps){ for(j in 1:N){ Eo[i,j] = lam0*k[i,j] pmean[i,j] = 1 -exp(-Eo[i,j]) # Pr(capture) = 1-exp(-lam0*k[i,j]) , which differs from Royle et al. (J. Appl. Ecol.) ### k[i,j] = exp(-d[i,j]^2 / sigma) } } #Now, make enounter history for each bear over T time steps nij=matrix(0, nrow=ntraps, ncol=N) for(i in 1:ntraps){ for(j in 1:N){ nij[i,j]=rbinom(1, J, pmean[i,j]) } } nij=t(nij) ##########Set up the variables for WinBUGS Madd=100 M=Madd+N aug=matrix(0, nrow=Madd, ncol=ntraps) Yarray=as.matrix(rbind(nij, aug)) trapmat=as.matrix(locs) xl=llx yl=llx yu=upx xu=upx data1<-list(y=Yarray, M=M, trapmat=trapmat, J=J, ntraps=ntraps, xu=xu, xl=xl, yu=yu, yl=yl) sink("SCRmodel.txt") cat(" model { sigma~dunif(0, 15) lam0~dgamma(.1,.1) psi~dunif(0, 1) for (i in 1:M){ z[i]~dbern(psi) SX[i]~dunif(xl, xu) SY[i]~dunif(yl, yu) for(j in 1:ntraps) { D2[i,j] <- pow(SX[i]-trapmat[j,1], 2) + pow(SY[i]-trapmat[j,2],2) Eo[i,j] <- lam0*exp(-D2[i,j]/sigma) pmean[i,j]<-1-(exp(-Eo[i,j])) tmp[i,j]<-pmean[i,j]*z[i] y[i,j]~dbin(tmp[i,j],J) } } N<-sum(z[1:M]) } ",fill=TRUE) sink() params = list('lam0', 'sigma', 'N', 'psi') sigma=5 psi=.6 SX=as.vector(runif(M,2, 10)) SY=as.vector(runif(M,2, 10)) z=as.vector(rbinom(M, 1, .5)) inits = function() {list(z=z, psi=psi, lam0=.2, SX=SX, SY=SY, sigma=sigma) } fit = bugs(data1, inits, params, model.file="SCRmodel.txt", debug=FALSE, n.chains=3, n.iter=1200, n.burnin=1000, n.thin=3) fit }