Barring a really strange coincidence, I suspect that wijw are working on the exact same problem, which I have continued to pursue without being able to get to the bottom of it. Having prepared the following description of the issue, I am going to go ahead an add it here in hopes that some further details will help:
I am studying recapture observations of released cohorts in a reintroduction program. I am not considering any natural recruitment. I am using the POPAN JS model and I am trying to fix the known number released for each cohort by fixing b[t] (i.e. pent) parameters and N (the superpopulation). However I am getting unexpected results: if the known number of releases just before occasion t is called R[t-1] then I expect Nhat[t] would be exactly R[t-1], which would be exactly B[t-1] which would be exactly N. Instead these quantities are B[t-1] = Nhat[t] = N + the number of encounter histories.
I have simulated some encounter histories, which are pasted at the end, for a simpler 1 cohort (= 1 release occasion) example. I simulated a set of 5-occasion encounter histories for 30 individuals released between the 1st and the 2nd occasion with constant p = 0.5 and constant phi = 0.7. After removing all-zero encounter histories there were 24 remaining. I set up a the PIM matrices for 5 parameters: (parm1): p(.), (parm2): phi(.), (parm3):b[1] = 1 (fixed), (parm4):b[2] = b[3] = b[4] = 0 (fixed), and (parm5): N = 30 (fixed). This gives incorrect results and two warnings:
** WARNING ** Error number 0 from VA09AD optimization routine.
** WARNING ** Invalid calculation occurred during processing derived parameters of this model.
The output shows that N = 30 but Nhat[2] =Bhat[1] = 54 (i.e. 30 + 24 encounter histories).
- Code: Select all
Parameter Estimate Standard Error Lower Upper
-------------------------- -------------- -------------- -------------- --------------
1:Phi 0.6832247 0.0950899 0.4768883 0.8361389
2:p 0.2571128 0.0509876 0.1701962 0.3686945
3:pent 1.0000000 0.0000000 1.0000000 1.0000000 Fixed
4:pent 0.0000000 0.0000000 0.0000000 0.0000000 Fixed
5:N 30.000000 0.0000000 30.000000 30.000000 Fixed
Grp. Occ. B-hat Standard Error Lower Upper
---- ---- -------------- -------------- -------------- --------------
1 1 54.000000 0.0000000 54.000000 54.000000
1 2 0.0000000 0.0000000 0.0000000 0.0000000
1 3 0.0000000 0.0000000 0.0000000 0.0000000
1 4 0.0000000 0.0000000 0.0000000 0.0000000
Grp. Occ. N-hat Standard Error Lower Upper
---- ---- -------------- -------------- -------------- --------------
1 1 0.0000000 0.0000000 0.0000000 0.0000000
1 2 54.000000 0.0000000 54.000000 54.000000
1 3 36.894136 5.1348557 28.122526 48.401673
1 4 25.206987 7.0165202 14.756587 43.058206
1 5 17.222037 7.1907901 7.8496743 37.784825
If I re-run the model with N fixed to 6 (i.e. the number of released but never-encountered individuals) the Nhat and Bhat values are correct and phi and p are close to the simulation values. However, this kludge does not make sense to me based on my understanding of what N represents. I understand that there is a b[0] value that is not shown and that is obtained as 1 - sum(b[1] ... b[n.occ-1]) but I don't see how that explains the result because the fixed values should ensure that sum(b[1] ... b[n.occ-1]) = 1 and therefore b[0] should be 0.
Possibly relevant: when I only fix parm 4 (i.e. b[2 ... 4]) to 0 the resulting estimate of b[1] is 1, which matches my understanding of the multinomial logit link BUT if I switch it around and fix b[1] = 1 and let b[2 ...4] be estimated I get all pent approximately = 1 (how is that possible if they share a mlogit(1) link?) negative abundances at early occasions, and 30 individuals "being born" every occasion.
Please help; what am I missing?
encounter histories to reproduce the analysis:
- Code: Select all
00001 1;
00010 1;
00011 1;
00100 5;
00101 3;
00111 1;
01000 9;
01010 1;
01100 1;
01101 1;