Robust Design (RDHuggins) weird population estimates

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

Robust Design (RDHuggins) weird population estimates

Postby wijw » Mon Jul 11, 2022 2:54 pm

I have a dataset with 483 unique individuals that are grouped into males, females, and unknown (juveniles). I conducted mark-recapture surveys from 2018 to 2021. My goal is to estimate abundance and survival probability.

I decided on the Robust Design model. I have four primary periods (2018 – 2021), and three secondary periods in each primary period (April-May, June-July, Aug-Sept). In 2020, I could not conduct surveys for the first secondary period (April-May). So, I created a dummy occasion to make sure the time intervals are equal. I also fixed the p and c for that occasion at 0.

Code: Select all
#Set time intervals between primary and secondary occasions
time.intervals = c(0,0,1,0,0,1,0,0,1,0,0)

#Create the processed dataframe and design data
rd = process.data(data = rd.txt,
                          model = "RDHuggins",
                          groups = "sex",
                          time.intervals = time.intervals)

ddl = make.design.data(rd)
 
#Fix the capture and recapture probability for 2020 THIRD PRIMARY AND FIRST SECONDARY occasion to zero (it is a dummy occasion).
ddl$p$fix=NA
ddl$p$fix[which(ddl$p$time==1 & ddl$p$session == 3)] = 0
ddl$c$fix=NA
ddl$c$fix[which(ddl$c$time==1 & ddl$c$session == 3)] = 0

In terms of survival estimates, I considered several models (constant/transience/sex/additive and interactive models). One thing to note is that in 2020, the population experienced a mass mortality event, and I found many unmarked dead female individuals. (Because they are unmarked I cannot add them to my analysis.)

Code: Select all
#Include the mass mortality year
ddl$S$mme = 0
ddl$S$mme[ddl$S$time == 2] = 1

I also found no evidence of trap dependence based on U-CARE. So my p and c estimates were shared.
Code: Select all
p.effort = list(formula = ~effort, share = TRUE)
p.season = list(formula = ~season, share = TRUE)
p.constant = list(formula = ~1, share = TRUE)

I also considered markovian, random, and no emigration.

When I ran different models, the best model, S(~1)Gamma''(~1)Gamma'(~1)p(~effort)c() and models that were not very well supported, gave me weird abundance estimates for 2020. There was an upward spike for females, and a downward spike for males and unknowns (basically contradicting what actually happened.) However, when I ran the code again without fixing the p and c for 2020 first secondary occasion, the results appear to resemble what occurred in the field. So I think the obtained estimates are likely a result of my recapture probability and the number caught for each group.

M(t+1):
Group 1 = sexFemale
66 56 54 47
Group 2 = sexMale
69 54 26 65
Group 3 = sexUnknown
38 22 18 49

Code: Select all
Real Function Parameters of { S(~1)Gamma''(~1)Gamma'(~1)p(~effort)c() }
                                                               95% Confidence Interval
     Parameter                Estimate       Standard Error     Lower           Upper
 --------------------------  --------------  --------------  --------------  --------------
     1:S gFemale c1 a0 t1     0.6003368       0.0570538       0.4852039       0.7053569                         
     2:p gFemale s1 t1        0.1102262       0.0137007       0.0860983       0.1400791                         
     3:p gFemale s3 t2        0.0829668       0.0138253       0.0595882       0.1144017                         
     4:Gamma'' gFemale c1 a   0.0000000       0.0000000       0.0000000       0.0000000      Fixed



Code: Select all
Estimates of Derived Parameters
   Population Estimates of { S(~1)Gamma''(~1)Gamma'(~1)p(~effort)c() }
                                                95% Confidence Interval
 Grp. Sess.     N-hat        Standard Error      Lower           Upper
 ---- -----  --------------  --------------  --------------  --------------
   1     1    223.29864       33.712874       169.83246       304.29601   
   1     2    189.80800       29.862809       142.85885       262.13419   
   1     3    340.29669       68.999878       233.70275       510.11873   
   1     4    205.85249       41.044656       143.51581       308.45058   
   2     1    233.87057       34.988066       178.26660       317.77048   
   2     2    183.02914       29.066215       137.42636       253.55946   
   2     3    163.84655       39.422568       105.56571       264.81735   
   2     4    284.68961       53.419421       202.33903       416.41887   
   3     1    128.79828       22.573274       94.185223       184.73482   
   3     2    74.567427       15.678404       51.663858       115.15492   
   3     3    113.43223       30.494711       69.793866       193.83762   
   3     4    214.61217       42.429045       150.03482       320.46475


Any advice is appreciated.
wijw
 
Posts: 23
Joined: Mon Oct 26, 2020 2:27 pm

Re: Robust Design (RDHuggins) weird population estimates

Postby jlaake » Mon Oct 17, 2022 2:12 pm

Look at your design data for p and c. There should not be a c for first occasion in a session because it is a recapture probability. See what your design data looks after setting fix values.

This is what robust example looks like with RDHuggins

Code: Select all
> data(robust)
> time.intervals=c(0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0)
> pd=process.data(robust,model="RDHuggins",time.intervals=time.intervals)
> ddl=make.design.data(pd)
> ddl$p
   par.index model.index group time session Time c
1          1          27     1    1       1    0 0
2          2          28     1    2       1    1 0
3          3          29     1    1       2    0 0
4          4          30     1    2       2    1 0
5          5          31     1    1       3    0 0
6          6          32     1    2       3    1 0
7          7          33     1    3       3    2 0
8          8          34     1    4       3    3 0
9          9          35     1    1       4    0 0
10        10          36     1    2       4    1 0
11        11          37     1    3       4    2 0
12        12          38     1    4       4    3 0
13        13          39     1    5       4    4 0
14        14          40     1    1       5    0 0
15        15          41     1    2       5    1 0
> ddl$c
   par.index model.index group time session Time c
1          1          42     1    2       1    0 1
2          2          43     1    2       2    0 1
3          3          44     1    2       3    0 1
4          4          45     1    3       3    1 1
5          5          46     1    4       3    2 1
6          6          47     1    2       4    0 1
7          7          48     1    3       4    1 1
8          8          49     1    4       4    2 1
9          9          50     1    5       4    3 1
10        10          51     1    2       5    0 1
>


If you can't locate the problem, then I'll need your code and data to get any further.

I'm not entirely sure why you had to add the dummy occasion. Unless I'm wrong I believe RDHuggins is assuming no losses during the session (April-Sept).

--jeff
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Robust Design (RDHuggins) weird population estimates

Postby jlaake » Fri Oct 21, 2022 5:08 pm

I worked offlist to diagnose any problems with regard to this post. There were no problems with the code but a misunderstanding about the model and the user's anticipation of the results given their knowledge that 50 females were found dead.

They expected the abundance estimate of females to be lower in sessions 3 and 4 and when that didn't happen they thought there was a problem. But if you look at the number captured below you can see that the number of females caught didn't decline and because the model had a lower capture probability for sessions 3 and 4 than 1 and 2, it caused abundance to increase.

Code: Select all
M(t+1):
Group 1 = sexFemale
66 56 54 47
Group 2 = sexMale
69 54 26 65
Group 3 = sexUnknown
38 22 18 49


Also since there were a fair number of individuals with unknown sex, potentially they could be mostly one sex or another.

When they removed the fixing of p=0 for the "added" occasion then it incorrectly estimated p for the session because it computed a constant probability of capture across 3 occasions rather than 2. Thus it estimated p* as 1-(1-p)^3 for session 3 rather than 1-(1-p)^2 which had the effect of inflating p* and lowering abundance estimate incorrectly because there were only 2 sampled occasions within the session.

Note also that while they needed to add

Code: Select all
ddl$c$fix=NA


so the parameters could be shared, they did not need

Code: Select all
ddl$c$fix[which(ddl$c$time==1 & ddl$c$session == 3)] = 0


because there isn't a recapture probability for the first occasion within a session.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to RMark

Who is online

Users browsing this forum: No registered users and 9 guests