Removal models accounting for temporary emigration

Summary Removal of protected species from sites scheduled for development is often a legal requirement in order to minimize the loss of biodiversity. The assumption of closure in the classic removal model will be violated if individuals become temporarily undetectable, a phenomenon commonly exhibited by reptiles and amphibians. Temporary emigration can be modeled using a multievent framework with a partial hidden process, where the underlying state process describes the movement pattern of animals between the survey area and an area outside of the study. We present a multievent removal model within a robust design framework which allows for individuals becoming temporarily unavailable for detection. We demonstrate how to investigate parameter redundancy in the model. Results suggest the use of the robust design and certain forms of constraints overcome issues of parameter redundancy. We show which combinations of parameters are estimable when the robust design reduces to a single secondary capture occasion within each primary sampling period. Additionally, we explore the benefit of the robust design on the precision of parameters using simulation. We demonstrate that the use of the robust design is highly recommended when sampling removal data. We apply our model to removal data of common lizards, Zootoca vivipara, and for this application precision of parameter estimates is further improved using an integrated model.


Web Appendix A: Parameter redundancy
A model is called parameter redundant when parameters in the model are confounded, so that the model could be reparameterized in terms of a small number of parameters (Catchpole and Morgan, 1997). In Section 3 of the main paper we discussed parameter redundancy of different removal models. Here in Section 1.1 we provide more detail on how to detect parameter redundancy for a fixed number of sampling occasions, giving an illustrative example in Section 1.1.1. Then in Section 1.2 we show how to generalize parameter redundancy results for any number of sampling occasions. Theoretically we can estimate all the parameters in a model that is not parameter redundant. However, good performance of the model is not guaranteed if the model is in same way close to a model that is parameter redundant; this is known as near redundancy. Near redundancy is discussed in Section 1.3. 1 1.1 Methods for detecting parameter redundancy results We have briefly described how to determine whether or not the model is parameter redundant in Section 3 of the main paper. This involves forming a derivative matrix, where θ is a vector of h parameters and where κ(θ) is a vector of parameter combinations that uniquely represent the model. This is known as an exhaustive summary (Cole et al., 2010). For removal models a suitable exhaustive summary is a vector of probabilities of individuals being removed at the kth occasion, for k = 1, . . . , K. This is similar to the exhaustive summary used in Cole et al. (2012).
The deficiency of the model, d, is calculated as h − r, where r is the rank of the derivative matrix. If d > 0, the model is parameter redundant. Otherwise, the model is termed full rank and is not parameter redundant. For a parameter redundant model, it is possible to show, if any, of the original parameters are estimable by solving α T D = 0, for a vector α(θ). There will be d solution to α T D = 0. If we denote entries of the d solutions as α i,j , for i = 1, . . . , h and j = 1, . . . , d, then any zero α i,j for all j corresponds to an estimable parameter. To find other estimable parameter combinations, we solve the system of linear first-order partial differential equations, v i=1 α i,j ∂f /∂θ i = 0, where j = 1, . . . , d (Catchpole et al., 1998;Cole et al., 2010). We demonstrate how this method works using an illustrative example in Section 1.1.1 below. Note that all of the parameter redundancy results in Tables 1 and 2 in the paper are based on the assumption that at least one animal is removed at every sampling occasion. This ensures that all of the exhaustive summary terms are present as we are interested in detecting parameter redundancy in the model itself instead of parameter redundancy due to the lack of data.

NNC example
The NNC model in Table 1 in the main paper considers a multievent framework allowing for temporary emigration but does not have a robust design structure. It contains a constant detection probability, constant transition probabilities and an initial state parameter. Therefore, the parameters in this case are θ = [p φ 12 φ 21 π], where p is the detection probability, φ 12 and φ 21 represents transition probabilities between state 1 and state 2, and π is the probability of being in the state 1 initially. Consider the model with 8 occasions of removal, and we define the exhaustive summary as a vector of probabilities of an individual being captured at the the kth occasion as the specification determines the model. We only show first three terms in the probability matrix below as the probabilities of removal dramatically get more complicated as the number of sampling times increase. All the Maple files are available in the supplementary material.
where L k,1 is the probability of animals being captured at the kth sampling occasion. We remove the more complicated probability that an animal is never captured to make the exhaustive summary simpler, as we assume we observe at least one individual for all possible other capture histories and all the probabilities sum up to one (Catchpole and Morgan, 1997). The derivative matrix is given by, and is found to have symbolic rank 3. However, there are 4 parameters, therefore the model is 2 parameter redundant with deficiency 1. In order to find if any of the original parameter can be estimated, we solve α T D = 0, which gives, As there are no zeros in α, none of the original parameters are individually estimable. To find the estimable parameter combinations, we solve the partial differential equation, The solution of which shows we can estimate pφ 21 , πp, and p(φ 12 − 1) − φ 12 − φ 21 .

Generalization of parameter redundancy results
In Section 1.1.1, we demonstrate how to determine the results of parameter redundancy for a fixed number of K sampling occasions with two secondary samples within each primary period. It is possible to generalize the results of parameter redundancy to any number of samples as described in Catchpole and Morgan (1997) and Cole et al. (2010). In this section, we show how to obtain general parameter redundancy results for our removal models with two secondary occasions within each primary period. Note that we could theoretically generalize the results to any number of secondary occasions. Because in this paper the data collected has two secondary occasions within each primary period, so we demonstrate the generalization of parameter redundancy results for data with two secondary occasions.
Furthermore, if there is only one extra parameter in the extended model, then D 2,2 (θ 2 ) is always full rank because it always has rank of one as a row vector. Therefore the extended model will be full rank if D 1,1 (θ 1 ) is full rank and there is only one new parameter in θ 2 . Therefore, the model is always full rank by induction. The use of extension theorem is demonstrated in the Maple file for the R-SR t C model in Table 2 in the paper.

Reparameterization theorem
If the original model is not full rank, then the extension theorem cannot be used for that model. Instead, we need to find a reparameterization of the original model which is full rank to begin with, then we can apply the extension theorem to the reparameterized model to determine its the general parameter redundancy result (Cole et al., 2010). We reparameterize our models in terms of the estimable parameter combinations. The deficiency and estimable parameter combinations are derived explicitly in Web Table 1 which provides additional detail on the estimable parameter combination for the parameter redundant models in the Table 2 in the main paper. We demonstrate how to use the reparameterization theorem for the model R-NN t C in Table 2 in the paper in the Maple file.

Near parameter redundancy
We have described near parameter redundancy for a full rank model in Section 3 of the paper. Following the method of detecting near redundancy in Catchpole et al. (2001), we calculated the expected information matrix in Maple and obtained the smallest standardized eigenvalues for the RMER models with constant capture probability in the paper, the standardized eigenvalues are computed as the absolute value of the ratio of the smallest eigenvalue and the largest eigenvalue. A very small eigenvalue in a full rank model indicates a problem with near redundancy. The results are showed in Web Tables 2-5, where we described Scenarios 1-2 and Scenarios 3-8 in the main paper and in Section 2 in the Web Appendix B respectively.
R-NNC has extremely small standardized eigenvalues across all scenarios in Web Table 2, which results in biased parameter estimates in the simulation studies given in Figures Table 3 we show that the standardized smallest eigenvalues for the R-SEC model are larger than those in the R-NEC model. Also, We observe similar results in the simulation study in Section 2.1.2, where the performance of R-SEC is better than R-NEC under each Scenario shown in Web Figure 2.
In Web Table 4, the R-SR t 2 C and R-SR t C have largest two standardized eigenvalues under each scenario, which suggests these two models behave better than the rest of the models. It matches what we found in Web Figures 5 and 6 in Section 2.2.1. Also, relatively smaller standardized eigenvalues for R-NR t C model in Web Table 4 explain its poor behaviours in Web Figures 5-6. Although we observe relatively large standardized eigenvalues for models in Web Table 5 especially under Scenarios 6 and 8, none of these models yield unbiased results for time-dependant transition probabilities in Web Figures 8.
We also notice there is no clear cut-off point between large and small standardized eigenvalues across all Scenarios in our cases. So, we suggest that any interpretations of standardized eigenvalues have to be made cautiously.
Web Table 2: Results of standardized smallest eigenvalues (S.Eig) for R-NNC, R-SNC, R-NRC and R-SRC models in Table 1 in the paper under Scenarios 1-4. All models are diagnosed with K = 8 sampling occasions, where there are T = 4 primary periods and k 1 = k 2 = k 3 = k 4 = 2 secondary samples for each primary period.

Scenario
Model Web Table 5: Results of standardised smallest eigenvalues (S.Eig) for R-NEC, R-SEC and R-NE 2 C models under Scenarios 5-8. All models are diagnosed with K = 8 sampling occasions, where there are T = 4 primary periods and k 1 = k 2 = k 3 = k 4 = 2 secondary samples for each primary period.
In this section, we provide more simulation results for our proposed models. Models with constraint "R" or "R t " are investigated with low capture probability under Scenarios 1-2 as described in the paper, while simulations are conducted for high capture probability under Scenarios 3-4. In addition, simulations are carried out for constraint "E" or "E t " under Scenarios 5-8 as shown below. The true values of capture probabilities p and transition probabilities φ 12 i for simulating the data under Scenarios 5-8 for the constraint "E" or "E t " are same with those under Scenarios 1-4 described in the paper. Web Table 6 summarizes the true values of parameters we used for simulations for all scenarios we considered.
In addition to investigating the performance of the MER, RMER and IRMER models using simulation, we also run simulations for geometric removal models under each simulation settings. In Sections 2.1.3, 2.2.3 and 2.3.2 we show the results obtained from simulations for geometric removal models using data simulated with robust design under simulation Setting 4.1 Setting 4.2 and Setting 4.3 respectively. We observe extremely large estimates of population sizes and very tiny estimates of constant capture probabilities compare to the true values of parameters across all Scenarios. Therefore, we show that the geometric removal model overestimates the population size and underestimates the constant capture probability if temporary emigration and the robust design are ignored.

Results under Setting 4.1
We have discussed the simulation results in Figure 1-5 together with the results of near parameter redundant in Section 1.3. The best performing model is the R-SRC model.

Results under Scenarios 1-4
We have shown the results for the estimates of population size N , capture probability p and transition probability φ 12 under Scenarios 1 and 2 in Figure Table 6: True values of parameters used for different simulation setting. Scenarios 1-4 are set for constraint "R". Scenarios 5-8 are set for constraint "E".  Figure 1: Results of estimated population size N (A), capture probability p (B) and transition probability φ 12 (C) for simulations with K = 10 and K = 20 sampling occasions displayed on the first and second columns respectively. The black horizontal lines are the true values used for simulating the data under different scenarios. These models are presented in the Table 1 in the paper.

Results under Scenarios 5-8
In this Section, we show simulation results for models with constraint "E" in Table 1 in the main paper under Scenarios 5-8. When the population has high mobility under Scenarios 5 and 7, individuals are exposed to the study area more frequently, so nearly all of the individuals are removed from the population by the end of the study. Therefore, we observe boundary estimates of population sizes as shown in Web Figure 2(A). In Web Figure 2(B) the results of estimated capture probabilities are unbiased apart from those obtained from the NEC model. In Web Figure  2(C), unbiased estimates of φ 12 are only obtained from the R-SEC model when we have K = 20 sampling occasions. However, the results in Section 2.2.2 suggest that biased estimated transition probabilities are obtained from the models incorporated with the constraint "E t " when timedependent transition probabilities are considered due to near parameter redundancy.  Figure 2: Results of estimated population size N (A) and capture probability p (B) and transition probability φ12 (C) for simulations with K = 10 and K = 20 sampling occasions displayed on the first and second columns respectively. These models are the corresponding models in Table 1 in the paper. The black horizontal lines are the true values for simulating the data under different scenarios.
11 2.1.3 Results for the GRM model We described the simulation Setting 4.2 in the main paper. We have discussed the simulation results in Figure 5-8 together with the results of near parameter redundancy in Section 3. The best performing model for time-dependent transition probabilities are the models with at least the constraint "S" and "R t ", i.e. R-SR t C and R-SR t 2 C models. None of the rest of model yield unbiased estimates for all of the parameters due to near parameter redundancy as discussed in Section 1.3.

Results under Scenarios 1-4
The results for the estimates of population size N , capture probability p and transition probability φ 12 under Scenarios 1 and 2 are illustrated in Figure 2 in the main paper. The results for the estimated population sizes N and capture probabilities p for the same models under Scenarios 3 and 4 are shown in Web Figure 5. In Web Figure 6 we show the estimates of time-varying transition probabilities for Scenarios 1-4.
Web Figure 5 illustrates that unbiased estimates of population sizes are obtained from both the R-SR t C and R-SR t 2 C models when there are 20 sampling occasions. However, the estimated population sizes are slightly biased low for R-SR t C for 10 sampling occasions. In addition, R-NR t C underestimates the population sizes due to near parameter redundancy as shown in Section 1.3. In Web Figure 6 we observe that the near parameter redundant R-NR t C model underestimates most of the transition probabilities. Furthermore, the estimated transition probabilities φ 12 between the last two primary periods obtained from the R-NR t C model are on the boundary across all Scenarios.
Therefore we conclude that when modelling time-dependant transition probabilities, we need to consider at least the combination of the constraint "S" and "R t ", otherwise the model is either parameter redundant or near parameter redundant.  Figure 5: Results of estimated population sizs N (A) and capture probability p (B) for simulations with K = 10 and K = 20 sampling occasions displayed on the first and second columns respectively. These are the corresponding models in Table 2 in the paper. The black horizontal lines are the true values used for simulating the data under Scenarios 3 and 4.

Results under Scenarios 5-8
In this Section, we present simulation results for the full rank models with constraint "E t " in Table  2 in the main paper, i.e. R-NE t C, R-SE t C and R-SE t 2 C. The boundary estimates of population sizes shown in Web Figure 7(A) is because we almost capture all of the individuals in the population during the study of length K = 20. The results in Web Figure 8 suggest that time-dependent transition probabilities obtained from the models incorporated with the constraint "E t " are reliable due to near parameter redundancy as discussed in Section 1.3. Therefore, the use of constraint "E t " is not recommended if we want to model time-dependent transition probabilities.   In this section, we present the simulation results under Scenarios 1-4 for the IRMER models in Table 2 in the paper. The detail for the simulation Setting 4.3 is given in the main paper.

Results under Scenarios 1-4
We have demonstrated that we need to consider at least the combination of the constraint "R t " and "S" in the model to get unbiased estimates of parameter in Section 2.2 as well as in the paper.
In this section, we show that the integrated modelling approach (Besbeas et al., 2002) cannot improve the biased estimates in the IR-NR t C model with the constraint "R t " only as shown in Web Figure 11-13. In Web Figures 11 and 13, the estimated population sizes N , M and transition probabilities are biased low for the IR-NR t C model across all Scenarios. Therefore, in the IRMER models the incorporation of both the constraint "R t " and "S" is still needed.   IR−NR t C IR−S 1,2 R t C IR−S 1,2 R a t C IR−S 1 R a t C IR−NR t C IR−S 1,2 R t C IR−S 1,2 R a t C IR−S 1 R a t C IR−NR t C IR−S 1,2 R t C IR−S 1,2 R a t C IR−S 1 R a t C Web Figure 12: Results of estimated capture probability p for simulations with K = 10 and K = 20 sampling occasions displayed on the first four and last four columns respectively. The black horizontal lines are the true values used for simulating the data under different scenarios. IR−NR t C IR−S 1,2 R t C IR−S 1,2 R a t C IR−S 1 R a t C IR−NR t C IR−S 1,2 R t C IR−S 1,2 R a t C IR−S 1 R a t C IR−NR t C IR−S 1,2 R t C IR−S 1,2 R a t C IR−S 1 R a t C Web Figure 13: 3 Web Appendix C: Investigation of the relaxation of constraint "R t Here we investigate the generalization of the constraint "R t " (φ 12 i +φ 21 i = 1). We assume φ 12 i +φ 21 i = v, where v ∈ [0, 2] (denoted as "V t " in Web Table 6). This constraint is motivated by the constraint of φ 12 i +φ 21 i = 1, which relaxes the assumption of the sum of two transition probabilities at primary period i being equal to 1. The free parameters are fully time-dependent φ 12 i and v which is constant over time. Note that the limits of possible values of φ 12 i are conditional on v. The upper and lower bounds for φ 12 i are shown in the following, Parameter redundancy results are in Web Table 6 If we regard v as a fixed number over time (not a free parameter), the parameter redundancy results are the same with those in Table 1 and 2 presented in the paper. Under this specification together with constraint "S", unbiased estimates can be obtained as shown in the simulation study for v = 1 in the paper. Web Table 7: Results of deficiency of RMER models, where φ 12 i are time-varying over time. All models are diagnosed with an even number of sampling occasions, where there are two secondary samples within each primary period. h is the number of parameters in the model, r is the rank of the derivative matrix, d is the deficiency of the model, calculated as h-r. "t" denotes time-dependence, "c" denotes a constant parameter, and "t+cov" denotes additive effect in terms of covariates. PRS represents the parameter redundancy status, where FR indicates that all parameters in the model are theoretically estimable, PR indicates that the model is parameter redundant, NR indicates that the model is full rank but near parameter-redundant for the scenarios we considered. EP represents the estimable combinations of parameters if the model is parameter redundant. All results hold for K ≥ 4 occasions.
Model code Models h r d PRS EP R-NV t C π, φ 12 (t), p(c) K K/2+1 K/2-1 PR p, π, f (p, π, φ 12 1 , φ 21 1 , v), · · · , f (p, π, φ 12 1 , · · · , φ 12 In this section we describe the non-parametric bootstrap procedure used for analysis of the common lizard data presented in Section 5 of the paper. Bootstrap samples are obtained by sampling with replacement from the original data. Bootstrap estimates are denoted by * on the superscript of the parameters, e.g. p * denotes the bootstrap estimates of the constant capture probability. The algorithm is shown below. • Step 1: rewrite the removal data in terms of individual removal records for all sampling occasions. If we removed n k individuals at the kth sampling occasion, we rewrite this with n k "k"s individual records. For example, the first three observations of the juvenile data set are 0, 5, 2. We rewrite them as five "2"s and two "3"s, i.e. 2, 2, 2, 2, 2, 3, 3 and repeat the same process for all of the observations. • Step 2: sample D individual removal records with replacement and rewrite them in terms of removal count data by calculating the resulting frequency table, where D is equal to the total number of individuals we observed in the original data. • Step 3: repeat Step 2 a large number of times (e.g. 500 resamples) to produce bootstrap samples. • Step 4: fit the model to each bootstrap sample and obtain the bootstrap estimates. • Step 5: Obtain the standard errors by calculating the standard deviation of the distribution of bootstrap estimates. The (1 − a)% confidence intervals are constructed by finding the (a/2%, 1 − a/2%) percentiles of the distribution of the bootstrap estimates. For example, the standard errors of p is calculated by the standard deviation of the distribution of p * and the 95% confidence interval is computed by (p * 0.025 , p * 0.975 ) where p * 0.025 and p * 0.975 are the 2.5% 97.5% percentiles of the bootstrapped p * respectively.
In Web Figure 18 we show 500 bootstrap samples obtained for the juvenile and adult common lizard data sets analyzed in Section 5 of the paper. The model is fitted to each of the bootstrap samples and parameter estimates are computed that are subsequently used to construct confidence intervals and calculate standard errors. In this section we present more results for the real data analyzed in Section 5 in the paper. In Web Table 8 We show all of the estimates apart from time-varying transition probabilities with standard errors and confidence interval calculated empirically from 500 bootstrap samples for the top three models in the Table 3 in the main paper. The fitted counts of individuals removed at each occasion for both juvenile and adult population are displayed in Web Figure 19. Furthermore, the estimates of time-dependent transition probabilities φ 12 i are shown in Web Figure 20. As we described in Section 3. the constraint "V t " (φ 12 where v ∈ [0, 2]) only works when individuals have high mobility. Moreover, if we treat v as a fixed constant (denoted as "V c ") rather than a free parameter. The parameter redundancy results are the same with those with the constraint "R t ". We fit the model IR-S ju,ad V c a Z to the real data over a grid fixed values of v. We show that the assumption of φ 12 i + φ 21 i = 1 results in the best fit for the common lizards data in Web Figure 21. Therefore we stick to the constraint "R" (φ 12 i + φ 21 i = 1) in our candidate models for the real data.
Web Table 8: Estimates and their standard errors (SE) and 95% confidence intervals (95% CI) for the three models with the lowest AIC in Table 3 in the paper for the common lizards data.n 0,ju and n 0,ad are the estimates of the number of juvenile and adult common lizards we failed to capture by the end of the study respectively.α andβ are the estimated logistic regression parameters (intercept and slope) for capture probability respectively.γ ad is the estimate of additive effect for adult common lizards on transition probability.π ad is the estimate of initial state parameter for the adult population.  i + φ 21 i = v, v ∈ (0, 2)) is a constant over time for common lizards data.
Web Table 9: Results of standardized smallest eigenvalues (S.Eig) and median relative bias in estimatesp andψ under different Scenarios, where true parameters are listed in p true and ψ true . All models are diagnosed with K = 8 sampling occasions where there are T = 4 primary periods and two secondary samples within each primary period. * indicates all individuals are assume to survive during the study.  Table 10: Median relative bias (MBias) in estimatesp,ψ andN from the R-SR t d C model under different Scenarios. Simulations are conducted for simulated N = 500 individuals with K = 20 sampling occasions where there two secondary samples within each primary period. True values of true parameters of p true and ψ true are listed below. True transition probabilities used are the same (0.8,0.7,0.8,0.3,0.6,0.7,0.8,0.6,0.6) for every Scenario.
This model R-SR t d C is determined to be full rank. We also investigate its eigenvalues under various Scenarios and we list the true values of parameters we used in Web Table 9. When capture probability is low its standardized smallest eigenvalue under Scenario A1, A2 or A3 is worse than those both R-NR t C and R-SR t C models under Scenario 1 as shown in Web Table 9. Similar conclusion is obtained when capture probability is high. Therefore it is not reliable to use as their standardized smallest eigenvalues are smaller than those from the R-NR t C model which is a near parameter redundant model.
We also conducted simulation study for a population of N = 500 individuals with K = 20 sampling occasions. We list the true parameters in Web Table 10 where φ 12 i are (0.8,0.7,0.8,0.3,0.6,0.7,0.8,0.6,0.6) for each of the Scenarios considered, i.e. individuals are assume to tend to stay off-site. The estimated population sizes N are negatively biased across all scenarios in Web Table 10 and the estimates of survival probabilityψ are over estimated apart from those under Scenario A1. The simulation results in Web Table 10 suggests the similar conclusion in Web Table 9, i.e. the R-SR t d C model with a constant survival probability cannot be used in practice.