R/sim2side.R

In benaug/SPIM: This package contains MCMC algorithms in R and (some) Rcpp for various Spatial Partial Identity Models. Prospective users should email Ben before using.

e2dist<-function (x, y)
{
  i <- sort(rep(1:nrow(y), nrow(x)))
  dvec <- sqrt((x[, 1] - y[i, 1])^2 + (x[, 2] - y[i, 2])^2)
  matrix(dvec, nrow = nrow(x), ncol = nrow(y), byrow = F)
}


#' Simulate data from camera trap SCR study with partial identity
#' @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 J x 3 matrix of trap locations and number of cameras at each trap. columns are X, Y, #cams (1 or 2)
#' @param buff the distance to buffer the trapping array in the X and Y dimensions to produce the state space
#' @param obstype observation type, either "bernoulli" or "poisson"
#' @param tf an optional J x K matrix of dynamic trap operation indicating when and where 0, 1, or 2 cams were deployed
#' @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 This function simulates data from a camera trap SCR study where some or all identities are partial.  The trap
#' matrix, X, should contain the X and Y locations of traps in columns 1 and 2 and the number of cameras at each
#' site in column 3 (either 1 or 2).  The extent of the state space is controlled by "buff", which buffers the
#' minimum and maximum X and Y extents.  Therefore, it probably only makes sense for square or rectangular grids.  Functionality
#' for user-provided polygonal state spaces will be added in the future.
#' @author Ben Augustine
#' @export

sim2side <-
  function(N=50,lam01=0.1,lam02=0.2,sigma=0.50,K=10,X=X,buff=3,obstype="bernoulli",tf=NA){
    library(abind)
    J<- nrow(X)
    #Check for user errors (add more later)
    if(ncol(X)!=3){
      stop("X must have 3 columns, X, Y, # cameras")
    }
    if(!is.null(dim(tf))){
      if(!all(dim(tf)==c(J,K))){
        stop("tf must be of dimension J x K")
      }
      if(!all(tf==1|tf==2|tf==0)){
        stop("elements of tf must be 0, 1, or 2")
      }
      twos=which(apply(tf,1,function(x){any(x==2)}))
      if(!all(twos%in%which(X[,3]==2))){
        stop("trap object X must have a 2 in 3rd column if tf indicates it had 2 cams operational at some point during survey")
      }
      if(any(twos%in%which(X[,3]==1))){
        stop("trap object must have a 1 in 3rd column if tf indicates it never had 2 cams operational")
      }
      usetf=TRUE
    }else{
      usetf=FALSE
      tf=matrix(rep(X[,3],K),ncol=K)
    }
    #######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)

    # 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,(tf[j,k]==1)*pd[i,j,1]+(tf[j,k]==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,(tf[j,k]==1)*pd[i,j,1]+(tf[j,k]==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])*(tf[j,k]==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,tf[j,k]*lamd[i,j,1]) #single side lambda multiplied by number of traps at site
            right[i,j,k]=rpois(1,tf[j,k]*lamd[i,j,1]) #single side lambda multiplied by number of traps at site
            both[i,j,k]=rpois(1,lamd[i,j,2])*(tf[j,k]==2)*1  #both side lambda multiplied by indicator for 2 traps at site
          }
        }
      }
    }else{
      stop("observation model not recognized")
    }

    ######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,2),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,2),sum)
      tr<- apply(right,c(1,2),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 3d
      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 3d
      left=array(left,dim=c(1,dim(left)))
    }
    if(length(dim(right))==2){#if there is 1 both guy need to keep 3d
      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)))
      }
    }
    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)
    if(usetf==TRUE){
      out<-list(left=left,right=right,both=both,y.obs=y.obs,s=s,X=X,tf=tf,IDknown=IDknown,n=n,nside=nside,Srecap=Srecap,ID_L=ID_L,ID_R=ID_R, K=K,buff=buff,obstype=obstype)
    }else{
      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,obstype=obstype)
    }
    return(out)
  }