#' Simulate data from camera trap SCR study with partial identity and simulate a random trap failure process
#' @param N a vector indicating the number of individuals to simulate
#' @param lam01 the single side detection function hazard rate
#' @param lam02 the double side detection function hazard rate
#' @param sigma the spatial scale parameter
#' @param K the number of capture occasions
#' @param X the K x 3 matrix of trap locations and number of cameras at each trap. columns are X, Y, #cams
#' @param buff the distance to buffer the trapping array in the X and Y dimensions to produce the state space
#' @param failprob the probability a camera or set of cameras at a site will fail on each occasion
#' @param faildur the number of occasions a camera remains inoperable
#' @param failtype 1 if cameras at a site fail together and 2 if they fail independently
#' @param obstype observation type, either "bernoulli" or "poisson"
#' @return a list containing the both, left, and right data sets, the both plus left only and both plus right
#' only data sets (stored in y.obs), the activity centers, the trap object, and several other data objects and summaries.
#' @description See the sim2side help file for a description of the capture process. The camera failure process operates
#' like this: cameras fail with probabilty failprob and remain disabled for faildur occasions. If failtype=1, traps at
#' single sites fail together (bear attack?) and if failtype=2, they fail independently. This function was written to investigate
#' the importance of using the 2D trap operation file.
#' @author Ben Augustine
#' @export
sim2sidetf <-
function(N=120,lam01=0.1,lam02=0.2,sigma=0.50,K=10,X=X,buff=3,failprob=0.05,faildur=2,failtype=2,obstype="bernoulli"){
library(abind)
#######Capture process######################
# # simulate a population of activity centers
s<- cbind(runif(N, min(X[,1])-buff,max(X[,1])+buff), runif(N,min(X[,2])-buff,max(X[,2])+buff))
D<- e2dist(s,X)
lamd<- abind(lam01*exp(-D*D/(2*sigma*sigma)),lam02*exp(-D*D/(2*sigma*sigma)),along=3)
J<- nrow(X)
if(failtype==1){#If traps fail together
#Simulate trap failure/correction process
onoff=matrix(NA,nrow=J,ncol=K)
onoff[,1]=1
offoccs=rep(0,J)
on=rep(1,J)
for(k in 2:K){
offoccs[offoccs>0]=offoccs[offoccs>0]+1 #if off, add another occasion
onoff[,k]=rbinom(J,on,1-failprob)
offoccs[(onoff[,k]==0)&(on==1)]=1 #Turn off stations that just failed
on[(onoff[,k]==0)&(on==1)]=0 #Turn off stations that just failed
backon=which(offoccs==faildur)
if(length(backon)>0){
on[backon]=1
offoccs[backon]=0
}
}
# Initialize the encounter history data for left and right sides
left <-right <- both <-array(0,dim=c(N,J,K))
if(obstype=="bernoulli"){
pd=1-exp(-lamd)
for(i in 1:N){
for(j in 1:J){
for(k in 1:K){
#P(A or B)=P(A)+P(B)-P(A and B)
left[i,j,k]=rbinom(1,1,((X[j,3]==1)*pd[i,j,1]+(X[j,3]==2)*(2*pd[i,j,1]-pd[i,j,1]^2))*onoff[j,k]) #single side p. two chances for capture with 2 cameras
right[i,j,k]=rbinom(1,1,((X[j,3]==1)*pd[i,j,1]+(X[j,3]==2)*(2*pd[i,j,1]-pd[i,j,1]^2))*onoff[j,k]) #single side p
both[i,j,k]=rbinom(1,1,pd[i,j,2]*onoff[j,k])*(X[j,3]==2)*1 #both side lambda multiplied by indicator for 2 traps at site
}
}
}
}else if(obstype=="poisson"){
for(i in 1:N){
for(j in 1:J){
for(k in 1:K){
left[i,j,k]=rpois(1,X[j,3]*lamd[i,j,1]*onoff[j,k]) #single side lambda multiplied by number of traps at site
right[i,j,k]=rpois(1,X[j,3]*lamd[i,j,1]*onoff[j,k]) #single side lambda multiplied by number of traps at site
both[i,j,k]=rpois(1,lamd[i,j,2]*onoff[j,k])*(X[j,3]==2)*1 #both side lambda multiplied by indicator for 2 traps at site
}
}
}
}else{
stop("observation model not recognized")
}
}else if(failtype==2){#If traps don't fail together
#Simulate trap failure/correction process. If 1 cam corrected, the other is corrected.
onoff=array(0,dim=c(J,K,2))
onoff[,1,1:2]=1
offoccs=matrix(0,nrow=J,ncol=2)
on=matrix(1,nrow=J,ncol=2)
singles=which(X[,3]==1)
doubles=which(X[,3]==2)
J2=length(doubles)
offoccs[singles,2]=0
on[singles,2]=0
onoff[singles,1,2]=0
for(k in 2:K){
offoccs[offoccs>0]=offoccs[offoccs>0]+1 #if off, add another occasion
onoff[,k,1]=rbinom(J,on[,1],1-failprob)
onoff[doubles,k,2]=rbinom(J2,on[doubles,2],1-failprob)
offoccs[(onoff[,k,1]==0)&(on[,1]==1),1]=1 #Turn off stations that just failed
offoccs[(onoff[singles,k,2]==0)&(on[singles,2]==1),2]=1 #Turn off stations that just failed
on[(onoff[,k,1]==0)&(on[,1]==1),1]=0 #Turn off stations that just failed
on[(onoff[singles,k,2]==0)&(on[singles,2]==1),2]=0 #Turn off stations that just failed
#Turn back on together
backon=which(apply(offoccs,1,max)==faildur)
if(length(backon)>0){
on[backon[backon%in%doubles],]=1
on[backon[backon%in%singles],1]=1
offoccs[backon,]=0
}
}
#number of cams on per trap/occasion
onoff2=onoff[,,1]+onoff[,,2]
# Initialize the encounter history data for left and right sides
left <-right <- both <-array(0,dim=c(N,J,K))
if(obstype=="bernoulli"){
pd=1-exp(-lamd)
for(i in 1:N){
for(j in 1:J){
for(k in 1:K){
if(onoff2[j,k]==2){
left[i,j,k]=rbinom(1,1,((X[j,3]==1)*pd[i,j,1]+(X[j,3]==2)*(2*pd[i,j,1]-pd[i,j,1]^2))) #single side p. two chances for capture with 2 cameras
right[i,j,k]=rbinom(1,1,((X[j,3]==1)*pd[i,j,1]+(X[j,3]==2)*(2*pd[i,j,1]-pd[i,j,1]^2))) #single side p
both[i,j,k]=rbinom(1,1,pd[i,j,2])*(X[j,3]==2) #both side lambda multiplied by indicator for 2 traps at site
}else if(onoff2[j,k]==1){
left[i,j,k]=rbinom(1,1,pd[i,j,1])
right[i,j,k]=rbinom(1,1,pd[i,j,1])
}
}
}
}
}else if(obstype=="poisson"){
for(i in 1:N){
for(j in 1:J){
for(k in 1:K){
left[i,j,k]=rpois(1,X[j,3]*lamd[i,j,1]*onoff2[j,k]) #single side lambda multiplied by number of traps at site
right[i,j,k]=rpois(1,X[j,3]*lamd[i,j,1]*onoff2[j,k]) #single side lambda multiplied by number of traps at site
both[i,j,k]=rpois(1,lamd[i,j,2]*onoff2[j,k])*(X[j,3]==2)*1 #both side lambda multiplied by indicator for 2 traps at site
}
}
}
}else{
stop("observation model not recognized")
}
}else{
warning("failtype must be 1 or 2")
}
######Process data#############
IDknown=which(apply(both,1,sum)>0)
Nknown=length(IDknown)
n=sum(apply(both+left+right,1,sum)>0)
nside=c(sum(apply(both,1,sum)>0),sum(apply(left,1,sum)>0),sum(apply(right,1,sum)>0))
#Count true spatial recaps
y2D=apply(both+left+right,c(1,2),sum)
scaps=rowSums(1*(y2D>0))
scaps[scaps>0]=scaps[scaps>0]-1
nscap=sum(scaps>0)
sumscap=sum(scaps)
Srecap=data.frame(nscap=nscap,sumscap=sumscap)
# relabel data sets if needed so that the left side always has more individuals than right after removing known IDs
nl.tmp<- sum(apply(left,1,sum)[apply(both,c(1),sum)==0]>0)
nr.tmp<- sum(apply(right,1,sum)[apply(both,c(1),sum)==0]>0)
if(nr.tmp > nl.tmp){
a<- left
left<- right
right<- a
}
# Move known individuals(after photo ID) to beginning
if(Nknown>0&Nknown<N){
IDunknown=which(is.na(match(1:N,IDknown)))
newleft<- abind(left[IDknown,,],left[IDunknown,,],along=1)
newright<- abind(right[IDknown,,],right[IDunknown,,],along=1)
newboth<- abind(both[IDknown,,],both[IDunknown,,],along=1)
s=s[c(IDknown,IDunknown),]
left<- newleft
right<- newright
both<- newboth
#Recalculate after sorting
IDknown=1:length(IDknown)
# Sort by number of captures in both data set.
#These are now defined to be in order by true IDs 1:N
tb<- apply(both,c(1,3),sum)
cap.both<- apply(tb,1,sum)[(Nknown+1):N]
oo1<- c(1:Nknown, Nknown + rev(order(cap.both)) )
both<- both[oo1,,]
left<- left[oo1,,]
right<- right[oo1,,]
s<- s[oo1,]
#Sort by number of captures in observed (R) left and right data sets.
#Not sure if this is necessary, inherited from Andy
tl<- apply(left,c(1,3),sum)
tr<- apply(right,c(1,3),sum)
cap.left<- apply(tl,1,sum)[(Nknown+1):N]
cap.right<- apply(tr,1,sum)[(Nknown+1):N]
oo2<- c(1:Nknown, Nknown+ rev(order(cap.left)) )
oo3<- c(1:Nknown, Nknown + rev(order(cap.right)) )
left<- left[oo2,,]
right<- right[oo3,,]
#Record order of true IDs
ID_L=(1:N)[oo2]
ID_R=(1:N)[oo3]
}else{
#Record order of true IDs
ID_L=(1:N)
ID_R=(1:N)
}
#Remove all 0 capture histories. Don't remove right and left all 0 histories for known IDs
both<- both[IDknown,,]
if(length(dim(both))==2){#if there is 1 both guy need to keep 4d
both=array(both,dim=c(1,dim(both)))
}
lcap=which(apply(left,1,sum)>0)
rcap=which(apply(right,1,sum)>0)
left<- left[sort(unique(c(lcap,IDknown))),,]
right<- right[sort(unique(c(rcap,IDknown))),,]
if(length(dim(left))==2){#if there is 1 both guy need to keep 4d
left=array(left,dim=c(1,dim(left)))
}
if(length(dim(right))==2){#if there is 1 both guy need to keep 4d
right=array(right,dim=c(1,dim(right)))
}
lorder=sort(unique(c(lcap,IDknown)))
rorder=sort(unique(c(rcap,IDknown)))
#Update true ID order for left and right for captured guys
ID_L=ID_L[lorder]
ID_R=ID_R[rorder]
#Build all possible observed data sets and count observed spatial recaps
B2D=apply(both,c(1,2),sum)
L2D=apply(left,c(1,2),sum)
R2D=apply(right,c(1,2),sum)
known=dim(both)[1]
if(known>0){
BLR=1*((both+left[1:known,,]+right[1:known,,])>0) #keep boths and lefts and rights for boths
if(length(dim(BLR))==2){
BLR=array(BLR,dim=c(1,dim(BLR)))
}
if(dim(left)[1]>known){
BLRL=abind(BLR,left[(known+1):dim(left)[1],,],along=1) #Add unknown lefts
}else{
BLRL=BLR
}
if(dim(right)[1]>known){
BLRR=abind(BLR,right[(known+1):dim(right)[1],,],along=1) #Add unknown rights
}else{
BLRR=BLR
}
}else{ #no known guys
BLR=both
if(dim(left)[1]>known){
BLRL=left #Add unknown lefts
}else{
BLRL=BLR
}
if(dim(right)[1]>known){
BLRR=right #Add unknown rights
}else{
BLRR=BLR
}
if(length(dim(BLRL))==2){
BLRL=array(BLRL,dim=c(1,dim(BLRL)))
}
if(length(dim(BLRR))==2){
BLRR=array(BLRR,dim=c(1,dim(BLRR)))
}
}
#BLR2D=apply(BLR,c(1,3),sum)
BLRL2D=apply(BLRL,c(1,2),sum)
BLRR2D=apply(BLRR,c(1,2),sum)
if(dim(B2D)[1]>0){
Bscaps=rowSums(1*(B2D>0))
Bscaps[Bscaps>0]=Bscaps[Bscaps>0]-1
Bnscap=sum(Bscaps>0)
Bsumscap=sum(Bscaps)
}else{
Bnscap=0
Bsumscap=0
}
if(dim(L2D)[1]>0){
Lscaps=rowSums(1*(L2D>0))
Lscaps[Lscaps>0]=Lscaps[Lscaps>0]-1
Lnscap=sum(Lscaps>0)
Lsumscap=sum(Lscaps)
}else{
Lnscap=0
Lsumscap=0
}
if(dim(R2D)[1]>0){
Rscaps=rowSums(1*(R2D>0))
Rscaps[Rscaps>0]=Rscaps[Rscaps>0]-1
Rnscap=sum(Rscaps>0)
Rsumscap=sum(Rscaps)
}else{
Rnscap=0
Rsumscap=0
}
if(dim(BLRL2D)[1]>0){
BLRLscaps=rowSums(1*(BLRL2D>0))
BLRLscaps[BLRLscaps>0]=BLRLscaps[BLRLscaps>0]-1
BLRLnscap=sum(BLRLscaps>0)
BLRLsumscap=sum(BLRLscaps)
}else{
BLRLnscap=0
BLRLsumscap=0
}
if(dim(BLRR2D)[1]>0){
BLRRscaps=rowSums(1*(BLRR2D>0))
BLRRscaps[BLRRscaps>0]=BLRRscaps[BLRRscaps>0]-1
BLRRnscap=sum(BLRRscaps>0)
BLRRsumscap=sum(BLRRscaps)
}else{
BLRRnscap=0
BLRRsumscap=0
}
Srecap$Bnscap=Bnscap
Srecap$Bsumscap=Bsumscap
Srecap$Lnscap=Lnscap
Srecap$Lsumscap=Lsumscap
Srecap$Rnscap=Rnscap
Srecap$Rsumscap=Rsumscap
Srecap$BLRLnscap=BLRLnscap
Srecap$BLRLsumscap=BLRLsumscap
Srecap$BLRRnscap=BLRRnscap
Srecap$BLRRsumscap=BLRRsumscap
dim(BLRR)
dim(BLRL)
y.obs=list(BLRL=BLRL,BLRR=BLRR)
out<-list(left=left,right=right,both=both,y.obs=y.obs,s=s,X=X,IDknown=IDknown,n=n,nside=nside,Srecap=Srecap,ID_L=ID_L,ID_R=ID_R, K=K,buff=buff,tf=onoff2,obstype=obstype)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.