R/mcmc.2sidetf.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.

mcmc.2sidetf <-
  function(data,niter=2400,nburn=1200, nthin=5, M = 200, inits=inits,
           swap=10,swap.tol=1,proppars=list(lam01=0.05,lam02=0.05,sigma=0.1,sx=0.2,sy=0.2),storeLatent=TRUE){
    library(abind)
    y.both<-data$both
    y.left.obs<-data$left
    y.right.obs<-data$right
    if(length(dim(y.both))==3){
      y.both=apply(y.both,c(1,2),sum)
    }
    if(length(dim(y.left.obs))==2){
      stop("Must use 3D single side data with 2D trap file")
    }
    if(dim(y.right.obs)[1]>dim(y.left.obs)[1]){
      storeL=y.left.obs
      storeR=y.right.obs
      dimL=dim(y.left.obs)
      y.left.obs=array(0,dim=dim(y.right.obs))
      y.right.obs=array(0,dim=dimL)
      y.left.obs=storeR
      y.right.obs=storeL
      warning("Right side data set larger than left so I switched them for convenience")
    }
    if(!("tf"%in%names(data))){
      stop("This function requires a 2D trap file. Tell Ben there is an error")
    }
    if(!is.matrix(data$tf)){
      stop("This is the wrong function for a 2D trap file. Tell Ben something went wrong")
    }
    tf=data$tf
    tf1=abind(tf,tf,along=0)
    for(i in 3:M){
      tf1=abind(tf1,tf,along=1)
    }
    tf2=rowSums(tf==2)
    tf2=matrix(rep(tf2,M),nrow=M,byrow=TRUE)
    X<-as.matrix(data$X)
    J<-nrow(X)
    K<- data$K
    IDknown<- data$IDknown
    Nfixed=length(IDknown)
    nleft<-dim(y.left.obs)[1]-Nfixed
    nright<-dim(y.right.obs)[1]-Nfixed
    if("vertices"%in%names(data)){
      vertices=data$vertices
      useverts=TRUE
      xlim=c(min(vertices[,1]),max(vertices[,1]))
      ylim=c(min(vertices[,2]),max(vertices[,2]))
    }else if("buff"%in%names(data)){
      buff<- data$buff
      xlim<- c(min(X[,1]),max(X[,1]))+c(-buff, buff)
      ylim<- c(min(X[,2]),max(X[,2]))+c(-buff, buff)
      vertices=rbind(xlim,ylim)
      useverts=FALSE
    }else{
      stop("user must supply either 'buff' or 'vertices' in data object")
    }
    ##pull out initial values
    psi<- inits$psi
    lam01<- inits$lam01
    sigma<- inits$sigma
    lam01<- inits$lam01
    lam02<- inits$lam02
    
    #Figure out what needs to be updated
    uplam01=uplam02=upIDs=TRUE
    if(lam01==0){
      uplam01=FALSE
      upIDs=FALSE
    }
    if(all(X[,3]==1)|lam02==0){
      uplam02=FALSE
    }
    if(upIDs==TRUE&(nleft==0&nright==0)){  #not sure what will happen if we only have left only or right only
      upIDs=FALSE
    }
    
    #augment 3 data sets
    y.both<- abind(y.both,array(0, dim=c( M-dim(y.both)[1],J)), along=1)
    y.left.obs<- abind(y.left.obs,array(0, dim=c( M-dim(y.left.obs)[1],J,K)), along=1)
    y.right.obs<-abind(y.right.obs,array(0,dim=c( M-dim(y.right.obs)[1],J,K)), along=1)
    
    #sort to minimize distance between initial matches. Skip if no single sides.
    # if(nleft>0|nright>0){
    if(nleft>0&nright>0){ #change 1/14/19
      y.left.tmp=apply(y.left.obs,c(1,2),sum)
      y.right.tmp=apply(y.right.obs,c(1,2),sum)
      IDs<- LRmatch(M=M,left=y.left.tmp, nleft=nleft, right=y.right.tmp, nright=nright, X, Nfixed=Nfixed)
      #Add unused augmented indivuals back in
      notusedL<- (1:M)[is.na(match(1:M,IDs$ID_L))]
      ID_L<-c(IDs$ID_L,notusedL)
      notusedR<- (1:M)[is.na(match(1:M,IDs$ID_R))]
      ID_R<-c(IDs$ID_R,notusedR)
    }else{
      ID_R=ID_L=1:M
    }
    
    #reoder left and right possibly true data sets
    y.left.true<- y.left.obs[order(ID_L),,]
    y.right.true<- y.right.obs[order(ID_R),,]
    
    #Make initial complete data set
    tmpdata<- y.both + apply(y.left.true+y.right.true,c(1,2),sum)
    z=1*(apply(tmpdata,1,sum)>0)
    z[sample(which(z==0),sum(z==0)*psi)]=1 #switch some uncaptured z's to 1.
    
    #Optimize starting locations given where they are trapped.
    s<- cbind(runif(M,xlim[1],xlim[2]), runif(M,ylim[1],ylim[2])) #assign random locations
    idx=which(rowSums(tmpdata)>0) #switch for those actually caught
    for(i in idx){
      trps<- X[tmpdata[i,]>0,1:2]
      trps<-matrix(trps,ncol=2,byrow=FALSE)
      s[i,]<- c(mean(trps[,1]),mean(trps[,2]))
    }
    if(useverts==TRUE){
      inside=rep(NA,nrow(s))
      for(i in 1:nrow(s)){
        inside[i]=inout(s[i,],vertices)
      }
      idx=which(inside==FALSE)
      if(length(idx)>0){
        for(i in 1:length(idx)){
          while(inside[idx[i]]==FALSE){
            s[idx[i],]=c(runif(1,xlim[1],xlim[2]), runif(1,ylim[1],ylim[2]))
            inside[idx[i]]=inout(s[idx[i],],vertices)
          }
        }
      }
    }
    known.vector<- c( rep(1,Nfixed), rep(0, M-Nfixed) )
    
    # some objects to hold the MCMC simulation output
    nstore=(niter-nburn)/nthin
    if(nburn%%nthin!=0){
      nstore=nstore+1
    }
    out<-matrix(NA,nrow=nstore,ncol=4)
    dimnames(out)<-list(NULL,c("lam01","lam02","sigma","N"))
    if(storeLatent){
      sxout<- syout<- zout<- ID_Lout<- ID_Rout<-matrix(NA,nrow=nstore,ncol=M)
    }
    idx=1 #for storing output not recorded every iteration
    
    D<- e2dist(s, X)
    lamd1<- lam01*exp(-D*D/(2*sigma*sigma))
    
    lamd2<- lam02*exp(-D*D/(2*sigma*sigma))
    pd1=1-exp(-lamd1)
    pd12=2*pd1-pd1*pd1
    pd2=1-exp(-lamd2)
    lamd1.cand=lamd1
    lamd2.cand=lamd2
    pd1.cand=pd1
    pd2.cand=pd2
    zero.guys<- apply(y.both,1,sum)+apply(y.left.true + y.right.true ,1,sum) == 0
    ones=1*(tf1==1)
    twos=1*(tf1==2)
    sumones=apply(ones,c(1,2),sum)
    sumtwos=apply(twos,c(1,2),sum)
    ll.y.both=dbinom(y.both,tf2,z*pd2,log=TRUE)
    ll.y.left=dbinom(y.left.true*ones,ones,z*pd1,log=TRUE)+
      dbinom(y.left.true*twos,twos,z*pd12,log=TRUE)
    ll.y.right=dbinom(y.right.true*ones,ones,z*pd1,log=TRUE)+
      dbinom(y.right.true*twos,twos,z*pd12,log=TRUE)
    if(!is.finite(sum(ll.y.both)))stop("Both side likelihood not finite. Make sure all camera stations recording both side captures have 2 cameras. Then try changing lam02 or sigma inits.")
    if(!is.finite(sum(ll.y.left)))stop("Left side likelihood not finite. Try changing lam01 or sigma inits.")
    if(!is.finite(sum(ll.y.right)))stop("right side likelihood not finite. Try changing lam01 or sigma inits.")
    
    ll.y.both.cand=ll.y.both
    ll.y.left.cand=ll.y.left
    ll.y.right.cand=ll.y.right
    
    for(iter in 1:niter){
      #Update lam01
      lly1sum=sum(ll.y.left)+sum(ll.y.right)
      if(uplam01){
        lam01.cand<- rnorm(1,lam01,proppars$lam01)
        if(lam01.cand > 0){
          lamd1.cand<- lam01.cand*exp(-D*D/(2*sigma*sigma))
          pd1.cand=1-exp(-lamd1.cand)
          pd12.cand=2*pd1.cand-pd1.cand*pd1.cand
          ll.y.left.cand=dbinom(y.left.true*ones,ones,z*pd1.cand,log=TRUE)+
            dbinom(y.left.true*twos,twos,z*pd12.cand,log=TRUE)
          ll.y.right.cand=dbinom(y.right.true*ones,ones,z*pd1.cand,log=TRUE)+
            dbinom(y.right.true*twos,twos,z*pd12.cand,log=TRUE)
          lly1candsum=sum(ll.y.left.cand)+sum(ll.y.right.cand)
          if(runif(1) < exp(lly1candsum - lly1sum)){
            lam01=lam01.cand
            lamd1=lamd1.cand
            pd1=pd1.cand
            pd12=pd12.cand
            ll.y.left=ll.y.left.cand
            ll.y.right=ll.y.right.cand
            lly1sum=lly1candsum
          }
        }
      }
      #Update lam02
      lly2sum=sum(ll.y.both)
      if(uplam02){
        lam02.cand<- rnorm(1,lam02,proppars$lam02)
        if(lam02.cand > 0){
          lamd2.cand<- lam02.cand*exp(-D*D/(2*sigma*sigma))
          pd2.cand=1-exp(-lamd2.cand)
          ll.y.both.cand <- dbinom(y.both,tf2,z*pd2.cand,log=TRUE)
          lly2candsum=sum(ll.y.both.cand)
          if(runif(1) < exp(lly2candsum-lly2sum)){
            lam02=lam02.cand
            lamd2=lamd2.cand
            pd2=pd2.cand
            ll.y.both=ll.y.both.cand
            lly2sum=lly2candsum
          }
        }
      }
      #Update sigma
      sigma.cand<- rnorm(1,sigma,proppars$sigma)
      if(sigma.cand > 0){
        lamd1.cand<- lam01*exp(-D*D/(2*sigma.cand*sigma.cand))
        lamd2.cand<- lam02*exp(-D*D/(2*sigma.cand*sigma.cand))
        pd1.cand=1-exp(-lamd1.cand)
        pd2.cand=1-exp(-lamd2.cand)
        pd12.cand=2*pd1.cand-pd1.cand*pd1.cand
        ll.y.left.cand=dbinom(y.left.true*ones,ones,z*pd1.cand,log=TRUE)+
          dbinom(y.left.true*twos,twos,z*pd12.cand,log=TRUE)
        ll.y.right.cand=dbinom(y.right.true*ones,ones,z*pd1.cand,log=TRUE)+
          dbinom(y.right.true*twos,twos,z*pd12.cand,log=TRUE)
        lly1candsum=sum(ll.y.left.cand)+sum(ll.y.right.cand)
        ll.y.both.cand <- dbinom(y.both,tf2,z*pd2.cand,log=TRUE)
        lly2candsum=sum(ll.y.both.cand)
        if(runif(1) < exp((lly1candsum+lly2candsum) -(lly1sum+lly2sum))){
          sigma<- sigma.cand
          lamd1=lamd1.cand
          lamd2=lamd2.cand
          pd1=pd1.cand
          pd12=pd12.cand
          pd2=pd2.cand
          ll.y.both=ll.y.both.cand
          ll.y.left=ll.y.left.cand
          ll.y.right=ll.y.right.cand
        }
      }
      
      ## Update latent ID variables
      ## Candidate: swap IDs of one guy with some other guy. The two guys to be swapped are
      ## chosen randomly from the z=1 guys
      if(any(z[!zero.guys]==0)){
        cat("some z = 0!",fill=TRUE)
        browser()
      }
      #Swap left guys first
      if(nleft>0){
        # User inputs how many swaps to make on each iteration my specifying "swap"
        #map lefts to boths
        #candmap used to remove disallowed candidates.NA any lefts that map back to z=0 indices and boths that are z=0 indices
        map=cbind(1:M,ID_L)
        candmap=map
        candmap[1:Nfixed,]=NA #Don't swap out IDknown guys
        candmap[z==0,1]=NA #Don't swap in z=0 guys.
        candmap[candmap[,2]%in%which(z==0),2]=NA #Don't choose a guy1 that will make you swap in a z=0 guy
        #These are the guys that can come out and go in
        OUTcands=which(!is.na(candmap[,2]))
        INcands=which(!is.na(candmap[,1]))
        for(up in 1:swap){
          ### In this code here "guy1" and "guy2" are indexing "left-side encounter histories" whereas
          ### "s.swap.in" and "s.swap.out" are indexing (both-side guys, s, z)
          guy1<- sample(OUTcands, 1)
          s.swap.out<- map[guy1,2]
          # to find candidates for swapping look in vicinity of it's current both side membership
          dv<-  sqrt( (s[s.swap.out,1]- s[INcands,1])^2 + (s[s.swap.out,2] - s[INcands,2])^2 )
          ## This is a list of both side activity centers that could be assinged to the right side guy
          ## under consideration based on how far away they are.  swap.tol is the spatial tolerance
          # if no one around, code swaps guy 1 with guy 1 (no swap, but does calculations)
          possible<- INcands[dv < swap.tol]
          jump.probability<- 1/length(possible)  # h(theta*|theta)
          # this is a particular value of s to swap for guy1.id
          if(length(possible)>1){
            s.swap.in <-  sample( possible, 1)
          }
          if(length(possible) ==1){
            s.swap.in <- possible
          }
          if(length(possible)==0) next #no guys close enough to swap

          #compute  h(theta|theta*)
          trash<-   sqrt( (s[s.swap.in,1]- s[INcands,1])^2 + (s[s.swap.in,2] - s[INcands,2])^2  )
          trash<-  INcands[trash < swap.tol]
          jump.back<-  1/length(trash)

          ##  Which left encounter history is currently associated with both guy s.swap.in?
          guy2<- which(map[,2]==s.swap.in)
          if(guy1==guy2) next #Same guy selected
          newID<-ID_L
          newID[guy1]<- s.swap.in
          newID[guy2]<- s.swap.out

          ## recompute 'data' and compute likelihood components for swapped guys only
          y.left.tmp<- y.left.obs[order(newID),,]#Reorder observed left side data set
          swapped=c(s.swap.out,s.swap.in)
          ll.y.left.tmp=dbinom(y.left.tmp[swapped,,]*ones[swapped,,],ones[swapped,,],pd1[swapped,],log=TRUE)+
            dbinom(y.left.tmp[swapped,,]*twos[swapped,,],twos[swapped,,],pd12[swapped,],log=TRUE)

          #MH step
          if(runif(1)<exp(sum(ll.y.left.tmp)-sum(ll.y.left[swapped,,]))*(jump.back/jump.probability) ){
            y.left.true=y.left.tmp #update left data
            ll.y.left[swapped,,]=ll.y.left.tmp
            ID_L<-newID #update data order
            map[c(guy1,guy2),2]=c(s.swap.in,s.swap.out)
            zero.guys<- apply(y.both,1,sum)+apply(y.left.true + y.right.true ,1,sum) == 0
            if( sum(y.both[guy1,] + y.left.tmp[guy1,,] )>0 & z[guy1]==0){
              cat("error on guy1",fill=TRUE)
              browser()
            }
            if( sum(y.both[guy2,] + y.left.tmp[guy2,,] )>0 & z[guy2]==0){
              cat("error on guy2",fill=TRUE)
              browser()
            }
          }
        }
        #Do any captured guys have z==0?
        if(any(z[!zero.guys]==0)){
          cat("coming out of ID_L update: some z = 0!",fill=TRUE)
          browser()
        }
      }

      #Repeat for rights
      if(nright>0){
        #Reuse left side of maps
        map[,2]=candmap[,2]=ID_R
        candmap[1:Nfixed,2]=NA
        candmap[candmap[,2]%in%which(z==0),2]=NA
        OUTcands=which(!is.na(candmap[,2])) #INcands the same
        for(up in 1:swap){
          ### In this code here "guy1" and "guy2" are indexing "right-side encounter histories" whereas
          ### "s.swap.in" and "s.swap.out" are indexing (both-side guys, s, z)
          guy1<- sample(OUTcands, 1)
          s.swap.out<- map[guy1,2]
          # to find candidates for swapping look in vicinity of it's current both side membership
          dv<-  sqrt( (s[s.swap.out,1]- s[INcands,1])^2 + (s[s.swap.out,2] - s[INcands,2])^2 )
          ## This is a list of both side activity centers that could be assinged to the right side guy
          ## under consideration based on how far away they are.  swap.tol is the spatial tolerance
          possible<- INcands[dv < swap.tol]
          jump.probability<- 1/length(possible)  # h(theta*|theta)
          # this is a particular value of s to swap for guy1.id
          if(length(possible)>1){
            s.swap.in <-  sample( possible, 1)
          }
          if(length(possible) ==1){
            s.swap.in <- possible
          }
          if(length(possible)==0) next #no guys close enough to swap

          #compute  h(theta|theta*)
          trash<-   sqrt( (s[s.swap.in,1]- s[INcands,1])^2 + (s[s.swap.in,2] - s[INcands,2])^2  )
          trash<-  INcands[trash < swap.tol]
          jump.back<-  1/length(trash)

          ##  Which right encounter history is currently associated with both guy s.swap.in?
          guy2<- which(map[,2]==s.swap.in)
          if(guy1==guy2) next #Same guy selected
          newID<-ID_R
          newID[guy1]<- s.swap.in
          newID[guy2]<- s.swap.out

          ## recompute 'data' and compute likelihood components for swapped guys only
          y.right.tmp<- y.right.obs[order(newID),,]#Reorder observed right side data set
          swapped=c(s.swap.out,s.swap.in)
          ll.y.right.tmp=dbinom(y.right.tmp[swapped,,]*ones[swapped,,],ones[swapped,,],pd1[swapped,],log=TRUE)+
            dbinom(y.right.tmp[swapped,,]*twos[swapped,,],twos[swapped,,],pd12[swapped,],log=TRUE)
          #MH step
          if(runif(1)<exp(sum(ll.y.right.tmp)-sum(ll.y.right[swapped,,]))*(jump.back/jump.probability) ){
            y.right.true=y.right.tmp #update right data
            ll.y.right[swapped,,]=ll.y.right.tmp
            ID_R<-newID #update data order
            map[c(guy1,guy2),2]=c(s.swap.in,s.swap.out)
            zero.guys<- apply(y.both,1,sum)+apply(y.left.true + y.right.true ,1,sum) == 0
            if( sum(y.both[guy1,] + y.right.tmp[guy1,,] )>0 & z[guy1]==0){
              cat("error on guy1",fill=TRUE)
              browser()
            }
            if( sum(y.both[guy2,] + y.right.tmp[guy2,,] )>0 & z[guy2]==0){
              cat("error on guy2",fill=TRUE)
              browser()
            }
          }
        }
        #Do any captured guys have z==0?
        if(any(z[!zero.guys]==0)){
          cat("coming out of ID_R update: some z = 0!",fill=TRUE)
          browser()
        }
      }
      #Update psi gibbs
      ## probability of not being captured in a trap AT ALL
      pbar1=(1-pd1)^(2*sumones)*(1-pd12)^(2*sumtwos) #can be a left or a right
      pbar2=(1-pd2)^tf2
      pbar0=pbar1*pbar2
      prob0<- exp(rowSums(log(pbar0)))
      fc<- prob0*psi/(prob0*psi + 1-psi)
      z[zero.guys & known.vector==0]<- rbinom(sum(zero.guys & (known.vector ==0) ), 1, fc[zero.guys & (known.vector==0) ])
      #update observation model ll
      ll.y.both=dbinom(y.both,tf2,z*pd2,log=TRUE)
      ll.y.left=dbinom(y.left.true*ones,ones,z*pd1,log=TRUE)+
        dbinom(y.left.true*twos,twos,z*pd12,log=TRUE)
      ll.y.right=dbinom(y.right.true*ones,ones,z*pd1,log=TRUE)+
        dbinom(y.right.true*twos,twos,z*pd12,log=TRUE)
      psi <- rbeta(1, 1 + sum(z), 1 + M - sum(z))

      # ## Now we have to update the activity centers
      for (i in 1:M) {
        Scand <- c(rnorm(1, s[i, 1], proppars$sx), rnorm(1, s[i, 2], proppars$sy))
        if(useverts==FALSE){
          inbox <- Scand[1] < xlim[2] & Scand[1] > xlim[1] & Scand[2] < ylim[2] & Scand[2] > ylim[1]
        }else{
          inbox=inout(Scand,vertices)
        }
        if (inbox) {
          dtmp <- sqrt((Scand[1] - X[, 1])^2 + (Scand[2] - X[, 2])^2)
          lamd1.cand[i,]=lam01*exp(-dtmp*dtmp/(2*sigma*sigma))
          lamd2.cand[i,]=lam02*exp(-dtmp*dtmp/(2*sigma*sigma))
          pd1.cand[i,]=1-exp(-lamd1.cand[i,])
          pd12.cand[i,]=2*pd1.cand[i,]-pd1.cand[i,]*pd1.cand[i,]
          pd2.cand[i,]=1-exp(-lamd2.cand[i,])
          ll.y.both.cand[i,]=dbinom(y.both[i,],tf2[i,],z[i]*pd2.cand[i,],log=TRUE)
          
          ll.y.left.cand[i,,]=dbinom(y.left.true[i,,]*ones[i,,],ones[i,,],z[i]*pd1.cand[i,],log=TRUE)+
            dbinom(y.left.true[i,,]*twos[i,,],twos[i,,],z[i]*pd12.cand[i,],log=TRUE)
          
          ll.y.right.cand[i,,]=dbinom(y.right.true[i,,]*ones[i,,],ones[i,,],z[i]*pd1.cand[i,],log=TRUE)+
            dbinom(y.right.true[i,,]*twos[i,,],twos[i,,],z[i]*pd12.cand[i,],log=TRUE)
          llysum=(sum(ll.y.both[i,])+sum(ll.y.left[i,,])+sum(ll.y.right[i,,]))
          llycandsum=(sum(ll.y.both.cand[i,])+sum(ll.y.left.cand[i,,])+sum(ll.y.right.cand[i,,]))
          if (runif(1) < exp(llycandsum-llysum)) {
            s[i, ]=Scand
            D[i, ]=dtmp
            lamd1[i, ]=lamd1.cand[i,]
            lamd2[i, ]=lamd2.cand[i,]
            pd1[i,]=pd1.cand[i,]
            pd12[i,]=pd12.cand[i,]
            pd2[i,]=pd2.cand[i,]
            ll.y.both[i,]=ll.y.both.cand[i,]
            ll.y.left[i,,]=ll.y.left.cand[i,,]
            ll.y.right[i,,]=ll.y.right.cand[i,,]
          }
        }
      }
      #Do we record output on this iteration?
      if(iter>nburn&iter%%nthin==0){
        if(storeLatent){
          sxout[idx,]<- s[,1]
          syout[idx,]<- s[,2]
          zout[idx,]<- z
          ID_Lout[idx,]<-ID_L
          ID_Rout[idx,]<-ID_R
        }
        out[idx,]<- c(lam01,lam02,sigma ,sum(z))
        idx=idx+1
      }
    }  # end of MCMC algorithm
    if(storeLatent==TRUE){
      list(out=out, sxout=sxout, syout=syout, zout=zout, ID_Lout=ID_Lout,ID_Rout=ID_Rout)
    }else{
      list(out=out)
    }
  }