R/simCatSPIM.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.

#' Simulate data from the categorical spatial partial identity model (SPIM)
#' @description Will document after I defend. Should be able to get the idea in the example found in the mcmc.CatSPIM help file.
#' @author Ben Augustine
#' @export
simCatSPIM <-
  function(N=120,lam0=0.1,sigma=0.50,K=10,X=X,buff=3,obstype="bernoulli",ncat=ncat,
           pID=pID,gamma=gamma,IDcovs=IDcovs){
    # 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<- lam0*exp(-D*D/(2*sigma*sigma))
    J<- nrow(X)
    
    #simulate IDcovs
    G.true=matrix(NA,nrow=N,ncol=ncat) #all IDcovs in population.
    for(i in 1:N){
      for(j in 1:ncat){
        G.true[i,j]=sample(IDcovs[[j]],1,prob=gamma[[j]])
      }
    }
    if(dim(unique(G.true))[1]!=N){
      print(paste("simulated",
                  length(unique(G.true[duplicated(G.true)]))+sum(duplicated(G.true)),"duplicated IDcovs in population"))
      }
    # Capture individuals
    y <-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){
            y[i,j,k]=rbinom(1,1,pd[i,j]) #single side lambda multiplied by number of traps at site
          }
        }
      }
    }else if(obstype=="poisson"){
      for(i in 1:N){
        for(j in 1:J){
          for(k in 1:K){
            y[i,j,k]=rpois(1,lamd[i,j]) #single side lambda multiplied by number of traps at site
          }
        }
      } 
    }else{
      stop("obstype not recognized")
    }

    #discard uncaptured inds and aggregate true IDcovs for all samples, keeping track of where they came from with A matrix (used to put observed data back together)
    caught=which(apply(y,c(1),sum)>0)
    y.true=y
    y=y[caught,,]
    if(K==1){
      y=array(y,dim=c(dim(y),1))
    }
    n=length(caught)
    n.samples=sum(y)
    G.cap=matrix(NA,nrow=n.samples,ncol=ncat)
    idx=1
    A=array(0,dim=c(dim(y),n.samples))  #configuration matrix: indicator matrix for which individual i occassion j  trap k corresponds to sample l. used to convert corrupt IDcovs to corrupt capture history
    for(i in 1:length(caught)){ #loop through inds (uncaptured already removed)
      for(j in 1:J){ #then traps
        for(k in 1:K){ #then occasions
          if(y[i,j,k]>0){ #is there at least one sample here?
            for(l in 1:y[i,j,k]){ #then samples
              G.cap[idx,]=G.true[caught[i],]
              A[i,j,k,idx]=1
              idx=idx+1
            }
          }
        }
      }
    }
    ycap=aperm(apply(A,c(2,3,4),sum),c(3,1,2))
    #Amplification failure
    G.drop=G.cap #n.samples x ncat
    for(j in 1:ncat){
      G.drop[which(rbinom(n.samples,1,pID[j])==0),j]=0 #0 is dropout
    }
    G.obs=unique(G.drop)
    nobs=nrow(G.obs)
    yobs=array(0,dim=c(nobs,J,K))
    ID=rep(NA,n.samples)
    for(i in 1:n.samples){
      for(j in 1:nobs){
        if(all(G.drop[i,]==G.obs[j,])){
          ID[i]=j
          next
        }
      }
    }
    yobs=array(0,dim=c(max(ID),J,K))
    for(i in 1:n.samples){
      map=which(A[,,,i]==1,arr.ind=TRUE)
      yobs[ID[i],map[2],map[3]]=yobs[ID[i],map[2],map[3]]+1
    }
    
    ycheck=array(0,dim=dim(y))
    for(i in 1:n.samples){
      idx2=which(A[,,,i]>0,arr.ind=TRUE)
      ycheck[idx2[1],,]=ycheck[idx2[1],,]+A[idx2[1],,,i]
    }
    if(!all(ycheck==y)){
      stop("not all y==ycheck")
    }
    #Make ID constraint matrix
    # n.samples=sum(y.obs)
    constraints=matrix(1,nrow=n.samples,ncol=n.samples)
    for(i in 1:n.samples){
      for(j in 1:n.samples){
        guys1=which(G.drop[i,]!=0)
        guys2=which(G.drop[j,]!=0)
        comp=guys1[which(guys1%in%guys2)]
        if(any(G.drop[i,comp]!=G.drop[j,comp])){
          constraints[i,j]=0
        }
      }
    }
    #check constraints
    a=which(constraints==1,arr.ind=TRUE)#consistent
    for(i in 1:nrow(a)){
      comp=G.drop[a[i,1],]>0&G.drop[a[i,2],]>0
      if(!all(G.drop[a[i,1],comp]==G.drop[a[i,2],comp])){
        stop("Error in constraint matrix")
      }
    }
    a=which(constraints==0,arr.ind=TRUE)#inconsistent
    if(length(a)>1){
      for(i in 1:nrow(a)){
        comp=G.drop[a[i,1],]>0&G.drop[a[i,2],]>0
        if(all(G.drop[a[i,1],comp]==G.drop[a[i,2],comp])){
          stop("Error in constraint matrix")
        }
      }
    }
    A=apply(A,c(1,4),sum)
    IDtrue=rep(NA,n.samples)
    for(i in 1:n.samples){
      IDtrue[i]=which(A[,i]==1)
    }
    out<-list(y=y,y.obs=ycap,y.true=y.true,G.true=G.true,G.cap=G.cap,G.obs=G.drop,IDlist=list(ncat=ncat,IDcovs=IDcovs),
              IDtrue=IDtrue,X=X,K=K,buff=buff,constraints=constraints,obstype=obstype,s=s,n=nrow(y))
    return(out)
  }