Introduction
Recall that, in this model, we assume that we only have marked individuals because the goal is to introduce the main parameters that govern the SCR model. More specifically, we want to estimate \(\lambda_0\), \(\sigma^2\), and the spatial coordinates of each activity center. To do this, we will rely on the data within zarray and the coordinates of the traps, which are known.
One statistical fact that is important to know is that if:
\[x_i \sim Poisson(\lambda_i)\]
then
\[\sum_i x_i \sim Poisson(\sum_i \lambda_i)\]
As a result of the statistical fact described above, we have that:
\[\sum_k z_{ijk} \sim Poisson(K \times \lambda_{ij}) \]
where K is the total number of occasions. Because of this, instead of working with zarray, we can sum over all occasions and work just with a regular matrix. We start by first importing the data and summing over all occasions.
rm(list=ls())
set.seed(1)
#import data
setwd('U:\\uf\\courses\\bayesian course\\group activities\\17 example scr\\mod1')
coord.traps=data.matrix(read.csv('coord traps.csv'))
ntraps=nrow(coord.traps)
tmp=read.csv('zarray.csv')
N=3
ntimes=5
zarray=array(tmp$V1,dim=c(N,ntraps,ntimes))
#sum across all occasions to generate summarized matrix
zmat=apply(zarray,c(1,2),sum)
rownames(zmat)=paste0('id',1:N)
colnames(zmat)=paste0('trap.id',1:ntraps)
zmat## trap.id1 trap.id2 trap.id3 trap.id4 trap.id5 trap.id6 trap.id7 trap.id8
## id1 0 0 0 0 0 0 0 0
## id2 0 5 4 1 0 0 8 13
## id3 0 0 0 0 0 0 0 0
## trap.id9 trap.id10 trap.id11 trap.id12 trap.id13 trap.id14 trap.id15
## id1 0 0 0 3 1 0 0
## id2 1 1 0 2 1 0 0
## id3 0 0 0 2 1 0 0
## trap.id16 trap.id17 trap.id18 trap.id19 trap.id20 trap.id21 trap.id22
## id1 1 7 4 0 0 4 15
## id2 0 0 0 0 0 0 0
## id3 0 1 3 7 0 0 1
## trap.id23 trap.id24 trap.id25
## id1 6 0 0
## id2 0 0 0
## id3 10 8 1Next, we fit the model to then determine if it was able to estimate all the parameters well.
library('jagsUI')
# MCMC settings
ni <- 1000 #number of iterations
nt <- 5 #interval to thin
nb <- 500 #number of iterations to discard as burn-in
nc <- 3 #number of chains
#parameters to monitor
params=c('lambda0','sig2','coord.acenter')
#basic settings
ntraps=nrow(coord.traps); #ntraps;
N=nrow(zmat)
dat1=list(zmat=zmat,coord.trap=coord.traps,
nind=N,ntimes=ntimes,ntraps=ntraps)
#create function that sets the initial values
inits=function(){
list(lambda0=rgamma(1,1,1), sig2=rgamma(1,1,1),
coord.acenter=cbind(runif(N,min=min(coord.traps[,'x']),max=max(coord.traps[,'x'])),
runif(N,min=min(coord.traps[,'y']),max=max(coord.traps[,'y']))))
}
#run JAGS
setwd('U:\\uf\\courses\\bayesian course\\group activities\\17 example scr\\mod1')
model1=jags(model.file="marked JAGS.R",inits=inits,
parameters.to.save=params,data=dat1,
n.chains=nc,n.burnin=nb,n.iter=ni,n.thin=nt,DIC=TRUE,
verbose=F)
summary(model1) #the Rhat results suggest that the algorithm has converged## Summary for model 'marked JAGS.R'
## Saved parameters: lambda0 sig2 coord.acenter deviance
## MCMC ran for 0.012 minutes at time 2023-03-17 13:56:16.
##
## For each of 3 chains:
## Adaptation: 100 iterations (sufficient)
## Burn-in: 500 iterations
## Thin rate: 5 iterations
## Total chain length: 1100 iterations
## Posterior sample size: 100 draws
##
## Successful convergence based on Rhat values (all < 1.1).
##
## DIC info: (pD = var(deviance)/2)
## pD = 7.5 and DIC = 122.972#look at model parameters
plot(model1$sims.list$lambda0,type='l',ylab='lambda0',xlab='Iterations')
abline(h=8,col='red')
plot(model1$sims.list$sig2,type='l',ylab='sigma2',xlab='Iterations')
abline(h=10,col='red')
These results show that it is straight-forward to estimate \(\lambda_0\) and \(\sigma^2\) if you have marked individuals.
Comments?
Send me an email at