{'fraction_non_alphanumeric': 0.030546380193339193, 'fraction_numerical': 0.013485257129003909, 'mean_word_length': 4.702436495593572, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling latent individual heterogeneity in mark-recapture data with Dirichlet process priors

The natural subgroups often seen in mark-recapture studies and the complexity of real mark-recapture data means that parametric and discrete style models can be insufficient. Non-parametric models avoid these often restrictive assumptions. We consider the non-parametric Dirichlet process for modelling latent individual heterogeneity in probability of observation and the probability of remaining in or out of a marine sanctuary. Simulation studies demonstrated accurate estimation of multiple groups of latent individual heterogeneity. Simulations were also used to identify the limits of the Dirichlet process. The ability of the Dirichlet process to pick up unimodal heterogeneity was explored in order to avoid potential spurious multimodality. In application to a subset of the data from the North Atlantic humpback whales we were able to estimate annual population-level variation in usage of the marine sanctuary and three measures of individual-level variation. With the Dirichlet process prior we were able to detect multimodality in each parameter.


Introduction
Heterogeneity is the rule rather than the exception in nature. Variation between individuals in behaviour, size, physiology or almost any other trait is fundamental to the evolution and dynamics of biological systems (Wilson and Nussey, 2010). Despite this truism, ecological models most often deal with the idealised average individual (Bolnick et al., 2003). While this abstraction from real populations has underpinned much of the key advances in understanding population dynamics, there are instances where it becomes necessary to consider the heterogeneity inherent of ecological systems.
A fundamental tool in understanding and assessing the demographics of real populations are mark-recapture studies. These involve following a sample of population members through time to infer abundance and/or survival rates. Individuals are captured, marked and released. At later instances in time, new samples of the population are obtained (e.g. re-captures or re-sightings) (Seber, 1982). Recaptures of previously captured individuals can be used to infer survival, and the ratios of recaptures to captures of new individuals can be used in estimating abundance. In this paper we restrict ourselves to the former within a multi-state recapture framework (Lebreton et al., 2009). Here individuals may transit between various latent or partially observed states. These may include states denoting survival status (i.e. dead, alive but unobserved, recaptured dead etc.) or in demographic states (e.g. life stages or age-classes). Such rates of transitions between stages form the basis of estimates used to populate classical population models (e.g. Leslie matrices and similar).
Typically, most mark-release recapture modeling treats individuals as homogeneous, conditional on state. Analysis approaches for mark-release recapture data which have explicitly attempted to account for individual heterogeneity, have most often involved either the use of a pre-set functional form (e.g. a Gaussian), or assignment of individuals to a prespecified number of groups (Pledger et al., 2003). Having to make assumptions about the number of groups a priori can result in model selection problems in determining the number of groups (Cubaynes et al., 2012). The use of any pre-set form is limited and limiting, as it enforces strict assumptions on the expected distribution of the population.
The subgroups often seen in mark-recapture studies and the complexity of real markrecapture data means that both parametric and discrete style models can be insufficient. This paper tackles this problem by considering a non-parametric approach, the Dirichlet process prior for modelling latent individual heterogeneity. The Dirichlet process prior is a flexible extension to a parametric model as it avoids assumptions about the functional form of the distribution, and it extends discrete style models to the infinite limit by avoiding any prespecifications about the number of groups (Dorazio et al., 2008;Navarro et al., 2006). Despite the appeal of the Dirichlet process prior, it has had little application in markrecapture analysis perhaps because of its complexity and somewhat confusing literature. One exception is Dorazio et al. (2008) who used the Dirichlet process prior to model animal abundance where heterogeneity in abundance between sites was poorly understood and not directly observable.
We present a Markov chain Monte Carlo (MCMC) sampler for a hierarchical Bayesian hidden Markov model, applied to mark-recapture data, which allows for individual heterogeneity in both the observation and process components of the model. In doing so we extend the approach presented by (Ford et al., 2012) into a fully Bayesian and more flexible approach. The methods we present therefore generalize the existing approaches to markrecapture and individual heterogeneity data and provide a new set of tools for understanding both accounting for individual heterogeneity in order to derive more robust inferences about populations and also for quantifying the nature and extent of individual heterogeneity in real populations.

The Dirichlet process prior
The Dirichlet process was first introduced by Ferguson (1973). Several well known methods for the representation of a Dirichlet process include the Polya urn scheme (Blackwell and MacQueen, 1973) or Chinese restaurant process (Pitman, 2006) and the stick-breaking prior (Sethuraman, 1994;Ishwaran and James, 2001). Following is a description of the Chinese restaurant process which is the basis of the algorithm used in this paper.
Consider the analogy of a Chinese restaurant with infinite seating capacity. The first customer enters the restaurant and sits at table one with probability one. Each subsequent customer entering the restaurant chooses a table with probability proportional to the number  of people already seated at the table, or a new table proportional to the precision parameter α (this parameter is described in detail below). Customers at the same table are served the same dish; customers at new tables are served a new dish at random. In this sense, individuals at each table receive the same parameter value (dish), and the table number indicates their cluster membership. In general terms, this means that the probability of seeing an already seen cluster is proportional to the number of individuals in that cluster, and the probability of seeing a new cluster is proportional to the precision parameter α (see Figure 1).
The Dirichlet process is a stochastic process defined as a distribution on distributions and is defined by two quantities: the base distribution G 0 , and the precision parameter α. Although the base distribution may be continuous, individual draws G from the Dirichlet process are discrete with probability one (Blackwell and MacQueen, 1973;Ferguson, 1973;Neal, 2000;Sethuraman, 1994). This means that draws from a Dirichlet process will be clustered on a countably infinite set of discrete values. The result is that values will be repeated, as individuals in the same cluster will have the same value. The lower α is, the more variability will be observed between individual realisations, and for any given realisation a small α will correspond to a smaller number of clusters (see Figure 1). The influence of α on the number of clusters can be seen in Figure 1,with the number of clusters increasing with α, along with the concentration of draws around G 0 for large α. The number of clusters will tend to ∞ with high values of α; conversely the number of clusters will tend to 1 with low values of α. In comparison to this non-parametric approach, finite mixtures must specify the number of clusters a priori. As such, as α tends to infinity, the Dirichlet process is the limit of the discrete groups approach which assumes a fixed number of groups. In this way α corresponds to the strength of prior belief in the base distribution and the number of groups, or clusters, which are likely to be sampled from it. Note that G 0 itself will generally be of specified parametric form, e.g. Normal, and will have unknown parameters which are updated separately to the Dirichlet process.
A generic Dirichlet process takes the form Here, we assume that data y i are independent conditional on θ i , and G is the mixing distribution over θ which has Dirichlet process prior DP (G 0 , α).
MCMC algorithms are the most common approach for inference in Dirichlet processes. Neal (2000) presented several algorithms which use the Chinese restaurant process approach to sample from the posterior distribution of the Dirichlet process. This paper incorporates one of the algorithms developed by Neal (2000). Alternative samplers include: blocked Gibbs sampling using the stick-breaking representation (Ishwaran and James, 2001); updates using a Metropolis-Hastings framework (Jain and Neal, 2004;Liang et al., 2007); and sequential Monte Carlo (Fearnhead, 2004).

Two-state model used for simulations
As usual with hidden Markov or multi-state models, the overall model is split into a process part and an observation part. For the process model, we assume that at time t an animal can be in either of two states S it : Here and Away, or H/A for short. Changes in the state over time are governed by a Markov process with transition matrix γ, so (omitting dependence on i for now) for any two states s and s * we have The four elements of γ can be written in terms of just two parameters γ HH and γ AA (respectively the probabilities of staying Here and staying Away), as follows: For the observation model, there are "capture attempts" (photo-ID expeditions) at each t, in which an animal may be seen if and only if it is Here. Our data for animal i are thus a time series X i,t 1i :T of 0 s (not seen) and 1 s (seen) where t 1i denotes the first observation of the animal (see below) and T the most recent expedition. If X it = 1 then we know S it = H, but if X it = 0 the state cannot be determined for certain. Formally, the probability of observation given state is expressed in terms of a parameter π by We start each animal's series at its first sighting of the given year, and condition on S t 1i i = 1. For synthetic data used in this paper, we assume no recruitment and simulate data with all animals present and seen on the first occasion.

North Atlantic humpback whale data
The methods developed here are applied to a mark-resight data set on a subpopulation of North Atlantic humpback whales sighted in the Stellwagen Bank National Marine Sanctuary (SBNMS), in the Gulf of Maine. Researchers from the Provincetown Centre for Coastal Studies began documenting North Atlantic humpback whales in the Gulf of Maine in 1975 and have to date individually identified over 1200. Humpback whales (Megaptera novaeangliae default) are distributed worldwide, with summer feeding ranges in mid to high-latitudes and winter breeding in low-latitude areas (Clapham and Mead, 1999). They can be uniquely identified by their natural markings: through the shape of their flukes and through patterns from natural pigmentation (Hammond, 1986).
The majority of the North Atlantic humpback whales breed over winter in the West Indies; a small number are thought to use the breeding grounds around the Cape Verde Islands (Stevick et al., 1998).
During summer, the whales disperse to six summer feeding regions. Although historically treated as a single stock, the six summer feeding regions in the North Atlantic hold relatively discrete subpopulations (Clapham and Mayo, 1987), with individuals demonstrating strong site fidelity to a particular feeding region over many years. Feeding sites include the Gulf of Maine, eastern Canada, west Greenland, Iceland and Norway (Katona and Beard, 1990), and patterns of movement suggest perhaps four distinct subpopulations (Stevick et al., 2006).
Individual humpback whales show high maternally directed site fidelity to these summer feeding ranges, as calves follow their mothers from breeding to feeding grounds (Clapham and Mayo, 1987).
The Gulf of Maine is the southern most summer feeding ground for the North Atlantic humpback whales. Individual humpback whales have been intensively studied in this region since the late 1970s. The SBNMS is one of several important feeding sites for North Atlantic humpback whales which summer in the Gulf of Maine. Due to the consistent aggregation of humpback whales and other marine life, the SBNMS was nominated as a national sanctuary in 1992. This area is not only an important feeding ground for the North Atlantic humpback whales, but is also a busy recreation and transportation area for humans with high levels of commercial and recreational vessel traffic. This overlap has resulted in many injuries to the whales from ship collision and entanglement in fishing gear (Robbins and Mattila, 2004).
Although both commercial and recreational fishing are allowed in the sanctuary, regulations have been established which prohibit various other activities such as sand and gravel mining. The sanctuary is a managed resource area equivalent to MPA Category VI (Hoyt, 2011;IUCN, 1994).
The SBNMS encompasses only a small part of the Gulf of Maine sub population's summer range, and although some individuals are seen regularly there during the summer, none are thought to remain permanently within its boundaries.

Three-state model used in application to real data
The two-state model above is extended to a three-state model for application to real data. A three-state hidden Markov model including death, developed in Ford et al. (2012), was applied to data from 237 mature North Atlantic humpback whales. In order to draw out a particular instance of heterogeneity in this population we considered only animals seen more than once after the first 8 seasons of data. This is because our primary aim here is to demonstrate heterogeneity within a set of seemingly alike individuals. Younger animals may well be better modeled by the inclusion of age or sex dependent covariates. The three-state model (including death) uses the same implementation as the two-state model described above.
With three states the nine elements of the transition matrix can be written in terms of just three parameters γ HH , γ AA , and γ D (respectively the probabilities of staying Here, staying Away and Dying in a week) as follows: where it is assumed that the probability of death (which is very low relative to the other transition rates) does not depend on whether the animal is Here or Away.
As sighting effort is focused in the middle of the year, we included all sightings from the 18th week of the year through to the 43rd week. The probability of survival, P surv , over the remaining 26 week period was calculated as An extra parameter q was introduced for the probability of being present in the marine sanctuary at the start of the season. We calculated the probability of each state in the first week of the new year to be: where P(S t−1 ) is the vector of state probabilities in the last week of the previous year.

Estimation
Given a series of observations X 1:T and prior distributions on π and γ, our aim is to estimate the posterior distribution using MCMC. The MCMC routine developed in this paper involves four main steps (five in application to real data).
Individual-level random effects were included on each of π, γ HH and γ AA and are updated using the Dirichlet process prior. We assume individual-level parameters to be consistent over time but have allowed for population-level annual variation (β yr ) in probability of remaining Here using logit-links: logitγ HH i,yr = β yr + γ HH i . Updates to β yr , γ D (death) and q are assumed to be fixed (not individually variable).
One iteration of the MCMC algorithm consists of the following steps: 1. Sampling the hidden state chain for all individuals.
2. Calculating summary statistics per individual conditional on its sampled states.
3. Updating the posteriors for individual-level parameters π i , γ HH i and γ AA i separately using Gibbs sampling from the Dirichlet process prior.

Updating the base distribution and precision parameter:
(a) Updating the base distribution G 0 using an Independent Metropolis-Hastings sampler with three proposal distributions whose parameters vary across iterations (b) Updating the precision parameter 5. Updating population-level fixed effects using an Independent Metropolis-Hastings sampler with a fixed proposal distribution: a multivariate t-distribution whose mean and variance are set using a preliminary fit from ADMB (see (Ford et al., 2012)).

Forward-Backward recursion
In order to update individual-level parameter values (θ) at each iteration, we require counts of successes and trials for each individual. These counts are obtained from the hidden state chains which are sampled using the Forward-Backward recursion scheme defined by Scott (2002) and described by Zucchini and MacDonald (2009). This recursion scheme starts by producing a forward probability vector α 2 , ..., α n , containing the probabilities of the underlying hidden states for each observation given all observed data up to time t. We calculate these forward probabilities, from 1 : T (T being the length of the observation history), for each state, given the observed data (X).
where P(X t |S t ) denotes the probability of the data given the state. Working backwards, we generate a sample path Z (T ) of the Markov chain in the order t = T, T − 1, T − 2, ..., 1, making use of the following proportionality argument: The second factor in equation 1 is simply a one-step transition probability in the Markov chain.

Counts of successes and trials per individual
Observations for an individual are assumed Binomial with probability π i . As the Beta prior for π is conjugate to the Binomial, the posterior is also Beta. For the probability of observation there is a trial whenever an animal is Here; the outcome is whether it was or wasn't seen. There is no trial when then animal is Away, since it is then guaranteed not to be seen. The counts of successes and trials for the transition probabilities (γ HH and γ AA ) are calculated from the sampled state chains. For γ HH , there is a trial whenever the animal was Here (excluding the final period); the outcome is whether it stayed Here or not. A similar scheme applies to γ AA .

Gibbs sampling via the Dirichlet process prior.
The individual-level random effects (π i , γ HH i and γ AA i ) are updated separately using a Dirichlet process prior which follows algorithm 8 by Neal (2000) (see algorithm 1).
The algorithms in Neal's paper (2000) work by assigning individuals to clusters. Due to the clustering property of the Dirichlet process, some of the individual parameter values θ i will be identical, and each θ i is associated with a cluster. Indicator variables c i are used to indicate the current cluster membership for each individual (which may change over the Algorithm 1 Algorithm 8 by Neal (2000) Let the state of the Markov chain consist of c = (c 1 , ..., c n ) and φ = (φ c : c ∈ c 1 , ..., c n ). Repeatedly sample as follows: • For i = 1, ..., n [where n indicates the number of individuals]: Let k − be the number of distinct c j for j = i, and let h = k − + m. Label these c j with values in 1, ..., k − . If c i = c j for some j = i, draw values independently from G 0 for those φ c for which k − < c ≤ h. If c i = c j for all j = i, let c i have the label k − + 1, and draw values independently from G 0 for those φ c for which k − + 1 < c ≤ h. Draw a new value for c i from 1, ..., h using the following probabilities: where n −i,c is the number of c j for j = i that are equal to c, and b is the appropriate normalizing constant. Change the state to contain only those φ c that are now associated with one or more observations.
• For all c ∈ c 1 , ..., c n : Draw a new value from φ c | y i such that c i = c, or perform some other update to φ c that leaves this distribution invariant.
course of the MCMC) and the clustering of individuals means that the number of active clusters will typically be much smaller than N ; K is used to refer to the number of active clusters. For k = 1, ..., K, each cluster c k will have associated parameter value φ k .
Algorithm 8 in Neal's paper (2000), the one implemented here, allows for efficient Gibbs sampling with a non-conjugate distribution. At each iteration, the algorithm temporarily includes m auxiliary components; these are new potential values for clusters, which may or may not actually get individuals assigned to them. For each individual, when updating c i , either an existing cluster is chosen or one of these m new components. The probability of joining an existing cluster will be proportional to the number of individuals in that cluster, and the probability of joining a new cluster will be proportional to α/m, the prior precision split equally among the auxiliary components. These auxiliary components are generated i.i.d from the base distribution and are discarded at each iteration if not used by the Gibbs sampler (i.e not chosen as a new cluster). The use of auxiliary components avoids the need to integrate with respect to the distribution G 0 as these auxiliary components represent the new possible components. This approach is similar to methods developed by MacEachern and Muller (1998) in that auxiliary components are used to update the model, with the difference that the auxiliary components exist only temporarily in Neal's algorithm.
Following Algorithm 8 in Neal's paper (2000) (see algorithm 1), individual parameter values for π i , γ HH i or γ AA i are updated by generating and assigning new clusters. For each class, c, the parameter φ c determines the associated probability for that class; the collection of all φ c is denoted by φ. In algorithm 1, F (y i , φ c ) is calculated as the density under a Binomial and c i indicates which latent class is associated with observation y i , where the numbering of c i is of no significance.

Updating hyper-parameters (a, b) for the base distribution G 0 and fixed effects
In order to update the hyper-parameters governing the base distribution G 0 and any populationlevel fixed effects, it suffices to use the machinery for the Independent Metropolis-Hastings sampler developed in Ford et al (in submission), which uses a proposal distribution derived from a logit-Normal approximation to the conditional posterior of (a,b). For reference we have included an appendix describing this method.

Updating the precision parameter α
Despite its importance, there is a lack of agreement in the literature outlining efficient methods to update the precision parameter (Dorazio, 2009;Kyung et al., 2010;Navarro et al., 2006;Escobar and West, 1995). A Gamma(a, b) prior is commonly used due to its conditional conjugacy property. However, the problem is knowing how to efficiently update the Gamma hyper-parameters (a, b). The most recent and concise work in this field is by Murugiah and Sweeting (2012) who propose values for the hyper-parameters which can be used in the presence or absence of information. They suggest that standard use of small a and b can result in high posterior weights for k = 1 and k = n,where k is the number of clusters. Instead they propose an alternative method which results in a = b = exp(−0.033n), giving a prior mean of unity with increasing standard deviation with larger n. The appeal of the method by Murugiah and Sweeting (2012) is that the prior gives less rigid adherence to G 0 with more data. In cases with small n it will generally be futile to search for, e.g. multimodality, so there is no gain in allowing overly flexible realisations of G 0 .
We combine work by Escobar and West (1995) and Murugiah and Sweeting (2012) to update the precision parameter: methods developed by Murugiah and Sweeting (2012) to update the hyper-parameters are incorporated into the sampling framework developed by Escobar and West (1995). Escobar and West (1995) describe how α can be updated by incorporating an auxiliary variable η into the Gamma prior. The formula for updating α is expressed as a mixture of two gamma posteriors, with the conditional mixing parameter for α and k, a simple Beta distribution. They found that p(α|k) is the marginal distribution from a joint distribution for α and continuous quantity, η, such that where η is sampled from a Beta distribution: (η|α, k) ∼ B(α + 1, n). Taking the conditional posteriors p(α|η, k) ∝ α a+k−2 (α + n)e −α(b−log(η)) ∝ α a+k−1 e −α(b−log(η)) + nα a+k−1 e −α(b−log(η)) when α > 0 this reduces to a mixture of two Gamma densities where π = (aα+k−1) (n(bα−log(η))+aα+k−1) .

Simulation testing
The two-state model was used to test the Dirichlet process using a synthetic data set with 30 individuals, each with 1000 length capture history. We assumed individuals came from (randomly) one of two discrete groups: π= 0.82 or 0.96; γ HH = 0.88 or 0.98; γ AA = 0.8 or 0.95. Three separate chains were run for 15000 iterations. The chains were arbitrarily thinned to every 2nd update and combined to form one chain of 22500 posterior samples. The chains were thinned to reduce any auto correlation between successive samples (Gilks et al., 1996).
The posterior density of log(α) (Figure 2) displays standard unimodal form. The posterior distribution of k (the number of clusters) indicated two clusters for each parameter. Figure  3 displays the posterior density for each of π i , γ HH i and γ AA i , with posterior estimates for each parameter clustered around the two true values used for simulation.

Limits of Dirichlet process prior
The following example is intended to highlight the potential limits of the Dirichlet process in identifying clusters. Data was simulated for 30 animals with 1000 length capture history and run for 15000 iterations, with the first 5000 discarded due to burn-in. Three groups were assumed for both π= 0.6, 0.85, 0.96 and γ HH = 0.5, 0.8, 0.95, and two groups for γ AA = 0.89, 0.97. Individuals were randomly assigned to a group for each parameter. Figure 4 indicates the inability of the Dirichlet process prior to distinguish between low π and low γ HH . The results indicate that the lowest true group in π (p = 0.6) could not be distinguished, and that the lowest estimated group in γ HH was lower than the actual true values used in data simulation (γ HH = 0.5). At higher probabilities the posterior density of parameters appeared to cluster around the true discrete values used in data simulation. This result is unsurprising due to uninformative data and the resulting inability to distinguish between not being present and not being seen.

Unimodal distributions
One concern with the use of Dirichlet process prior is the potential for spurious multimodality when in fact none is present. To investigate whether this is likely to be a problem we generated 10 synthetic data sets of 30 individuals each with 1000 length capture history. For each parameter (π, γ HH and γ AA ), synthetic data was simulated from a Normal distribution with low variance, N (2, 0.1). The MCMC algorithm was run for 10000 iterations. There was no evidence of bi-modality in the results ( Figure 5). The results of this simulation experiment therefore suggest that spurious multimodality given a truly unimodal distribution is unlikely.

North Atlantic humpback whale data analysis
One chain was run for 25000 iterations with the first 5000 discarded to burn-in. Density plots for the log of α the precision parameter, indicate expected unimodal density (Figure6). The posterior distribution of the number of clusters indicated more variation for π compared to both γ HH and γ AA . Figure 7 indicates some multimodality for each of π, γ HH and γ AA . In each case, low probabilities corresponded to whales seen only a few times. With such an uninformative data history it is difficult to distinguish between not being seen and not being present. As expected, more observations corresponded to higher probability of observation and presence in the marine sanctuary.

Discussion
In some studies, covariates may adequately explain the majority of individual heterogeneity present in the data. However, in some cases and for certain species (for example cetaceans as considered here), it is unrealistic to expect to be able to collect all necessary covariates, or even to know which covariates would likely explain the heterogeneity. Nonetheless,  Figure 7: Density for π, γ HH and γ AA for 25000 iterations for 237 whales.
adequately capturing latent heterogeneity is important to ensure accurate analysis of markrecapture data. The use of the hidden Markov model combined with the Dirichlet process prior, provides a powerful tool for capturing latent individual heterogeneity.
Using simulations, our results show the Dirichlet process prior was able to accurately capture multimodality in three measures of individual heterogeneity: probability of observation, probability of remaining in the marine sanctuary and probability of remaining away. Through simulation studies we were able to explore the accuracy, and limits, of the Dirichlet process prior to distinguish multiple groups using this framework. We found in certain areas of parameter space, the Dirichlet process prior was capable of capturing up to three distinct groups. However, the model was not capable of distinguishing between low probability of observation and low probability of remaining in the marine sanctuary. This aliasing is unsurprising and is due to the lack of information contained in the capture histories (Ford et al., 2012), rather than due to the Dirichlet process prior.
In application to North Atlantic humpback whales, we found evidence of multimodality apparent in each parameter. As expected, we found low posterior probabilities corresponded to whales seen only a few times. However, with uninformative data histories identifiability issues are expected as the model cannot discern between individuals not being seen or simply not being present. The variation in both the state transition probabilities implies substantial differences in proportion of time spent in the marine sanctuary. This estimate has implications in the ability to predict the long term usage of the marine sanctuary and for population survival and growth. Whilst this extra uncertainty may have implications for the understanding the population's usage of the marine sanctuary, it is worth considering what would be inferred from a fixed effect model under similar circumstances. In this case it would be likely to be overestimated (Ford et al., 2012) compared to the results from the model here with the Dirichlet process prior.
There are several extensions and applications of the Dirichlet process which were not explored here but are important considerations and interesting areas for future exploration. Additional interesting applications could involve further exploration of correlations between individual random effects: for example, the multiple behavioural modes indicated that individuals who were often away were more likely to be infrequently observed. With the addition of random effects onto arrival time each year it would be interesting to see the correlation between arrival and departure, and arrival and length of stay in the marine sanctuary. In future research, it would be worthwhile investigating the ability of the Dirichlet process to model this, or other, correlations in behaviour.
In comparison to parametric distributions, the Dirichlet process allows for multiple modes in both the observation and state process. Heterogeneity in detection in mark-recapture data has been a hurdle in the accurate estimation of abundance. With the potential to identify multiple modes in the probability of observation, the Dirichlet process has the potential to give more accurate estimates of abundance. The Dirichlet process also has important application to more effective marine spatial planning as it provides a method to more accurately capture the individual behaviour, which translates into more accurate estimations of proportion of time spent in the marine sanctuary.
The development of Bayesian hierarchical models has been the focus of much effort in mark-recapture research (King, 2012). Despite this, non-parametric approaches have received little attention. We have presented a hierarchical hidden Markov model which allows for both process and observation error and have incorporated the Dirichlet process prior to account for individual heterogeneity on both the observation and process components. We anticipate that this powerful addition to mark-recapture analysis will be useful in application to other problems by allowing for accurate estimation of multiple behavioural modes.