R/RMark simulation

[SoConfused]

Hello,

I'm trying to simulate a release of fish over 5 non-consecutive years (2002, 2010, 2012, 2013, 2014).
The model is a simple constant phi and p. The results of the simulation are nowhere near my initial values. I get very similar results running the simulation directly in MARK, so there's either something very stupid I do upfront, or there's something fundamental I'm misunderstanding in the process. Any help would be appreciated...

Code: Select all

library(RMark)
# create simulated data (code borrowed from here
#https://sites.google.com/site/wild8390/software/simulate)

simul.cjs <- function(phi, p, marked){
   n.occasions <- length(p) + 1
   Phi <- matrix(phi, n.occasions-1, nrow=sum(marked), byrow=T)
   P <- matrix(p, n.occasions-1, nrow=sum(marked), byrow=T)

   CH <- matrix(0, ncol = n.occasions, nrow = sum(marked))

   #define a vector with marking occasion
   mark.occ <- rep(1:length(marked), marked[1:length(marked)])

   #fill in CH
   for(i in 1:sum(marked)){
      CH[i,mark.occ[i]] <- 1
      if(mark.occ[i] == n.occasions)
         next
      for(t in (mark.occ[i]+1):n.occasions){
              #survive?
              sur <- rbinom(1, 1, Phi[i,t-1])
              if(sur == 0)
            break #move to next
         #recaptured?
         rp <- rbinom(1, 1, P[i,t-1])
            if(rp == 1)
               CH[i,t] <- 1
                             } #t
                 } #i
   return(CH)
                  }


###function to create capture history character strings (need for input to RMARK)
pasty <- function(x){
   k <- ncol(x)
   n <- nrow(x)
   out <- array(dim=n)
   for(i in 1:n){
      out[i] <- paste(x[i,], collapse="")
             }
   return(out)      }


# simulate release of 10,000 fish every year - the number is crazy high to make sure there's no sample size issue in the estimates
marked <- rep(10000, 5)

data <- simul.cjs(rep(0.3, 4), rep(0.2, 4), marked)
data <- pasty(data)
data <- data.frame(ch = data)

data.proc <- process.data(data, model="CJS", begin.time=2002, time.intervals = diff(c(2002, 2010, 2012, 2013, 2014)), initial.ages = 0)

data.ddl <- make.design.data(data.proc)

data.analysis=function(){   
   Phi.1 = list(formula=~1)   
   p.1=list(formula=~1)
   cml=create.model.list("CJS")
   mark.wrapper(cml,data=data.proc,ddl=data.ddl,output=FALSE, retry = 0, delete = TRUE)
           }

data.results <- data.analysis()
data.results$Phi.1.p.1$results$real



[jlaake]

It appears you have your solution but I wanted to mention that the function simHMM in my marked package provides simulation capability for any of the models in that package. Below is an example which simulates the data with what I understood to be your situation and then fits the model with RMark and marked versions of the CJS model. Note that RMark and marked have functions with the same name but do not function exactly the same, so I've used the R notation package::function to run specific functions from each package.

Code: Select all

library(marked)
library(RMark)
#
# simulate data with simHMM in marked package
#
# create 5 cohorts at 2002,2010,2012,2013,2014
df=data.frame(ch=c("1000000000000","0000000010000","0000000000100","0000000000010","0000000000001"),freq=rep(1000,5),stringsAsFactors=FALSE)
# note: marked uses 4 values and throws out the last release which doesn't contain any information for cjs
dp=marked::process.data(df,model="HMMcjs",begin.time=2002)
ddl=marked::make.design.data(dp)
# fix p to 0 for 2003-2009,2011
ddl$p$fix=ifelse(ddl$p$time%in%(c(2003:2009,2011)),0,NA)
# simulate data
data=simHMM(dp,ddl,initial=list(Phi=-0.8,p=0.5))
#
# fit model in RMark
#
# RMark does not accept "," in capture history which marked uses as defaults so remove those
markdata=data
markdata$ch=apply(marked::splitCH(markdata$ch),1,paste,collapse="")
# process data and create design data with RMark functions
dp=RMark::process.data(markdata,begin.time=2002)
ddl=RMark::make.design.data(dp)
# fix p=0 for the years with no sampling
ddl$p$fix=ifelse(ddl$p$time%in%(c(2003:2009,2011)),0,NA)
# fit model and show coefficients
model=RMark::mark(dp,ddl)
coef(model)

# fit model with marked
dp=marked::process.data(data,begin.time=2002)
ddl=marked::make.design.data(dp)
# fix p=0 for the years with no sampling
ddl$p$fix=ifelse(ddl$p$time%in%(c(2003:2009,2011)),0,NA)
model=crm(dp,ddl)
coef(model)