Multipopulation Mark-Recapture
Matt Falcy
February 8, 2019
The Situation
There are 83 stream units (hereafter “population units”) where two-sample mark-recapture studies were performed. Juvenile steelhead were captured for marking and recaptures by blockneting a population unit and electrofishing. These 83 studies span 3 years. Each population unit occurs in either a Pool or Fast water. All population units in Pools also have a snorkel count. Two covaraites, Cover and Clarity, are recorded at each population unit.
Questions
- What proportion of the total abundance is observed by snorkelers?
- Does detection probability (from mark-recapture, not snorkel) differ in pool vs. fast?
- Does the ratio of abundnces in Pools vs. Fast change across years?
Approach
I’m going to fit a simple mark-recapture model called “M0” (sensu Otis, 1982). Model “M0” assumes that detection probablity on the first and second sampling occasions are identical. It also assumes no behavioral response (“trap happy”) and no heterogeneity (all individuals have same capture probability).
Several aspects extend beyond the simple “M0”. First, I will use a data augmentation approach (Kery and Schaub 2012). Second, I wil fit all 83 populations simultaneously. Third, I evaluate question 1 by including the snorkel count as a binomial proportion of the total abudnance in each population unit.
Load the data and clean it up
setwd("//Kalawatseti/home/falcym/docs/Projects/People/Constable/Detection")
D<-read.csv('Sth.csv')
D<-D[D$SthdMarked>0,] # Remove all sites where no fish are marked.
D[,7]>D[,5]# Verify that there aren't more fish recaped with marks than initial marks
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
D<-D[-16,] #Remove "TRUE"
Y2016<-which(D$Year==2016)
Y2017<-which(D$Year==2017)
Y2018<-which(D$Year==2018 | D$Year==1900)#1900 is actually 2018 as per R. Constable email 2/7/2019.
pool<-which(D$UnitType=='P')
fast<-which(D$UnitType=='F')
UnitType<-vector()
UnitType[pool]<-1
UnitType[fast]<-2
D.pool<-D[pool,]
D.fast<-D[fast,]
D[1:10,]
## Year Site Unit UnitType SthdMarked SthdCap SthdRecap
## 1 1900 9000018 1 P 12 1 10
## 2 1900 9000018 2 F 3 2 2
## 3 1900 9000018 3 P 19 28 10
## 4 1900 9000018 4 F 21 5 17
## 5 1900 9000018 5 P 14 17 8
## 6 1900 9000018 6 P 14 5 7
## 7 1900 14768018 1 P 1 3 0
## 8 1900 14768018 2 F 5 1 3
## 9 1900 14768018 3 P 9 4 6
## 10 1900 14768018 4 F 5 5 2
## SnorkelCountSthd Cover Clarity
## 1 7 1 3
## 2 0 0 2
## 3 29 2 2
## 4 0 2 2
## 5 26 1 2
## 6 15 1 3
## 7 1 1 3
## 8 0 0 1
## 9 6 2 2
## 10 0 1 2
Make Capture Histories
I want to make capture histories for each individual from the data contained in “D”. Here are the relations: [1 1]. Captured on both occasions. This is “Recap” [1 0]. Captured on the first occasion but not the second. This is “Marked” - “Recap” [0 1]. Not captured on the first occasion but captured on the second. This is “Cap” [0 0]. Never seen. This is what we want to know. Using data augmentation, I’ll created a bunch of these histories then estimate the fraction of these fake individuals that must be real. Because I will assume constant detection probablity across events, the histories [1 0] and [0 1] will collapse into “1” for observed once over two occasions.
While im in the population loop below, I’ll also process the covaraites and the snorkel count. The output, “y” contains two columns. The first column is number of times the individual was seen. The second column is the population unit ID. “”
###
h.11<-vector()
h.10<-vector()
h.01<-vector()
Y.11<-data.frame()
Y.10<-data.frame()
Y.01<-data.frame()
Y.00<-data.frame()
Y<-array()
y<-array()
C<-vector()
Cov1.p<-vector()
Cov2.p<-vector()
Cov.N<-vector()
for (p in 1:dim(D)[1]){
h.11[p]<-D$SthdRecap[p]
h.10[p]<-D$SthdMarked[p]-D$SthdRecap[p]
h.01[p]<-D$SthdCap[p]
Y.11<-matrix(rep(c(1,1),each=h.11[p]),c(h.11[p],2))
Y.10<-matrix(rep(c(1,0),each=h.10[p]),nrow=h.10[p])
Y.01<-matrix(rep(c(0,1),each=h.01[p]),nrow=h.01[p])
Y.00<-matrix(rep(c(0,0),each=100),nrow=100)
Y<-rbind(if(length(Y.11)>0){Y.11}else{},if(length(Y.10)>0){Y.10}else{},if(length(Y.01)>0){Y.01}else{},if(length(Y.00)>0){Y.00}else{})
y<-rbind(y,cbind(apply(Y,1,sum),rep(p,h.11[p]+h.10[p]+h.01[p]+100)))
C[p]<-h.11[p]+h.10[p]+h.01[p]+100# sum of real + augmented indiduals.
Cov1.p[p]<-(D$Clarity[p]-2) #Covariate 1
Cov2.p[p]<-D$Cover[p] #Covariate 2
Cov.N[p]<-D$SnorkelCountSthd[p]
}
Cov.N[fast]<-NA #zeros in fast are actually NA
y<-y[!is.na(y[,1]),]
y[1:10,]
## [,1] [,2]
## 2 1
## 2 1
## 2 1
## 2 1
## 2 1
## 2 1
## 2 1
## 2 1
## 2 1
## 2 1
JAGS Model
JAGS function that answers the questions above with these data. Variables “densisity1” and “density2”, can be uncommented if WAIC is desired.
M0.pop<-function(){
N[1]<-sum(z[1:C[1]])#Derived
cs[1]<-C[1]#cumulative sum
for (j in 2:pop){
cs[j]<-cs[j-1]+C[j]#cumulative sum
N[j]<-sum(z[(cs[j-1]+1):cs[j]])#Derived
}
b0~dnorm(0,0.001)#Prior
b1~dnorm(0,0.001)#Prior
b2~dnorm(0,0.001)#Prior
pN~dbeta(1,1)
for (j in 1:pop){
psi[j]~dunif(0,1)#Prior
logit(p[j])<-b0+b1*Cov1.p[j]+b2*Cov2.p[j]
Cov.N[j]~dbin(pN,N[j])#Snorkle count~ dbin(p, N)
#density2[j]<-dbin(Cov.N[j],pN,N[j])
}
for (i in 1:total){
z[i]~dbin(psi[y[i,2]],1)
mu[i]<-z[i]*p[y[i,2]]
y[i,1]~dbin(mu[i],2)
#density1[i]<-dbin(y[i,1],mu[i],2)
}
mu.pool<-sum(N[pool])/npool
mu.fast<-sum(N[fast])/nfast
N.tot16<-sum(N[Y2016])
N.tot17<-sum(N[Y2017])
N.tot18<-sum(N[Y2018])
#Estimate Abundance with snorkel
for (i in 1:10){
Ni[i]<-(10*i)/pN
}
}#close function
Run JAGS
data<-list("y"=y,"C"=C,'total'=dim(y)[1],'pop'=max(y[,2]),"Cov.N"=Cov.N,"Cov1.p"=Cov1.p,"Cov2.p"=Cov2.p,"Y2016"=Y2016,"Y2017"=Y2017,"Y2018"=Y2018,"pool"=pool,"fast"=fast,"npool"=length(pool),"nfast"=length(fast))
Nchains<-3#Number of Markov chain Monte Carlo simulations
inits<-function(){
list("b0"=rnorm(1,0,0.1),"b1"=rnorm(1,0,0.1),"b2"=rnorm(1,0,0.1),"psi"=runif(max(y[,2])),"z"=rep(1,dim(y)[1]))
}
params=c("p","N","b0","b1","b2","psi","pN","Ni","N.tot16","N.tot17","N.tot18","mu.pool","mu.fast")
library(R2jags)
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
##
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
##
## traceplot
jags.out<-jags.parallel(model.file=M0.pop,data=data,inits=inits,parameters.to.save=params,n.chains=3,n.iter=10000,n.burnin=2000,n.thin=5)
Results
(jags.out$BUGSoutput$mean)
## $N
## [1] 15.475417 7.211042 61.311667 34.408750 45.140208 24.404167 4.575417
## [8] 10.552083 16.252292 13.758542 6.651042 3.884792 3.100625 3.745208
## [15] 1.382917 4.537500 4.509583 6.323958 3.804167 7.070208 36.715625
## [22] 9.993958 27.271875 78.676458 12.261667 3.036042 4.119167 1.753958
## [29] 10.766667 7.257917 8.546667 3.900625 1.194792 4.128333 5.161458
## [36] 1.442500 2.460417 2.636250 3.801458 2.365417 2.691667 1.378542
## [43] 8.452083 1.366667 6.177292 3.008125 4.927083 1.618958 3.500625
## [50] 7.163958 3.507292 3.749375 5.181667 4.512500 3.677500 2.596042
## [57] 3.514583 1.739375 6.050000 1.436458 12.969583 9.995417 4.708333
## [64] 5.138958 5.025000 2.419167 8.567500 7.169583 11.360000 11.773542
## [71] 7.822917 5.033333 2.665625 6.086667 8.602917 1.366875 1.421250
## [78] 2.394583 5.538750 3.878958 6.058542 3.880208 2.669375
##
## $N.tot16
## [1] 99.82063
##
## $N.tot17
## [1] 122.9148
##
## $N.tot18
## [1] 495.6523
##
## $Ni
## [1] 20.8579 41.7158 62.5737 83.4316 104.2895 125.1474 146.0053
## [8] 166.8632 187.7211 208.5790
##
## $b0
## [1] -0.07291058
##
## $b1
## [1] 0.3902466
##
## $b2
## [1] 0.06824644
##
## $deviance
## [1] 1679.027
##
## $mu.fast
## [1] 8.065712
##
## $mu.pool
## [1] 9.006743
##
## $p
## [1] 0.5947356 0.4819197 0.5158037 0.5158037 0.4988380 0.5947356 0.5947356
## [8] 0.3886941 0.5158037 0.4988380 0.4988380 0.5781546 0.4988380 0.5947356
## [15] 0.5947356 0.6105324 0.4819197 0.4988380 0.5947356 0.4988380 0.5781546
## [22] 0.5781546 0.5947356 0.4819197 0.4988380 0.4988380 0.4988380 0.4819197
## [29] 0.5158037 0.4819197 0.5158037 0.5781546 0.5947356 0.4988380 0.5781546
## [36] 0.5781546 0.5781546 0.5781546 0.5781546 0.6252975 0.6105324 0.6105324
## [43] 0.5947356 0.6105324 0.5158037 0.4988380 0.4988380 0.5158037 0.5947356
## [50] 0.4988380 0.5947356 0.6105324 0.5781546 0.4819197 0.4819197 0.5947356
## [57] 0.5947356 0.4819197 0.5947356 0.5781546 0.6105324 0.4819197 0.6105324
## [64] 0.5158037 0.5324738 0.5947356 0.5781546 0.4988380 0.4819197 0.4988380
## [71] 0.4988380 0.5324738 0.5781546 0.5947356 0.5947356 0.5947356 0.5781546
## [78] 0.6105324 0.5158037 0.5781546 0.5947356 0.5781546 0.4988380
##
## $pN
## [1] 0.481211
##
## $psi
## [1] 0.14364293 0.07648651 0.41838481 0.27678708 0.34727395 0.20967959
## [7] 0.05258222 0.10643040 0.15068392 0.13195291 0.07179411 0.04631217
## [13] 0.03930981 0.04552258 0.02306227 0.05245849 0.05262118 0.06898257
## [19] 0.04598341 0.07425949 0.28637683 0.09935499 0.22635903 0.50112686
## [25] 0.11997538 0.03944763 0.04837943 0.02745244 0.10793250 0.07805918
## [31] 0.08866230 0.04704694 0.02130965 0.04847427 0.05824117 0.02362251
## [37] 0.03350116 0.03508260 0.04504834 0.03231284 0.03543987 0.02310230
## [43] 0.08751349 0.02276168 0.06836657 0.03895404 0.05556979 0.02591273
## [49] 0.04280126 0.07650374 0.04268798 0.04553175 0.05819253 0.05295163
## [55] 0.04481039 0.03463697 0.04338904 0.02630913 0.06559116 0.02357345
## [61] 0.12364184 0.10110641 0.05410764 0.05837357 0.05640457 0.03244744
## [67] 0.08750081 0.07643114 0.11190139 0.11525467 0.08357987 0.05687035
## [73] 0.03483656 0.06636731 0.08780293 0.02291833 0.02325571 0.03241836
## [79] 0.06221212 0.04666682 0.06643450 0.04629686 0.03565050
densplot(as.mcmc(jags.out$BUGSoutput$sims.matrix[,1]))
HPDinterval(as.mcmc(jags.out$BUGSoutput$sims.list$N[,1]),0.95)#
## lower upper
## var1 13 18
## attr(,"Probability")
## [1] 0.95
densplot(as.mcmc(jags.out$BUGSoutput$sims.list$pN),show.obs=FALSE,xlab='Proportion Observed by Snorkelers',ylab='Probability Density',main='Steelhead')
par(mfrow=c(2,1))
densplot(as.mcmc(jags.out$BUGSoutput$sims.list$b1),show.obs=FALSE,xlab='Clarity',ylab='Probability Density',main='Steelhead')
densplot(as.mcmc(jags.out$BUGSoutput$sims.list$b2),show.obs=FALSE,xlab='Cover',ylab='Probability Density')
mu.fast<-vector()
for (i in 1:(length(fast))){
k<-fast[i]
mu.fast[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$p[,k]))#
}
mu.pool<-vector()
for (i in 1:(length(pool))){
k<-pool[i]
mu.pool[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$p[,k]))#
}
library(vioplot)
## Loading required package: sm
## Package 'sm', version 2.2-5.6: type help(sm) for summary information
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
vioplot(mu.pool,mu.fast, names=c("Fast","Pool"),col="green")
title("Steelhead Detection Probability")
F2016<-intersect(Y2016,fast)
P2016<-intersect(Y2016,pool)
F2017<-intersect(Y2017,fast)
P2017<-intersect(Y2017,pool)
F2018<-intersect(Y2018,fast)
P2018<-intersect(Y2018,pool)
mu.F2016<-vector()
for (i in 1:(length(F2016))){
k<-F2016[i]
mu.F2016[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$N[,k]))#
}
mu.P2016<-vector()
for (i in 1:(length(P2016))){
k<-P2016[i]
mu.P2016[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$N[,k]))#
}
mu.F2017<-vector()
for (i in 1:(length(F2017))){
k<-F2017[i]
mu.F2017[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$N[,k]))#
}
mu.P2017<-vector()
for (i in 1:(length(P2017))){
k<-P2017[i]
mu.P2017[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$N[,k]))#
}
mu.F2018<-vector()
for (i in 1:(length(F2018))){
k<-F2018[i]
mu.F2018[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$N[,k]))#
}
mu.P2018<-vector()
for (i in 1:(length(P2018))){
k<-P2018[i]
mu.P2018[i]<-mean(as.mcmc(jags.out$BUGSoutput$sims.list$N[,k]))#
}
vioplot(mu.F2016,mu.P2016,mu.F2017,mu.P2017,mu.F2018,mu.P2018,names=c("Fast 2016","Pool 2016","Fast 2017","Pool 2017","Fast 2018","Pool 2018"))
title("Steelhead Abundance")