Simulation of a Northern fur seal CR data set

NFS_CR_simulation.R

library("marked")

#### Some preliminaries ###
ilogit = function(x){log(x/(1-x))}
min.adult.age=5
additional.resight=3
adult.age="stable"

set.seed(111)
n.rel.pup=c(138,57,38,250,rep(100,4))
n.rel.ad=c(31,94,44,40, rep(20,4))
  # Release schedule for pups
  release.schedule=diag(length(n.rel.pup))
  ch.char=apply(release.schedule, 1, paste, collapse="")
  additional.resight=floor(additional.resight)
  if(additional.resight>0) release.schedule=cbind(release.schedule, matrix(0,nrow(release.schedule), additional.resight))
  simd=data.frame(ch=ch.char, freq=n.rel.pup, initial.age=0, stringsAsFactors=FALSE)
  # Release data for adults
  if(adult.age=="stable"){
    stable.dist=c(12,11,10,9,9,8,7,6,6,5,4,4,3,2,2,1,.5,.4,.3,.2,.1)
    ad.ages=sample(c(5:25),sum(n.rel.ad), prob=stable.dist, replace=TRUE)
    simd.ad=data.frame(ch=rep(ch.char, n.rel.ad), freq=1, initial.age=ad.ages, stringsAsFactors=FALSE)
  } else {
    simd.ad=data.frame(ch=ch.char, freq=n.rel.ad, initial.age=min.adult.age, stringsAsFactors=FALSE)
  }
  simd=rbind(simd, simd.ad)
  sd=process.data(simd, model="hmmCJS")
  ddl=make.design.data(sd)
  ddl$Phi$ageClass=cut(ddl$Phi$Age, c(-999,2,3,4,999), 
                       include.lowest=TRUE, right=FALSE,
                       labels=c("pup/juv1","juv2","juv3","ad"))
  ddl$p$ageClass=cut(ddl$p$Age, c(-999,2,3,4,999), 
                     include.lowest=TRUE, right=FALSE,
                     labels=c("juv1","juv2","juv3","ad"))
  ### FIX AGE 1 DETECTION PROB TO ZERO ###
  ddl$p$fix = ifelse(ddl$p$ageClass=="juv1",0,NA)
  ########################################
  initial=list(
    Phi=c(ilogit(sqrt(0.25)), ilogit(0.64), ilogit(0.73), ilogit(0.77)), 
    p=c(ilogit(0.2), ilogit(0.67), ilogit(0.92))
    #p=c(ilogit(0.001), ilogit(0.999), ilogit(0.999), ilogit(0.999))
  )
  model.spec=list(p=list(formula=~ageClass-1), Phi=list(formula=~ageClass-1))
  rel=simHMM(sd, ddl, model.parameters=model.spec, initial=initial)
  model.spec
  
  #Analyze simulated data
  reld=process.data(rel, model="CJS")
  ddl.rel=make.design.data(reld)
  ddl.rel$Phi$ageClass=cut(ddl.rel$Phi$Age, c(-999,2,3,4,999), 
                           include.lowest=TRUE, right=FALSE,
                           labels=c("pup/juv1","juv2","juv3","ad"))
  ddl.rel$p$ageClass=cut(ddl.rel$p$Age, c(-999,2,3,4,999), 
                         include.lowest=TRUE, right=FALSE,
                         labels=c("juv1","juv2","juv3","ad"))
                         
  ### FIX AGE 1 DETECTION PROB TO ZERO ###
  #ddl.rel$p$fix = ifelse(ddl.rel$p$ageClass=="juv1",0,NA)
  ########################################
                                   
  fit=crm(reld, ddl=ddl.rel, model.parameters=model.spec, hessian=TRUE)
  fit$results$reals