Pradel f with time intervals

Author

Hines

Background

The goal is to determine how Mark handles the recruitment parameter, \(f\) with respect to the intervals between occasions. As with survival and seniority, we would expect recruitment to be affected by the time between capture occasions. To test how Mark handles the recruitment parameter, we’ll generate expected value data for 6 occasions, assuming equal time intervals. Then we’ll throw out one of the occasions and analyze the resulting data with a time-specific model for survival (\(\phi\)) and recruitment (\(f\)).

Generate expected value data

Here is a recursive function to generate expected-value capture-histories

cap <- function(ch, n, per) {
   ncap=n*p; notcap=n-ncap
   if (per<nper) {
      nalive=notcap*S; ndead=notcap-nalive;
      if (is.na(nx[ch])) nx[ch]<<-0
      nx[ch] <<- nx[ch]+ndead
      if (nalive>0) cap(ch,nalive,per+1)
      nalive=ncap*S; ndead=ncap-nalive; substr(ch,per,per)='1'
      if (is.na(nx[ch])) nx[ch]<<-0
      nx[ch]<<-nx[ch]+ndead
      if (nalive>0) cap(ch, nalive, per+1)
   } else {
     if (is.na(nx[ch])) nx[ch]<<-0
     nx[ch]<<-nx[ch]+notcap
     substr(ch,per,per)='1'
     if (is.na(nx[ch])) nx[ch]<<-0
     nx[ch]<<-nx[ch]+ncap
   }
}

Here is the code to set the desired parameters and call the function to generate capture-histories.

N=100               # number of individuals at start (occasion 1)
S=0.6               # survival rate
f=0.2               # recruitment (no. new animals in t+1 per old animal in t)
nper=6              # number of capture occasions
p=0.4               # capture probability

ch0=paste(rep(0,nper),collapse='');        # initialize capture-history to '000000'
nx=nx1=nx2=nx3=nx4=nx5=nx6=rep(0,2^nper)   # initialize cap-hist frequencies

#  name the cap-hist frequencies with cap-histories ('000000','000001','000010',...)
nxnames=sapply(0:(2^nper-1),
    function(i)  paste(rev(as.integer(intToBits((i)))[1:nper]),collapse=''))
names(nx)=nxnames

cap(ch0,N,1);       # generate histories for initial cohort of N animals
nx1=nx              # function allways saves cap-hist frequencies as nx, 
                    #   so save frequencies as nx1
B0=N                # save number entering before 1st occasion as B0

Next, we’ll call the function again to generate histories for new individuals entering the population between occasions 1 and 2.

#  new animals in occ2 = B1 = N*f
nx=rep(0,2^nper); names(nx)=nxnames;    # re-initialize cap-hist freq
B1=N*f                                  # set no. new individuals = N * f
cap(ch0, B1, 2); nx2=nx;                # generate histories for new animals (starting @ occ 2)

Next, we’ll do the same for occasions, 3-6.

#  new animals in occ3 = B2 = (N*S + B1)*f
nx=rep(0,2^nper); names(nx)=nxnames; 
B2= (N*S + B1)*f
cap(ch0, B2, 3); nx3=nx

# new animals in occ4 = B3 = (N*S*S  + B1*S   + B2) * f
nx=rep(0,2^nper); names(nx)=nxnames; 
B3= (N*S^2 + B1*S+ B2)*f
cap(ch0, B3, 4); nx4=nx

# new animals in occ5 = (N*S*S    + B1*S   + B2) * f
nx=rep(0,2^nper); names(nx)=nxnames
B4= (N*S^3 + B1*S^2 + B2*S + B3)*f
cap(ch0, B4, 5); nx5=nx

# new animals in occ5 = (N*S*S    + B1*S   + B2) * f
nx=rep(0,2^nper); names(nx)=nxnames
B5=(N*S^4 + B1*S^3 + B2*S^2 + B3*S + B4)*f
cap(ch0, B5, 6) ; nx6=nx

Now, we can add the cap-hist frequencies from all cohorts and print a table of the number of animals alive at each occasion (col) from each cohort (row).

nx=nx1+nx2+nx3+nx4+nx5+nx6;   # overall capture-history frequencies
ch=names(nx);

x=rbind(                       # table of cohort survival
  N*S^(0:5),                #  number alive from initial cohort at each time
  c(0,B1*S^(0:4)),          #  number alive from 2nd cohort at each time
  c(0,0,B2*S^(0:3)),        #  number alive from 3rd cohort at each time
  c(0,0,0,B3*S^(0:2)),      #  number alive from 4th cohort at each time
  c(0,0,0,0,B4*S^(0:1)),    #  number alive from 5th cohort at each time
  c(0,0,0,0,0,B5)              #  number alive from 6th cohort at each time
)
x=rbind(x,colSums(x))
rownames(x)=c('N',paste0('B',1:5),'totN')
cat('cohort survival:\n'); print(x)
cohort survival:
     [,1] [,2] [,3] [,4]  [,5]   [,6]
N     100   60   36 21.6 12.96  7.776
B1      0   20   12  7.2  4.32  2.592
B2      0    0   16  9.6  5.76  3.456
B3      0    0    0 12.8  7.68  4.608
B4      0    0    0  0.0 10.24  6.144
B5      0    0    0  0.0  0.00  8.192
totN  100   80   64 51.2 40.96 32.768

Looking at the table…

  • From the initial cohort of 100 animals, 60 survived to occasion 2, 36 survived to occasion 3, 21.6 survived to occasion 4.
  • From the 2nd cohort of 20 new animals in occasion 2, 12 survived to occasion 3, 7.2 survived to 4.
  • From the 3rd cohort of 16 new animals in occasion 3, 9.6 survived to occasion 4, 5.76 survived to 5.

Just to check, let’s compute recruitment from the table. Recruitment is the number of new individuals at time \(i\) per old animal at time \(i-1\). The number of new animals is the diagonal of this matrix, and the number of old animals is the column sums.

cat('true f:',diag(x)[-1]/x[7,-6],'\n ')
true f: 0.2 0.2 0.2 0.2 0.2 
 

Run Mark with all data and all intervals = 1

If we run Mark with these data, the estimates should match the input parameters.

library(RMark)
This is RMark 3.0.0
 Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
i=grep('1',ch); ch=ch[i]; nx=nx[i] # remove '00000' histories

pd=process.data(
  data=data.frame(ch=ch, freq=nx, stringsAsFactors=F),
  model="Pradrec", time.intervals=rep(1,nper-1))
dd=make.design.data(pd)
m1=mark(pd,dd,model.parameters=list(
  Phi=list(formula=~time), f=list(formula=~time),
  p=list(formula=~1)), silent=TRUE
)

Output summary for Pradrec model
Name : Phi(~time)p(~1)f(~time) 

Npar :  11
-2lnL:  565.7238
AICc :  589.6711

Beta
                     estimate        se        lcl       ucl
Phi:(Intercept)  4.054663e-01 0.5954453 -0.7616065 1.5725391
Phi:time2       -2.487561e-06 0.9138028 -1.7910561 1.7910511
Phi:time3        4.033106e-08 0.8427574 -1.6518045 1.6518046
Phi:time4        6.537285e-08 0.9166545 -1.7966429 1.7966430
Phi:time5       -4.073389e-06 1.0670394 -2.0914013 2.0913932
p:(Intercept)   -4.054653e-01 0.3159397 -1.0247071 0.2137765
f:(Intercept)   -1.609439e+00 0.9267610 -3.4258903 0.2070129
f:time2          1.176884e-06 1.3458317 -2.6378290 2.6378313
f:time3         -3.673460e-07 1.1851781 -2.3229494 2.3229487
f:time4          7.137739e-07 1.2273062 -2.4055196 2.4055210
f:time5          4.603908e-07 1.2848658 -2.5183366 2.5183375


Real Parameter Phi
         1         2         3         4         5
 0.6000003 0.5999997 0.6000003 0.6000003 0.5999993


Real Parameter p
   1   2   3   4   5   6
 0.4 0.4 0.4 0.4 0.4 0.4


Real Parameter f
         1         2         3   4         5
 0.1999998 0.2000001 0.1999998 0.2 0.1999999

Parameter estimates of recruitment f match the input of 0.2.

Try with data omitting occasion 4

To make the time-intervals different (one of them at least), we’ll pretend we didn’t sample at occasion 4. To do this, we simply remove the 4th column from each capture-history.

###############  remove occasion 4 *************************
ch=paste0(substr(ch,1,3),substr(ch,5,6)) 
#############################

i=grep('1',ch); ch=ch[i]; nx=nx[i] # remove '00000' histories

Now, run Mark with interval lengths correctly specified and check recruitment estimates.

pd=process.data(
  data=data.frame(ch=ch, freq=nx, stringsAsFactors=F),
  model="Pradrec", time.intervals=c(1,1,2,1))
dd=make.design.data(pd)
m2=mark(pd,dd,model.parameters=list(
  Phi=list(formula=~time), f=list(formula=~time),
  p=list(formula=~1)), silent=TRUE
)

Output summary for Pradrec model
Name : Phi(~time)p(~1)f(~time) 

Npar :  9
-2lnL:  446.3803
AICc :  465.9176

Beta
                     estimate        se        lcl       ucl
Phi:(Intercept)  4.054649e-01 0.6727744 -0.9131729 1.7241027
Phi:time2        2.092065e-07 0.9759202 -1.9128035 1.9128039
Phi:time3        2.925480e-07 0.7802033 -1.5291982 1.5291988
Phi:time5       -2.457663e-07 1.1161124 -2.1875806 2.1875801
p:(Intercept)   -4.054650e-01 0.4246352 -1.2377500 0.4268199
f:(Intercept)   -1.609437e+00 1.0158644 -3.6005317 0.3816570
f:time2         -4.089511e-07 1.3657436 -2.6768580 2.6768571
f:time3         -2.760293e-07 1.1109656 -2.1774930 2.1774924
f:time5         -2.348393e-06 1.3767900 -2.6985109 2.6985062


Real Parameter Phi
         1   2   3         5
 0.5999999 0.6 0.6 0.5999999


Real Parameter p
   1   2   3   5   6
 0.4 0.4 0.4 0.4 0.4


Real Parameter f
         1   2         3         5
 0.2000001 0.2 0.2000001 0.1999996

With time-intervals correctly specified, notice the estimates of phi are all 0.6, matching the input parameter. Also, all estimates of recruitment are equal 0.2, matching input. This tells us that Mark is producing estimates of phi and f which are scaled to the time-interval lengths (ie., survival rate is survival per unit time, not survival per interval, and recruitment is per unit time, not per interval.)

Since we omitted occasion 4, the estimate of f(3) is the estimated number of new animals in occasion 5 per old animal in occasion 4, which is assumed to be equal to the number of new animals in occasion 4 per old animal in occasion 3.

Now, run Mark with all time-intervals = 1.

##  re-run with time.intervals all = 1
pd=process.data(
  data=data.frame(ch=ch, freq=nx, stringsAsFactors=F),
  model="Pradrec", time.intervals=rep(1,nper-2))
dd=make.design.data(pd)
m3=mark(pd,dd,model.parameters=list(
  Phi=list(formula=~time), f=list(formula=~time),
  p=list(formula=~1)), silent=TRUE
)

Output summary for Pradrec model
Name : Phi(~time)p(~1)f(~time) 

Npar :  9
-2lnL:  446.3803
AICc :  465.9176

Beta
                     estimate        se        lcl       ucl
Phi:(Intercept)  4.054655e-01 0.6727736 -0.9131707 1.7241017
Phi:time2       -5.779834e-07 0.9759185 -1.9128009 1.9127997
Phi:time3       -9.808297e-01 0.8423538 -2.6318432 0.6701837
Phi:time4       -7.685012e-07 1.1161089 -2.1875742 2.1875726
p:(Intercept)   -4.054651e-01 0.4246354 -1.2377505 0.4268203
f:(Intercept)   -1.609438e+00 1.0158666 -3.6005369 0.3816602
f:time2          5.083897e-07 1.3657423 -2.6768545 2.6768555
f:time3          3.364726e-01 1.1262982 -1.8710719 2.5440172
f:time4          9.013057e-07 1.3767946 -2.6985166 2.6985184


Real Parameter Phi
         1   2    3         4
 0.6000001 0.6 0.36 0.5999999


Real Parameter p
   1   2   3   4   5
 0.4 0.4 0.4 0.4 0.4


Real Parameter f
         1   2    3         4
 0.1999999 0.2 0.28 0.2000001

Here, the survival and recruitment estimates are per interval, so the 3rd interval survival is 0.36, or the square of 0.6 since the interval length is 2. The 3rd estimate of recruitment is 0.28 meaning that the number of new animals in occasion 5 per animal alive in occasion 3 is higher than the other recruitment estimates due to the longer interval length.