MARK simulations - CJS data

[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. I ask for an identity link of the simulated data and a sin link of the analysis, with 50 sims and an output of the real values and their SE. My release numbers are 20600; 9117; 3981; 4328; 14.

Depending on the phi and p I choose, the results of the simulation are nowhere near my initial values.

For example - Phi = 0.7, p = 0.02; real estimated Phis range between 0.941 and 0.96, real estimated p range between 0.012 and 0.015.
Phi = 0.3, p = 0.02; real estimated Phis range between 0.862 and 0.908, real estimated p range between 0.005 and 0.007.
Phi = 0.3, p = 0.1; real Phis range between 0.871 and 0.906, real p values range between 0.031 and 0.042...

I posted this on the RMark forum as well, but had no replies, so reposting as a MARK-only question (viewtopic.php?f=21&t=3124)

Any help would be appreciated...

[cooch]

Hello,

I'm trying to simulate a release of fish over 5 non-consecutive years (2002, 2010, 2012, 2013, 2014).

Which is undoubtedly the root of the problem -- you have fairly large, unequal intervals between sampling occasions. You should consult the - sidebar - in Chapter 4 (pp. 24-28) before going much further. Time-invariant models are somewhat complicated with unequal intervals.

[jhines]

You're simulating data assuming 1 year between samples, but analyzing data assuming 8,2,1,1 years between samples. To correct the simulation:

Code: Select all

simul.cjs <- function(phi, p, marked){
  n.occasions <- length(p) + 1
  Phi <- matrix(phi, n.occasions-1, nrow=sum(marked), byrow=T)
  intervals=diff(c(2002,2010,2012,2013,2014))
  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]^intervals[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)
}

[SoConfused]

Holy, must feel great to be right all the time

I reran the sims with 1-year intervals, and the numbers came back perfectly. Question is, where do I go from here? The idea of this thing was to estimate variability of abundance estimates (based on # of released hatchery fish in each year and estimated survival) as a function of assumed phi and p.

A naive question - can I simulate the data using 1-year intervals (since all I care about is the spread of survival estimates anyway), and then apply the survival estimates using the # of year intervals? So, if I simulate survival of 0.3 with 1-year intervals, I would then calculate 20,600 (number released in 2002) times 0.3^8 for the first step?

Alternatively, do I run the simulation using 3 phi values, where the first one would be constant phi/8, then phi/2, and then phi, and run from there?

Any intelligent suggestions would be most welcome...

[SoConfused]

You're simulating data assuming 1 year between samples, but analyzing data assuming 8,2,1,1 years between samples.

Thanks so much for pointing that out and altering the code, you're perfectly right. I reran the RMark bit, and am getting good estimates now...

Thank you for your help!

[jlaake]

Sounds like you have a solution. It is always good to simulate a model to understand it. If others want to simulate data from cjs and some other models with an existing function in R, see my reply to the original post on the RMark subforum.