Introduction
Recall that, in this model, we assume that we only have unmarked individuals. In other words, we only know the total number of individuals seen in each trap and each occasion. For the moment being, we will assume that we know the overall number of individuals that are present. The goal is to show that this model struggles to converge even in this very optimistic scenario because of identifiability issues regarding estimating the individual parameters \(\lambda_0\) and \(\sigma^2\) as well as the center of each activity center.
Before we jump in, recall that:
\[n_j =\sum_{i=1}^n \sum_{k=1}^k z_{ijk} \sim Poisson(K \times \sum_{i=1}^n \lambda_{ij}) \]
where \(n_j\) is the number of individuals in trap \(j\). This is going to be the main expression in our SCR model.
Simulating data for model 2
Here is how we simulate data for this model.
rm(list=ls())
# library('proxy')
library('ggplot2')
set.seed(1)
#parameter values
sig2=20
lambda0=8
#number of unmarked individuals
nind=10
#spatial extent
xmax=ymax=100
#coordinates of activity center
coord.acenter=data.frame(x=runif(nind,min=0,max=xmax),
y=runif(nind,min=0,max=ymax))
#traps and their coordinates
ntraps=25 #number of traps
tmp=seq(from=1,to=xmax,length.out=5)
coord.traps=expand.grid(x=tmp,
y=tmp)
#encounter times
ntimes=5
#calculate distances
dist1=matrix(NA,ntraps,nind)
for (j in 1:ntraps){
for (i in 1:nind){
x2=(coord.traps$x[j]-coord.acenter$x[i])^2
y2=(coord.traps$y[j]-coord.acenter$y[i])^2
dist1[j,i]=sqrt(x2+y2)
}
}
# dist1a=proxy::dist(x=coord.traps,y=coord.acenter)
# unique(dist1-dist1a)
#generate z
n=matrix(NA,ntraps,ntimes)
zarray=array(NA,dim=c(nind,ntraps,ntimes))
for (t1 in 1:ntimes){
z=matrix(NA,ntraps,nind)
for (i in 1:nind){
g=exp(-dist1[,i]/sig2)
z[,i]=rpois(ntraps,lambda0*g)
zarray[i,,t1]=z[,i]
}
n[,t1]=apply(z,1,sum)
}
#sum over time steps
N=apply(n,1,sum)
coord.traps$N=N
ggplot()+
geom_tile(data=coord.traps,aes(x=x,y=y,fill=N))+
geom_point(data=coord.acenter,aes(x=x,y=y),col='red')
Estimating parameters
To estimate the model parameters, we rely on the code below
library('jagsUI')
# MCMC settings
ni <- 10000 #number of iterations
nt <- 5 #interval to thin
nb <- 5000 #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;
nind=10
nvec=coord.traps$N
dat1=list(nvec=nvec,coord.trap=coord.traps,
nind=nind,ntimes=ntimes,ntraps=ntraps)
#run JAGS
setwd('U:\\uf\\courses\\bayesian course\\group activities\\17 example scr\\mod2')
model2=jags(model.file="unmarked JAGS.R",
parameters.to.save=params,data=dat1,
n.chains=nc,n.burnin=nb,n.iter=ni,n.thin=nt,DIC=TRUE,
verbose=F)
summary(model2) #the Rhat results suggest that the algorithm has converged## Summary for model 'unmarked JAGS.R'
## Saved parameters: lambda0 sig2 coord.acenter deviance
## MCMC ran for 0.429 minutes at time 2022-06-15 14:16:08.
##
## For each of 3 chains:
## Adaptation: 100 iterations (sufficient)
## Burn-in: 5000 iterations
## Thin rate: 5 iterations
## Total chain length: 10100 iterations
## Posterior sample size: 1000 draws
##
## **WARNING** Rhat values indicate convergence failure.
##
## DIC info: (pD = var(deviance)/2)
## pD = 22.2 and DIC = 189.143# model2
#look at model parameters
plot(model2$sims.list$lambda0,type='l',ylab='lambda0',xlab='iterations')
abline(h=8,col='red')
plot(model2$sims.list$sig2,type='l',ylab='sigma2',xlab='iterations')
abline(h=20,col='red')
Notice that we have increased the number of iterations relative to what we used for model 1 and we are still having a hard time making the model converge even though we are assuming, rather optimistically, that the overall number of individuals is known. This is partly because we have a very strong correlation between \(\lambda_0\) and \(\sigma^2\).
#look at correlation between these parameters
plot(model2$sims.list$lambda0,model2$sims.list$sig2,
ylab='sigma2',xlab='lambda0')
Furthermore, if you look at the Rhat and n.eff values for the coordinates of the activity centers, you will see a lot of issues there as well. This makes sense. You can probably tell that there are more activity centers in certain areas but not in others just based on which traps had more individuals. However, because none of the individuals are marked, you cannot tell to which individuals these activity centers belong to.
#look at Rhat and n.eff
model2## JAGS output for model 'unmarked JAGS.R', generated by jagsUI.
## Estimates based on 3 chains of 10000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 5000 iterations and thin rate = 5,
## yielding 3000 total samples from the joint posterior.
## MCMC ran for 0.429 minutes at time 2022-06-15 14:16:08.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat
## lambda0 7.743 1.510 5.221 7.629 11.020 FALSE 1 1.024
## sig2 20.909 2.728 16.317 20.654 26.790 FALSE 1 1.029
## coord.acenter[1,1] 57.556 31.574 1.920 67.272 98.220 FALSE 1 1.308
## coord.acenter[2,1] 69.322 24.775 6.928 75.092 98.621 FALSE 1 1.097
## coord.acenter[3,1] 64.885 28.217 5.012 72.850 98.730 FALSE 1 1.100
## coord.acenter[4,1] 50.910 32.846 1.584 60.160 97.329 FALSE 1 1.026
## coord.acenter[5,1] 44.914 32.799 1.141 36.529 96.464 FALSE 1 1.010
## coord.acenter[6,1] 53.602 32.130 2.248 62.562 98.038 FALSE 1 1.072
## coord.acenter[7,1] 51.974 31.770 1.997 60.832 97.778 FALSE 1 1.144
## coord.acenter[8,1] 62.220 29.710 2.603 71.968 98.406 FALSE 1 1.084
## coord.acenter[9,1] 61.892 29.214 4.980 71.272 98.030 FALSE 1 1.027
## coord.acenter[10,1] 61.449 31.002 2.025 71.145 98.415 FALSE 1 1.101
## coord.acenter[1,2] 54.983 25.205 10.051 56.981 96.284 FALSE 1 1.106
## coord.acenter[2,2] 54.713 25.999 13.068 50.765 98.055 FALSE 1 1.022
## coord.acenter[3,2] 51.503 25.187 10.209 47.551 96.562 FALSE 1 1.009
## coord.acenter[4,2] 57.809 25.635 10.352 65.285 97.329 FALSE 1 1.099
## coord.acenter[5,2] 55.804 27.076 6.873 65.516 96.095 FALSE 1 1.011
## coord.acenter[6,2] 54.256 26.622 9.137 58.895 96.674 FALSE 1 1.019
## coord.acenter[7,2] 51.214 27.962 7.148 51.908 96.580 FALSE 1 1.068
## coord.acenter[8,2] 55.854 25.466 10.404 58.317 97.244 FALSE 1 1.005
## coord.acenter[9,2] 52.367 26.919 8.049 49.223 97.139 FALSE 1 1.051
## coord.acenter[10,2] 58.408 24.661 18.078 62.411 98.122 FALSE 1 1.026
## deviance 166.931 6.716 155.554 166.141 182.070 FALSE 1 1.017
## n.eff
## lambda0 146
## sig2 109
## coord.acenter[1,1] 11
## coord.acenter[2,1] 43
## coord.acenter[3,1] 33
## coord.acenter[4,1] 86
## coord.acenter[5,1] 223
## coord.acenter[6,1] 33
## coord.acenter[7,1] 18
## coord.acenter[8,1] 31
## coord.acenter[9,1] 92
## coord.acenter[10,1] 29
## coord.acenter[1,2] 24
## coord.acenter[2,2] 94
## coord.acenter[3,2] 234
## coord.acenter[4,2] 27
## coord.acenter[5,2] 225
## coord.acenter[6,2] 112
## coord.acenter[7,2] 35
## coord.acenter[8,2] 374
## coord.acenter[9,2] 44
## coord.acenter[10,2] 81
## deviance 128
##
## **WARNING** Rhat values indicate convergence failure.
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 22.2 and DIC = 189.143
## DIC is an estimate of expected predictive error (lower is better).
Comments?
Send me an email at