SCR_mod3.knit

Introduction

Seeing how model 1 performed so much better than model 2 strongly suggests that combining data from marked and unmarked individuals should be very helpful, at least in estimating the global parameters \(\lambda_0\) and \(\sigma^2\). This will require that we have two likelihoods, one for the marked and one for the unmarked individuals.

Recall that we are still assuming that we know the true number of individuals. Furthermore, incorporating marked individuals will probably not solve the identifiability issue associated with estimating the activity centers for the unmarked individuals.

Simulating data for model 3

Here is how we simulate data for this model.

rm(list=ls())
library('ggplot2')
set.seed(1)

sig2=5
lambda0=10
xmax=ymax=100
nind=20
nind.marked=5
nind.unmark=nind-nind.marked
coord.acenter=data.frame(x=runif(nind,min=25,max=75),
                         y=runif(nind,min=25,max=75))
plot(y~x,data=coord.acenter,xlim=c(0,xmax),ylim=c(0,ymax))

ntraps=25 #number of traps
tmp=seq(from=1,to=xmax,length.out=5)
coord.traps=expand.grid(x=tmp,
                        y=tmp)
points(y~x,data=coord.traps,col='red',pch=19,cex=0.5)

ntimes=10 #encounter times

#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 data
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]
  }
}

#separate marked data
ind.marked=1:nind.marked
ind.unmarked=(nind.marked+1):nind
zarray.marked=zarray[ind.marked,,]
zarray.unmarked=zarray[ind.unmarked,,]

#calculate summary for marked and unmarked
nvec.unmarked=apply(zarray.unmarked,2,sum)
zmat.marked=apply(zarray.marked,c(1,2),sum)

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')

#get data
nind=nind.marked+nind.unmark

#data for jags
dat1=list(zmat=zmat.marked,coord.trap=coord.traps,
          nind.marked=nind.marked,ntimes=ntimes,ntraps=ntraps,
          nvec=nvec.unmarked,nind=nind)

#run JAGS
setwd('U:\\uf\\courses\\bayesian course\\group activities\\17 example scr\\mod3')
model3=jags(model.file="marked unmarked JAGS.R",
            parameters.to.save=params,data=dat1,
            n.chains=nc,n.burnin=nb,n.iter=ni,n.thin=nt,DIC=TRUE)

## 
## Processing function input....... 
## 
## Converting data frame 'coord.trap' to matrix.
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 150
##    Unobserved stochastic nodes: 42
##    Total graph size: 3699
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 5000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 5000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

summary(model3) #the Rhat results suggest that the algorithm has converged

## Summary for model 'marked unmarked JAGS.R'
## Saved parameters: lambda0 sig2 deviance 
## MCMC ran for 0.904 minutes at time 2022-06-15 14:16:47.
## 
## 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
## 
## Successful convergence based on Rhat values (all < 1.1). 
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 24.8 and DIC = 182.329

#look at model parameters
plot(model3$sims.list$lambda0,type='l',ylab='lambda0',xlab='iterations')
abline(h=10,col='red')

plot(model3$sims.list$sig2,type='l',ylab='sigma2',xlab='iterations')
abline(h=5,col='red')

Comments?

Send me an email at