1 2 3 4 5 6 7 8 9
20210307_Fencing_RGB_40m 0.350 0.326 0.293 0.301 0.303 0.172 0.013 0.302 0.346
20210307_Fencing_RGB_60m 0.246 0.195 0.193 0.203 0.002 0.009 0.130 0.141 0.185
10 11 12 13 14 15 16 17 18
20210307_Fencing_RGB_40m 0.340 0.345 0.352 0.371 0.31 0.188 0.101 0.024 0.021
20210307_Fencing_RGB_60m 0.265 NA NA NA NA NA NA NA NA
19 20
20210307_Fencing_RGB_40m 0.014 0.005
20210307_Fencing_RGB_60m NA NADetailed model formulation, analysis, and results
Supporting Information for: Estimating red deer (Cervus elaphus) population density using drones in a steep and rugged terrain (DOI: 10.1002/wlb3.01533)
This document presents a detailed model description, details about the model results, including convergence diagnostics, posterior predictive checks, tables and figures. The document is written in Quarto, and the source file containing all executable R code is available in the GitHub-repository accompanying the paper, https://github.com/torbjore/dronedeer.
1 Detailed model formulation
Deer densities were estimated using hierarchical state space models fitted by the use of Bayesian Markov Chain Monte Carlo (MCMC) sampling in NIMBLE ver. 1.3.0 (Perry de Valpine et al. 2017, 2024). We combined a double-observer observation model with an overdispersed Poisson process model describing the variation in deer densities among sampling sites (focal images). We will here first derive the observation model and the process model, using a hierarchical formulation, and then provide details about the analysis. The NIMBLE model code, together with the data, are provided on GitHub (https://github.com/torbjore/dronedeer).
To estimate abundance and population density, one must account for the fact that not all individuals that are present will be observed (i.e., imperfect detection). Detection failure can be due to either individuals being hidden (e.g, under trees) and hence not available for registration in the images, or because they were missed by the observers while inspecting the images even though they were present (available for detection) in the images (Brack, Kindel, and Oliveira 2018). By comparing registrations from multiple overlapping images, we found that it was very unlikely that an individual would be completely hidden in all overlapping images (see main text), and hence we did not account for imperfect availability. However, it was not uncommon that an individual registered by one observer was missed by the other observer. Hence, we used a double observer protocol to estimate detection probabilities (Nichols et al. 2000; Williams, Nichols, and Conroy 2002):
At each site (focal image) \(i\) we have records of how many individuals that were seen by only observer 1 (\(y_{10,i}\)), how many individuals were seen by only observer 2 (\(y_{01,i}\)), and how many individuals were seen by both observers (\(y_{11,i}\)). We denote the observed data from each site by a vector \(\mathbf{y}_i = \left[y_{01,i},y_{10,i},y_{11,i} \right]\), and the number of individuals seen by at least one observer is the sum of the elements in this vector, \(\mathbf{Y}_i = y_{01,i} + y_{10,i} + y_{11,i}\). The observation model defined in the first section below describes the probability (likelihood) of getting these observed data, conditional on the number of individuals (\(N_i\)) present at each site; i.e. \(L(\mathbf{y}_i |N_i)\). The number of individuals at each site, \(N_i\), is unknown, and the stochastic model for this is described in the next section.
1.1 Observation model
Given that the two observers detect individuals independently with probabilities \(p_1\) and \(p_2\) respectively, the probability that an individual is detected by only observer 1 is \(P_{10} = p_1 (1-p_2)\), the probability of only being detected by observer 2 is \(P_{01} = (1-p_1)p_2\), and the probability of being detected by both observers is \(P_{11} = p_1 p_2\). The probability of being detected by at least one observer is one minus the probability of not being detected by any of the observers, \(p^* = 1-(1-p_1)(1-p_2)\). Note that this is equivalent to \(p^* = P_{10} + P_{01} + P_{11}\).
The key to derive a model for the observed data (\(L(\mathbf{y}_i│N_i)\)), is to first write the probability of the observed data conditional on the number of individuals seen by at least one observer (\(Y_i\)), and then multiply this with the probability of observing \(Y_i\). I.e., \(L(\mathbf{y}_i│N_i) = L(\mathbf{y}_i|Y_i)L(Y_i|N_i)\). The first of these components is a multinomial likelihood, and the second component is a binomial likelihood. Hierarchically, we write this as
\[ \begin{align} \mathbf{y}_i &\sim \text{Multinomial} (\pi,Y_i) \\ Y_i &\sim \text{Binomial}(p^*,N_i) \end{align} \]
where \(\pi=[P_{10},P_{01},P_{11} ]/p^*\) (i.e., the multinomial probability vector \(\pi\) contains the probabilities of the three detection categories conditional on detection). Note that all observed individuals must belong to one, and only one, of the three detection categories, and hence the elements of \(\pi\) sum to one. Note also that when \(Y_i=0\), the multinomial probability of observing \(\mathbf{y}_i=[0,0,0]\) becomes 1 regardless of the parameter values.
This initial observation model has only two top level parameters: the detection probabilities of the two observers, \(p_1\) and \(p_2\) (these could have been made site specific, and then been modelled as a function of covariates on e.g. a logit scale). Posterior predictive checks (see below) indicated, however that the two independent observers had sufficiently variable detection probabilities among sites to call for a model with random detection probabilities – see details in Section 2.2.
1.2 Process model
The number of individuals present at a site was modelled as an overdispersed Poisson process. Hence, we assume that number of individuals present at site \(i\), \(N_i\), is drawn from a Poisson distribution with expectation \(\lambda_i a_i\), where \(\lambda_i\) is the expectation per unit area (we used hectares) and \(a_i\) is the area of the site,
\[ N_i \sim \text{Poisson}(\lambda_i a_i) \]
Overdispersion was accommodated by modelling \(\lambda_i\) as a random variable. As we anticipated that density estimates may be somewhat sensitive to the choice of the distribution of \(\lambda_i\), we examined three alternative random distributions; log-normal, exponential, and gamma.
1.2.1 Log-normal expectation
For the log-normal model we used a standard log-linear log-normal model,
\[ \begin{align} \log(\lambda_i) &\sim N(\mu_i, \sigma^2) \\ \mu_i &= x_i \beta \end{align} \]
where \(x_i\) is a vector of site-covariates and \(\beta\) is a vector of the linear coefficients (including an intercept). Mean deer density in the log-normal model is then \(e^{\mu_i+\sigma^2/2}\) and the variance is \((e^{\sigma^2}-1) e^{2\mu_i+\sigma^2}\).
1.2.2 Exponential expectation
The exponential distribution is a one-parameter distribution with probability density function \(f(\lambda) = \Lambda e^{-\Lambda\lambda}\), expectation \(1/\Lambda\), and variance \(1/\Lambda^2\). We modeled the expectation as a log-linear function of covariates (i.e., assuming multiplicative effects), and hence formulated the model hierarchically as
\[ \begin{align} \lambda_i &\sim \text{Exponential}(\Lambda_i) \\ \Lambda_i &= e^{-\mu_i} \\ \mu_i &= \mathbf{x}_i \boldsymbol{\beta} \end{align} \]
Mean deer density in the exponential model thus becomes \(e^{\mu_i}\) and the variance is \(e^{2\mu_i}\).
1.2.3 Gamma expectation
The gamma distribution parameterized with a shape parameter \(a\) and a rate parameter \(b\) has probability density function \(f(\lambda) = \frac{b^a}{\Gamma(a)} \lambda^{a-1} e^{-\beta\lambda}\), expectation \(\frac{a}{b}\), and variance \(\frac{a}{b^2}\) (the variance is proportional to the mean). We initially tried keeping the rate parameter constant and using a log-linear covariate model for the expectation (i.e., \(a_i = be^{\mathbf{x}_i\boldsymbol{\beta}}\)), but this model had very poor mixing, even when using blocked mcmc-sampling (perhaps because \(b\) is a parent of \(a\)). Hence, we instead kept the shape parameter constant and let the rate parameter depend on \(\mu_i\) and \(\sigma\) in such a way that the variance-mean relationship is the same as in the log-normal model (the variance is proportional to the square of the mean),
\[ \begin{align} \lambda_i &\sim \text{Gamma}(a,b_i) \\ b_i &= ae^{-(\mu_i+\sigma^2/2)} \\ a &= \frac{1}{e^{\sigma^2} - 1} \\ \mu_i &= \mathbf{x}_i \boldsymbol{\beta} \end{align} \]
Hence, mean and variance of deer density in the gamma model has the same expressions as in the log-normal model, respectively \(\frac{a}{b} = e^{\mu_i+\sigma^2/2}\) and \(a/b^2 = \left( e^{\sigma^2}-1 \right) e^{2\mu_i+\sigma^2 }\). Compared to the log-normal distribution with the same mean and variance, the gamma distribution has a higher density close to zero and a lower density far out in the tail of the distribution.
1.3 Modelling of covariate effects
In our application, we included the distance from the agricultural field as the only site covariate in \(\mathbf{x}_i\). The covariate effect was the same in all study area-month combinations, but the intercept was allowed to vary across area-months. Denoting the index for the area-month that site \(i\) belongs to as \(j(i)\), we can formulate the covariate model as,
\[ \mu_i = \mu_{0,j(i)} + x_i \beta \]
In principle, \(\mathbf{x}_i\) can include both covariates that vary among sites within study areas and covariates that vary only among, and not within, study areas. However, it is important to realize that if \(\mathbf{x}_i\) includes study area covariates, random effects of study areas would need to be included in order to avoid pseudoreplication (sites within study areas is not a random sample of sites with the same study area covariate value).
2 Modeling focus and strategy
2.1 Specification of priors and starting values
For the observation model parameters \(p_1\) and \(p_2\) (detection probabilities) we initially used \(\text{Uniform}(0,1)\) priors. However, as posterior predictive checks indicated that random effects on detection probabilities were required, we decided to use logit-normal priors for detection probabilities (see Section 2.2). We also used data from trials at a deer enclosure in a natural habitat to inform the priors for these parameters (see Section 2.3).
The covariate (distance from field) was normalized by subtracting the mean and dividing by the standard deviation prior to model fitting, and we used a \(\text{N}(0,\sigma=2)\) prior for the covariate effect, \(\beta\).
One approach to the modelling could have been to specify somewhat informative priors for the intercepts \(\mu_{0,j(i)}\) based on prior knowledge of reasonable deer densities in the area. However, due to the methodological nature of our study, we decided to use essentially uninformative priors for \(\mu_{0,j(i)}\), \(\text{N}(0,\sigma=10)\).
To find reasonable starting values for the MCMC-chains, we first computed estimates of \(p_1\) and \(p_2\) as the proportion of those individuals that were seen by observer 2 that were also seen by observer 1 and vice versa. Simplified estimates of population densities were then computed from the deer counts, the estimated areas of the sites, and the initial estimates of \(p_1\) and \(p_2\) (see code for details <link to GitHub repository will be inserted after peer review>).
For the log-normal and the gamma models, the prior distribution for \(\sigma\) was specified as
\[ \sigma \sim \text{Gamma}(\text{shape} = 0.01, \text{ rate} = 0.01) \]
In Section 4.3.8, we present computed overlap between prior and posterior distributions for all parameters.
2.2 Posterior predictive checks and model modifications
We performed posterior predictive checks (Gelman et al. 2004; Conn et al. 2018) with respect to both \(Y_i\)’s and \(\mathbf{y}_i\)’s by comparing the sum of the squared Pearson residuals (discrepancy measure) for observed and simulated data in every step of the MCMC-chains. The discrepancy measures for \(\mathbf{y}_i\) were conditional on the observed values of \(Y\), and hence a lack of fit with respect to \(Y_i\)’s will not lead to a lack of fit with respect to \(\mathbf{y}_i\). Bayesian p-values (proportion of the discrepancy measure that are larger for the simulated data than for the observed data) did not indicate any lack of fit with respect to \(Y_i\)’s (\(p > 0.4\)), but there appeared to be a lack of fit with respect to \(\mathbf{y}_i\)’s (\(p < 0.03\)), indicating that there were more variation in the counts from the two independent observers than expected if detection probability was constant across sites for each observer. To accommodate this variation in the model, we introduced random detection probabilities, assuming the logit-transformed detection probabilities to be normal. Hence, in our final observation model, we used
\[ \begin{align} p_{1,i} &= \frac{1}{1+e^{-(\eta_1 + \epsilon_{1,i})}} \\ p_{2,i} &= \frac{1}{1+e^{-(\eta_2 + \epsilon_{2,i})}} \\ \epsilon_{1,i} &\sim \text{N}(0,\sigma_p)\\ \epsilon_{2,i} &\sim \text{N}(0,\sigma_p)\\ \end{align} \]
A model with a wide prior for \(\sigma_p\) resulted in unreasonable high upper values; the high prevalence of \(Y_i = 0\) in the data could be explained by either a high number of sites having very low expected densities (which is reasonable) or a high number of sites having very low detection probabilities (which is not reasonable). Hence, to facilitate model fitting, we used a somewhat informative prior for \(\sigma_p\); a Gamma-prior truncated at an upper bound of 0.59,
\[ \sigma_p \sim\text{T} \left[ \text{Gamma}(1,0.05), 0, 0.59 \right] \]
The upper bound of \(0.59\) was set based on a prior assessment that detection at the 97.5 percentile site should have a maximum of 10 times higher odds than in the 2.5 percentile site (\(\log(10)/(1.96 \times 2) \approx 0.59\)). In Section 4.3.9.4 we investigate the sensitivity of our results with respect to the choice pf prior for \(\sigma_p\).
Bayesian p-values did not indicate any lack of fit of this final model with respect to \(Y_i\)’s (\(p = 0.50\)) or \(\mathbf{y}_i\)’s (\(p = 0.19\)) (see Section 4.3.6).
2.3 Using data from the deer enclosure to inform priors on detectability
The enclosure surrounded a natural habitat with a similar forest type as at the study areas. It is therefore reasonable to assume that detection probabilities inside and outside the enclosure are similar. It is, however, not reasonable to assume that the variation in deer densities among sites (focal images) are similar because the enclosure had much higher density, much more females than males, many young individuals, and there were feeding stations where the deer tended to aggregate.
Initial analysis showed that there was no significant difference in detection probability in the enclosures for the two observers, and we hence simplified the observation model to have the same mean logit-scale detection probability for the two observers, \(\eta_1 = \eta_2 = \eta\). The model was then fitted with the following uninformative prior:
\[ \eta \sim \text{N}(0,2) \\ \]
When fitting the full model to the data from the study areas, we used as priors
\[ \begin{align} \eta_1 &\sim \text{N}\left(\text{mean}(\boldsymbol{\eta}), 1.5\text{sd}(\boldsymbol{\eta})\right) \\ \eta_2 &\sim \text{N}\left(\text{mean}(\boldsymbol{\eta}), 1.5\text{sd}(\boldsymbol{\eta})\right) \\ \end{align} \] where \(\text{mean}(\boldsymbol{\eta})\) and \(\text{sd}(\boldsymbol{\eta})\) refer to the mean and standard deviation of the posterior samples for \(\eta\) in the enclosure model. The posterior standard deviations were conservatively increased by 50% to account for the fact that detection probabilities inside and outside of the enclosure could be different. See Section 4.2.3 for details.
2.4 Implementation and model selection
The models were fitted by the use of Bayesian MCMC sampling in NIMBLE (v1.3.0; (P. de Valpine et al. 2024; Perry de Valpine et al. 2017)). Convergence was assessed on three independent MCMC chains (with different start values) with the Gelman and Rubin (1992) diagnostics using the ‘gelman.diag-function’ in the ‘coda’ package ver. 0.19-4.1 (Plummer et al. 2006) in R ver. 4.3.3 (R Core Team 2024).
For each of the three models for the random Poisson expectation (log-normal, exponential and gamma), we fitted both linear and quadratic models for \(\mu_i\) as a function of distance from the field, as well as a model with no covariate effects; hence obtaining 9 models. We selected the best model according to the Watanabe-Akaike (Widely Applicable) Information Criterion, WAIC (Watanabe 2010) for presentation of estimates and prediction (the criterion seeks to minimize the prediction error variance as an optimal trade-off between accuracy and precision of the fitted model predictions).
We ran MCMC-sampling on each of the nine models with 3 chains with random and independent starting values for 180,000 iterations after a burn-in of 20,000 iterations and stored every 6th posterior sample (obtaining a total of 90,000 samples for inference). Satisfactory convergence was confirmed by the (Gelman and Rubin 1992) ‘potential scale reduction factor’ (psrf) being <1.05 for all parameters (see Section 4.3.5).
2.5 Estimating densities at the specific sites included in each survey
The estimators for mean density presented in Section 1.2 assume that the sites are random samples of a large number of sites of the same characteristics, and the estimates apply to this statistical population of sites. We also computed estimates of density (at the time of flight) at the specific sites included in each survey as the estimated number of individuals present divided by the total area covered, \(\bar{d}_s =\sum_{i=1}^n \hat{N}_i/\sum_{i=1}^n a_i\) (the sums are over all \(n\) sites in the survey and the subscript \(s\) in \(\bar{d}_s\) indicates that the mean only applies to the specific sites sampled).
2.6 Estimating number of deer in a specific area
When estimating number of deer in a specific area, such as in the deer enclosures, simply computing the posterior of \(N_{\text{Tot}} = \bar{d}_s/q\) were \(\bar{d}_s\) is the mean density at the surveyed sites (Section 2.5) and \(q\) is the proportion of the area included in survey, will underestimate the sampling variance (uncertainty) of the total number of deer in the area since density in the sampled and unsampled part of the area would not be exactly the same.
To compute the sampling variance of mean number of deer at sample units, we must account for the fact that we sample from a finite area without replacement. It is not straightforward to do this in a Bayesian setting (Banerjee 2024), and our case is further complicated by the fact that sample sites vary in size (particularly in the deer enclosures, where parts of the focal images were outside the enclosures). Hence, since we only have a finite study area in the deer enclosure, and we only estimate number of deer in the enclosure as a proof of concept, we took a more frequentist approach to the problem, ignoring the fact that sites are of different sizes:
If \(n\) sites are sampled from an area covering \(K\) sites, the standard expression for the sampling variance of the mean number of deer in the \(K\) sites is \[ \text{Var}(\bar{N}) = \frac{K-n}{K} \frac{s^2}{n} \] where \(s^2\) is is the empirical variance in number of individuals per site. The first factor is the proportion of the \(K\) sites that are not included in the sample, and the variance becomes zero when we have counts of the entire area (when \(n = K\)).
The total number of deer in the area is \(N_{\text{Tot}} = K\bar{N}\) with sampling variance \(\text{Var} \left( N_{\text{Tot}} \right) = K^2 \text{Var} \left( \bar{N} \right)\). Replacing with the proportion of the area covered in the sample, \(q = n/K\), and inserting the expression for \(\text{Var}(\bar{N})\) above, we get
\[ \text{Var} \left( N_{\text{Tot}} \right) = \left( \frac{n}{q} - n \right) \frac{s^2}{q} \] (see also Williams, Nichols, and Conroy (2002) chap. 12.4.1).
In our situation, we do not have exact counts of number of deer per site (due to imperfect detection), but instead we have posterior samples of \(N_i\) at each site \(i\). To assess the uncertainty in \(N_{\text{Tot}}\), we computed posterior samples for lower and upper 95% confidence limits assuming a normal distribution (by the central limit theorem) at the log-scale, using the delta-method approximation,
\[ \text{CI} = \exp(\log(N_{\text{Tot}}) \pm 2 N_{\text{Tot}}^{-1}s) \]
In the analysis of the deer enclosure data, we were not interested in describing the variation among sites, and therefor instead specified uninformative priors for the true number of deer at each site hierarchically as
\[ \begin{align} N_i &\sim \text{Poisson}(\lambda_i) \\ \lambda_i &\sim \text{Uniform}(0, 50) \\ \end{align} \]
3 Settings of UVA and apps
Before take-off, the UAV was calibrated and set to positioning flight mode (“P-mode”). The DJI 4.0-app was used. Settings of the DJI 4.0-app were kept on default mode except for some options that were changed to the following:
Focus: infinity
Shutter priority: S
ISO: auto ON
Shutter: 1/1000
Image size: 3:2
Image format: JPEG
White balance: auto ON
Apart from selection of UAV model, the settings of the Map Pilot-app were identical for the ‘DJI Mavic 2 Pro’ drone and the ‘DJI Mavic 2 Enterprise Dual’ drone. The changes made to the default settings were:
Model: ‘DJI Mavic 2 Pro’ or ‘DJI Mavic 2 Zoom/Enterprise’
Image format: RAW image OFF
Speed class rating: Class 3
Enable radius guide: OFF
4 Detailed analysis and results
4.1 Organization of the data
All data are in the data folder of the repository. The data-files included here are:
SiteData.rda: The full data from the study sites (excluding enclosure). SeeSiteDataInfo.htmlfor description of the variables.CountData.rda: Based onSiteData.rda, but includes only counts in a long format.CountData_Fence.rda: Counts in a long format from the enclosure.nimbleData.rda: Data from the study sites organized for NIMBLE models. File is produced by the script indata\Make_nimbleData.R.nimbleData_Fence.rda: Data from the enclosure organized for NIMBLE models. File is produced by the script indata\Make_nimbleData_Fence.R.
4.2 Enclosures
The enclosure contained a total of 117 red deer in a 5 hectare area of natural forest. Two survey flights were conducted over the enclosures; one at 40 m above ground, and one at 60 meter above ground. We used these data to first test the method, and then to obtain information about detection probabilities, which we used to inform priors for detection probabilities in the models fitted to the wild deer data.
4.2.1 Overview of data
Area of the sampling sites in hectare (ha); rows represent survey at 40 m and 60 m above ground respectively:
The reason why some sites cover a much smaller area than other sites is that large parts of these sites (focal images) were outside the enclosure. The survey from 40 m above ground contained data from 20 non-overlapping sites, covering a total of 4.48 ha (90% of the 5 ha enclosure). The survey from 60 m above ground contained data from 10 non-overlapping sites, covering a total of 1.57 ha (31% of the enclosure). The reason for the lower overall coverage from the 60 m survey is that larger proportions of the focal images from this survey covered ground outside of the enclosure.
Counts at 40 m, \(\mathbf{y}_i\) for each site (column) \(i\) (rows: first number of deer detected by only observer 1, then number of deer detected by only observer 2, and finally number of deer detected by both observers):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
only_ab 0 0 1 5 4 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
only_jb 0 0 2 4 6 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0
both 0 0 2 5 25 12 0 0 0 0 0 0 0 2 0 0 0 0 0 0Counts at 60 m:
1 2 3 4 5 6 7 8 9 10
only_ab 1 1 3 0 0 0 0 0 5 0
only_jb 1 1 0 0 0 0 0 0 1 1
both 0 2 5 0 0 0 1 3 12 44.2.2 Deer numbers
The NIMBLE-model to estimate total number of deer in the enclosure (see Section 2.6) is defined in R/nimble_models/enclosure_model_N.q and the model is fitted with the script in R/run_enclosure_model_N.R (see the script for mcmc trace plots, convergence diagnostics and posterior predictive checks).
Summary of posterior distributions:
Iterations = 1:20000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 20000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
N_sum[1] 74.1866 2.2919 0.0093568 0.031415
N_sum[2] 42.7034 1.6195 0.0066117 0.017952
N_tot[1] 82.8033 2.5581 0.0104436 0.035064
N_tot[2] 136.1121 5.1621 0.0210741 0.057220
N_lower[1] 59.2949 2.2073 0.0090111 0.030687
N_lower[2] 68.1236 3.7530 0.0153214 0.040290
N_upper[1] 115.6390 2.9567 0.0120706 0.038901
N_upper[2] 272.0958 8.4520 0.0345050 0.073830
eta[1] 1.2741 0.3082 0.0012581 0.007189
eta[2] 1.3994 0.3568 0.0014565 0.006676
sigma_p 0.4091 0.1348 0.0005504 0.004798
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
N_sum[1] 71.0000 73.0000 74.0000 75.0000 80.0000
N_sum[2] 41.0000 41.0000 42.0000 43.0000 47.0000
N_tot[1] 79.2466 81.4789 82.5950 83.7111 89.2919
N_tot[2] 130.6828 130.6828 133.8702 137.0576 149.8071
N_lower[1] 56.2978 57.5210 58.9315 60.4565 64.6580
N_lower[2] 64.4909 64.4909 67.2939 69.9641 77.8174
N_upper[1] 111.5499 113.4792 115.1536 117.2773 122.7774
N_upper[2] 264.8124 264.8124 269.6509 276.2757 293.5380
eta[1] 0.6781 1.0686 1.2709 1.4771 1.8914
eta[2] 0.7069 1.1539 1.3904 1.6445 2.1011
sigma_p 0.1037 0.3198 0.4377 0.5211 0.5837Here, the indices 1 or 2 in brackets represent survey at 40 m and 60 m respectively. N_sum is the sum of individuals in the surveyed sites (see section Section 2.5). N_tot is the estimated total number of deer in the enclosure, \(N_\text{Tot}\), as defined in section Section 2.6. The posterior distribution for this random variable is too narrow as it assumes that density in the parts of the enclosure that has not been surveyed is exactly the same as in the surveyed part of the enclosure. N_lower and N_upper are the confidence limits derived in Section 2.6. Finally, the three last parameters are the mean and standard deviation of the logit scale detection probabilities.
Judging by the posterior median, the approximate frequentist confidence intervals derived in Section 2.6 marginally fail to cover the true known number of deer in the enclosure (117 deer) in the 40 m above ground survey (confidence interval range from about 59 to 115 individuals). The true value is however included in the much wider confidence interval for the 60 m above ground survey (confidence interval range from about 67 to 270 individuals). Even though the credible intervals for the number of deer at the surveyed sites (posterior distributions for N_sum) were not discouragingly wide, the confidence intervals for the total number of deer in the enclosure (based on N_lower and N_upper) were rather wide, particularly for the 60 m above ground survey. This is due to the clustered distribution of deer within the enclosure, combined with somewhat limited sampling effort. Estimates from multiple surveys could have been combined to gain precision (not included here).
4.2.3 Detection probabilities
The surveys of the sampling sites in the wild population were conducted with flights at 60 m above ground. Hence, we used the posterior distribution of eta[2] to inform priors for detection probabilities in the models fitted to these data. Transforming eta[2] to the probability scale gives a posterior median detection probability of 0.80 and a 95% credible interval ranging from 0.67 to 0.89 (see also Figure 2).
The posterior distribution of eta[2] is plotted in Figure 1 (black and red) and the corresponding posterior distribution for detection probability at the median site is plotted in Figure 2. As a prior for eta[2] in the analysis of the wild deer data from the sampling sites (blue lines), we used a normal distribution with the same mean and 1.5 times the standard deviation of the posterior samples from the model fitted to the enclosure data.
eta[2] (black), the normal distribution with the same mean and standard deviation as the posterior samples (red), and the normal distribution with the same mean and 1.5 times the standard deviation of the posterior samples (blue). The latter distribution (in blue) was used as a prior for eta[2] in the models fitted to the survey data from the wild populations.
Warning: package 'logitnorm' was built under R version 4.5.2
eta[2] transformed to the probability scale (i.e., detection probabiity at the median site). Colour codes are the same as in Figure 1.
The posterior distribution of sigma_p is shown in Figure 3 (note that there is a prior upper bound at 0.59; see Section 2.2 for prior specification). In Figure 4, this posterior distribution is transformed to the posterior distribution for the odds-ratio for detection at the 97.5 percentile site versus the 2.5 percentile site.
sigma_p. The range of the posterior samples is shown as a thick line in the density plot.
Including random effects on detection probability was deemed necessary by posterior predictive checks (Section 2.2), but there is obviously little information in the data to estimate heterogeneity in detection probability among sites. We also expect that the sigma_p parameter may be confounded with heterogeneity in actual density among sites (Section 2.2).
4.3 Wild deer
4.3.1 Overview of data
The full data are in the SiteData.rda files, with variables described in SiteDataInfo.html. Table 1 gives an overview of the data from the four study areas and two sampling periods (March and April).
| Sample area, month | Flight altitude (m) | Number of sites | Mean elevation (masl) | Mean area of sites (ha) | SD of site area | Total area (ha) | Number seen only by obs. 1 | Number seen only by obs. 2 | Number seen by both observers | Number seen in total | Simple estimate of density (ha\(^{-1}\)) |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Haugen, March | 60 | 43 | 533 | 0.65 | 0.07 | 28.07 | 1 | 2 | 11 | 14 | 0.511 |
| Haugen, April | 60 | 24 | 530 | 0.66 | 0.07 | 15.75 | 0 | 1 | 11 | 12 | 0.780 |
| Raa, March | 40 | 35 | 129 | 0.30 | 0.08 | 10.58 | 0 | 0 | 0 | 0 | 0.000 |
| Raa, April | 40 | 24 | 133 | 0.22 | 0.03 | 5.36 | 0 | 0 | 0 | 0 | 0.000 |
| Søre Bjørkum, March | 60 | 87 | 219 | 0.75 | 0.21 | 64.88 | 0 | 11 | 12 | 23 | 0.363 |
| Søre Bjørkum, April | 40 | 40 | 258 | 0.59 | 0.18 | 23.40 | 2 | 1 | 1 | 4 | 0.175 |
| Sprakehaug, March | 60 | 31 | 581 | 0.60 | 0.06 | 18.64 | 0 | 1 | 1 | 2 | 0.110 |
| Sprakehaug, April | 60 | 17 | 559 | 0.60 | 0.05 | 10.20 | 0 | 0 | 0 | 0 | 0.000 |
Most areas where surveyed by drone flights at 60 meter altitude above ground (pre-programmed based on a digital elevation model), but some areas where surveyed by flights at 40 meters above ground. Considerable variation in size of the sampling sites, among and within study areas, also occurred due to the steep and rugged terrain (Figure 5).
Deer were seen at 27 of the 301 sites, and the number of deer per site were distributed as follows
table(nimbleData$data$Y)
0 1 2 3 5 7
274 10 12 3 1 1 (I.e, there were 274 sites with zero deer seen, 10 sites with 1 deer, and so on).
The mean number of deer per site was 0.18 and the variance was 0.5 (the variance-mean ratio is 2.72). This shows that there is considerable overdispersion in number of deer per site relative to a Poisson distribution.
4.3.2 Simple estimates for MCMC staring values
Simple point estimate of detectability based on how many of those seen by observer 1 is seen by observer 2 and vice versa can be computed as follows:
colsumy <- apply(nimbleData$data$y, 3, sum, na.rm=TRUE)
p1hat <- colsumy[3]/(colsumy[2]+colsumy[3]) # = 0.692
p2hat <- colsumy[3]/(colsumy[1]+colsumy[3]) # = 0.923
(phat_simple <- 1 - (1-p1hat)*(1-p2hat)) both
0.9763314 This is the probability that at least one of the two observers registered a given deer in the picture (obs1 has probability 0.7 and obs2 has probability 0.9).
Simple density estimates based on number of deer seen divided by total area, corrected for detection probability <1 (number of deer per hectare (1/100 km\(^2\))):
(lambdahat <- (Y_per_sam/phat_simple)/area_per_sam) Haugen (April) Haugen (March) Raa (April)
0.7804950 0.5107558 0.0000000
Raa (March) Søre Bjørkum (April) Søre Bjørkum (March)
0.0000000 0.1750525 0.3631129
Sprakehaug (April) Sprakehaug (March)
0.0000000 0.1098906 4.3.3 WAIC model selection
Table 2 shows the WAIC values for the nine candidate models (see Section 1.2). The log-normal and gamma models give a significantly better fit to the data than the exponential models. The gamma models are slightly better (lower WAIC) than the log-normal models, but there is very little difference between the three gamma models. We will focus our inference on the gamma_0 (no covariate effect) and gamma_2 (a quadratic covariate effect) models.
| model | WAIC |
|---|---|
| exponential_0 | 122.84 |
| exponential_1 | 122.28 |
| exponential_2 | 121.48 |
| lognormal_0 | 113.68 |
| lognormal_1 | 113.64 |
| lognormal_2 | 113.79 |
| gamma_0 | 112.32 |
| gamma_1 | 112.40 |
| gamma_2 | 112.31 |
4.3.4 Summary of the WAIC-best models
Summary of the posterior samples from the gamma_0 model:
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu0[1] -1.0704 0.7483 0.0024943 0.003031
mu0[2] -1.7099 0.5662 0.0018874 0.002955
mu0[3] -10.1353 5.7643 0.0192144 0.179622
mu0[4] -10.8694 5.7181 0.0190603 0.208475
mu0[5] -2.8287 0.7452 0.0024839 0.003891
mu0[6] -2.0824 0.4167 0.0013890 0.002950
mu0[7] -10.2581 5.4858 0.0182860 0.128556
mu0[8] -3.3936 0.9820 0.0032734 0.005852
sigma 1.4934 0.1040 0.0003467 0.001906
eta1 1.0726 0.2385 0.0007949 0.007906
eta2 1.8295 0.2890 0.0009633 0.012143
sigma_p 0.3862 0.1547 0.0005156 0.011353
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu0[1] -2.4131 -1.5708 -1.111 -0.6261 0.5313
mu0[2] -2.8087 -2.0870 -1.717 -1.3437 -0.5791
mu0[3] -23.9575 -13.4812 -8.942 -5.6965 -2.3891
mu0[4] -24.4236 -14.0446 -9.717 -6.5018 -3.2194
mu0[5] -4.3346 -3.3146 -2.814 -2.3321 -1.3956
mu0[6] -2.9169 -2.3564 -2.076 -1.8024 -1.2802
mu0[7] -23.5083 -13.3408 -9.121 -6.0887 -2.8083
mu0[8] -5.4705 -4.0037 -3.346 -2.7326 -1.5899
sigma 1.2824 1.4254 1.495 1.5643 1.6926
eta1 0.6226 0.9086 1.069 1.2311 1.5524
eta2 1.3038 1.6235 1.820 2.0322 2.3944
sigma_p 0.0610 0.2789 0.422 0.5169 0.5835Summary of the posterior samples from the gamma_2 model:
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu0[1] -1.1067 0.7291 0.0024305 0.002946
mu0[2] -1.6812 0.5543 0.0018478 0.002838
mu0[3] -8.9378 5.3741 0.0179137 0.157768
mu0[4] -9.6603 5.5136 0.0183785 0.208530
mu0[5] -3.2823 0.8189 0.0027296 0.005188
mu0[6] -1.9036 0.5562 0.0018538 0.004154
mu0[7] -10.2412 5.3528 0.0178427 0.128269
mu0[8] -3.7547 1.0429 0.0034762 0.007175
beta[1] 0.9337 0.4436 0.0014788 0.004359
beta[2] -0.3039 0.2745 0.0009152 0.002705
sigma 1.4751 0.1045 0.0003482 0.001906
eta1 1.0757 0.2493 0.0008310 0.008396
eta2 1.7942 0.2830 0.0009432 0.010332
sigma_p 0.4192 0.1358 0.0004525 0.007925
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu0[1] -2.4300 -1.5932 -1.1431 -0.6668 0.4428
mu0[2] -2.7541 -2.0489 -1.6909 -1.3232 -0.5628
mu0[3] -21.4400 -12.2023 -7.9519 -4.8066 -1.3833
mu0[4] -22.5264 -12.9091 -8.5306 -5.4130 -2.1854
mu0[5] -4.9540 -3.8125 -3.2621 -2.7311 -1.7371
mu0[6] -2.9876 -2.2724 -1.9030 -1.5407 -0.8046
mu0[7] -23.4464 -13.2376 -9.1630 -6.1994 -2.9120
mu0[8] -5.9628 -4.4033 -3.7076 -3.0438 -1.8489
beta[1] 0.1254 0.6309 0.9123 1.2137 1.8693
beta[2] -0.8517 -0.4833 -0.3013 -0.1233 0.2322
sigma 1.2633 1.4065 1.4765 1.5458 1.6742
eta1 0.5925 0.9076 1.0724 1.2402 1.5750
eta2 1.2474 1.5995 1.7905 1.9880 2.3457
sigma_p 0.1146 0.3292 0.4558 0.5321 0.58484.3.5 Trace plots, posterior distribution plots, and convergence ststistics
Here, showing trace-plots and posterior density plots for the gamma_2 model only:
Gelman and Rubin (1992) convergence diagnostics of gamma_2 model:
Potential scale reduction factors:
Point est. Upper C.I.
mu0[1] 1.00 1.00
mu0[2] 1.00 1.00
mu0[3] 1.00 1.00
mu0[4] 1.01 1.02
mu0[5] 1.00 1.00
mu0[6] 1.00 1.00
mu0[7] 1.01 1.03
mu0[8] 1.00 1.00
beta[1] 1.00 1.00
beta[2] 1.00 1.00
sigma 1.00 1.00
eta1 1.01 1.03
eta2 1.00 1.01
sigma_p 1.13 1.38
Multivariate psrf
1.064.3.6 Posterior predictive checks
Posterior predictive checks with Bayesian p-values comparing simulated and observed values of \(Y\) and \(\mathbf{y}\) (see Section 1 and Section 2.2) are presented in Figure 6. Although the posterior median sum of squared Pearson residuals (ssPr) for \(Y\) simulated from the model are close that of the observed data (Bayesian p-value = 0.62), the distribution of the ssPr have a higher tail towards higher values for the simulated data compared to the observed data (Figure 6, top-right panel). This indicates that the upper credible limits from the gamma model may be somewhat conservative (note however that the exponential model (which has narrower credible intervals) has considerably less support by WAIC; see Table 8). There is no indication of serious lack of fit with respect to \(\mathbf{y}\) (Bayesian p-value = 0.20).
gamma2 model. Left column: Sum of the squared Pearson residuals for simulated data (y-axis) are plotted against the sum of the squared pearson residuals for the original data (x-axis) for every posterior sample (n = 90,000). The proportions of the points that are above the 1-to-1 lines in red are the Bayesian p-values, printed above the plots. Top panels represent the total counts by the two observers (\(Y\)), and bottom panels represent the counts broken down to ‘seen by only observer 1’, ‘seen by only observer 2’ and ‘seen by both’ (\(\mathbf{y}\) conditional on the observations of \(Y\)). The distributions of the descrepancy statistics are shown in the right column.
4.3.7 Predictions from best model
4.3.7.1 Deer density at the sample sites
While the main focus of the modelling was to estimate density where the sampled sites are assumed to be a random sample of a large number of sites of the same characteristics, we may also estimate deer density at the specific sites sampled (see Section 2.5). These estimates (posterior summaries) are presented in Table 3 from the gamma0 model, and in Table 4 from the gamma2 model. Uncertainty in these estimates are only due to imperfect detection.
gamma0 model) for deer densities at the sampled sites at each sample area and month. Units are number of deer per hectare (1/100 km\(^2\)). The three rightmost columns show the 50%, the 2.5% and the 97.5% quantiles (in the three cases where no deer were observed (at Raa and at Sprakehaug in April) the mean is higher than the 97.5% quantile due to the low median and high variance of the posterior distribution).
| Sample area and month | Mean (deer/ha) | Median (deer/ha) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 0.5184 | 0.4987 | 0.4987 | 0.6055 |
| Haugen (April) | 0.7909 | 0.7620 | 0.7620 | 0.8890 |
| Raa (March) | 0.0004 | 0.0000 | 0.0000 | 0.0000 |
| Raa (April) | 0.0008 | 0.0000 | 0.0000 | 0.0000 |
| Søre Bjørkum (March) | 0.3697 | 0.3699 | 0.3545 | 0.4162 |
| Søre Bjørkum (April) | 0.1783 | 0.1709 | 0.1709 | 0.2136 |
| Sprakehaug (March) | 0.1118 | 0.1073 | 0.1073 | 0.1609 |
| Sprakehaug (April) | 0.0004 | 0.0000 | 0.0000 | 0.0000 |
gamma2 model) for deer densities at the sampled sites at each sample area and month. Units are number of deer per hectare (1/100 km\(^2\)). The three rightmost columns show the 50%, the 2.5% and the 97.5% quantiles (in the three cases where no deer were observed (at Raa and at Sprakehaug in April) the mean is higher than the 97.5% quantile due to the low median and high variance of the posterior distribution).
| Sample area and month | Mean (deer/ha) | Median (deer/ha) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 0.5188 | 0.4987 | 0.4987 | 0.6055 |
| Haugen (April) | 0.7918 | 0.7620 | 0.7620 | 0.8890 |
| Raa (March) | 0.0004 | 0.0000 | 0.0000 | 0.0000 |
| Raa (April) | 0.0008 | 0.0000 | 0.0000 | 0.0000 |
| Søre Bjørkum (March) | 0.3704 | 0.3699 | 0.3545 | 0.4162 |
| Søre Bjørkum (April) | 0.1785 | 0.1709 | 0.1709 | 0.2136 |
| Sprakehaug (March) | 0.1120 | 0.1073 | 0.1073 | 0.1609 |
| Sprakehaug (April) | 0.0004 | 0.0000 | 0.0000 | 0.0000 |
4.3.7.2 Density estimates
Estimated densities as a quadratic function of distance from the field (gamma2 model) are shown in Figure 7, and estimates at the mean distance from the field are shown in Table 5. Table 6 shows the estimates from the gamma0 model.
gamma2 model). Blue lines are the posterior medians, and the 95% credible intervals are shown with the shaded area. Stippled vertical lines are plotted at the overall mean distance from the field (162 meters). Note the widely different ranges on the y-axis (logarithmic scale).
gamma2 model) for estimated deer densities at mean distance from the field (162 meters; stippled lines in Figure 7) at each of the study areas and months (when sites are considered to be a random sample of a large number of sites of the same characteristics). Units are number of deer per km\(^2\).
| Sample area and month | Mean (deer/km\(^2\)) | Median (deer/km\(^2\)) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 65.43 | 54.11 | 20.49 | 175.99 |
| Haugen (April) | 138.90 | 93.14 | 27.71 | 499.73 |
| Raa (March) | 3.80 | 0.06 | 0.00 | 33.60 |
| Raa (April) | 9.35 | 0.10 | 0.00 | 75.87 |
| Søre Bjørkum (March) | 52.31 | 43.57 | 16.08 | 139.92 |
| Søre Bjørkum (April) | 15.60 | 11.32 | 2.24 | 53.49 |
| Sprakehaug (March) | 11.70 | 7.26 | 0.79 | 47.65 |
| Sprakehaug (April) | 3.07 | 0.03 | 0.00 | 16.26 |
gamma0 model) for estimated deer densities at each of the study areas and months (when sites are considered to be a random sample of a large number of sites of the same characteristics). Units are number of deer per km\(^2\).
| Sample area and month | Mean (deer/km\(^2\)) | Median (deer/km\(^2\)) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 65.81 | 54.00 | 20.10 | 178.82 |
| Haugen (April) | 162.10 | 98.51 | 29.19 | 559.85 |
| Raa (March) | 1.37 | 0.02 | 0.00 | 12.24 |
| Raa (April) | 3.31 | 0.04 | 0.00 | 28.07 |
| Søre Bjørkum (March) | 41.40 | 37.84 | 18.21 | 85.90 |
| Søre Bjørkum (April) | 24.04 | 18.14 | 4.20 | 78.61 |
| Sprakehaug (March) | 16.42 | 10.70 | 1.32 | 63.98 |
| Sprakehaug (April) | 2.70 | 0.03 | 0.00 | 18.53 |
4.3.7.3 Process variation among sites
To interpret the parameter \(\sigma\) (sigma) which describes the variation in deer densities (\(\lambda\)) among sites, we plot quantiles of the fitted gamma-distribution in Figure 8. Since the variance is high compared to the mean, the median of the distribution is close to zero.
4.3.8 Overlap with prior
Degrees of overlap between posterior and prior distributions for the parameters in the gamma2 model are shown in Table 7. See Section 2 for prior justification. Note particularly that detection probability parameters eta1, eta2 and sigma_p had priors informed by posterior detection probabilities from the enclosure data (see Section 4.2.3), and hence there are relatively high overlap between posterior and prior distributions for these parameters. Posterior distributions for the ‘mu’ parameters (see Section 1.2.3) where no deer has been observed, “Raa” and “Sprakehaug (April)” (see section Section 4.3.1), also have an overlap >0.4.
Warning: package 'LaplacesDemon' was built under R version 4.5.2gamma2 model (multiply by 100 to get %).
| Parameter | Study area and month | Prior distribution | Overlap |
|---|---|---|---|
| beta[1] | N(0, sd = 2) | 0.349 | |
| beta[2] | N(0, sd = 2) | 0.262 | |
| eta1 | N(1.4, sd = 0.36) | 0.573 | |
| eta2 | N(1.4, sd = 0.36) | 0.533 | |
| mu0[1] | Haugen (April) | N(0, sd = 10) | 0.158 |
| mu0[2] | Haugen (March) | N(0, sd = 10) | 0.124 |
| mu0[3] | Raa (April) | N(0, sd = 10) | 0.471 |
| mu0[4] | Raa (March) | N(0, sd = 10) | 0.442 |
| mu0[5] | Søre Bjørkum (April) | N(0, sd = 10) | 0.165 |
| mu0[6] | Søre Bjørkum (March) | N(0, sd = 10) | 0.124 |
| mu0[7] | Sprakehaug (April) | N(0, sd = 10) | 0.412 |
| mu0[8] | Sprakehaug (March) | N(0, sd = 10) | 0.199 |
| sigma | Gamma(shape = 0.01, rate = 0.01) | 0.005 | |
| sigma_p | T[Gamma(shape = 1, rate = 0.05), 0, 0.59] | 0.684 |
4.3.9 Some sensitivity assessment
4.3.9.1 Choice of random distribution for \(\lambda\)
To assess the sensitivity of the model predictions to the choice of random distribution of the Poisson expectation, \(\lambda\), across sites, predictions from each of the nine candidate models at mean distance from field at Haugen in March are shown in Table 8. The exponential models give the most precise estimates, but have a substantially lower support by WAIC (Section 4.3.3). The log-normal models, which have marginally less support by WAIC than the gamma models, give substantially higher posterior mean and 97.5 percentile (upper 95% credible limit) than the gamma models (the posterior medians, and to a smaller degree the 2.5 percentile, are also somewhat higher). This is because the log-normal distribution has a heavier tail than the gamma-distribution with the same median and variance. The similarities in the WAIC values reflect that we have little information about the exact shape of the random distribution in the data (other than having evidence that the exponential distribution fits the data less well compared to the other two candidate distributions). In conclusion, we recognize that the upper 95% credible limits are highly model dependent, but have higher confidence in the lower 95% credible limits and the posterior medians.
| Model | WAIC | Mean (deer/km\(^2\)) | Median (deer/km\(^2\)) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|---|
| exponential_0 | 122.84 | 53.06 | 50.91 | 26.96 | 91.47 |
| exponential_1 | 122.28 | 49.16 | 47.20 | 24.61 | 85.37 |
| exponential_2 | 121.48 | 49.44 | 47.43 | 24.57 | 85.91 |
| lognormal_0 | 113.68 | 151.46 | 81.22 | 23.81 | 561.65 |
| lognormal_1 | 113.64 | 174.19 | 73.54 | 21.45 | 512.14 |
| lognormal_2 | 113.79 | 165.16 | 76.92 | 21.83 | 597.10 |
| gamma_0 | 112.32 | 65.81 | 54.00 | 20.10 | 178.82 |
| gamma_1 | 112.40 | 63.06 | 52.41 | 19.92 | 168.83 |
| gamma_2 | 112.31 | 65.43 | 54.11 | 20.49 | 175.99 |
4.3.9.2 Naive model without overdispersion
We did not consider a model without overdispersion as we did not think this would be biologically reasonable and simple preliminary analysis indicated that our data indeed was overdispersed relative to a Poisson distribution. However, since models without overdispersion (particularly for N-mixture models) are often used in the literature, we here also present results from a model without overdispersion and compare it with the WAIC-best model.
The summary from the model without overdispersion and a quadratic covariate effect is printed below. This model had significantly higher WAIC-value (128.07) than the models with overdispersion (see Table 8) and the posterior predictive goodness-of-fit assessment based on Bayesian p-values indicates a significant lack of fit with respect to the total counts \(Y\) (Bayesian p-value = 0.020; Figure 9). As a result, when ignoring overdispersion, credible intervals for predictions for the area represented by the samples are too narrow due to the overdispersion that has not been accounted for (Figure 10 and Table 9).
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu0[1] -0.3655 0.2987 0.0009955 0.0011138
mu0[2] -0.8007 0.2782 0.0009272 0.0010317
mu0[3] -8.9068 5.8388 0.0194627 0.0225293
mu0[4] -9.4107 5.7483 0.0191609 0.0224928
mu0[5] -2.1461 0.5503 0.0018344 0.0021084
mu0[6] -0.9655 0.2616 0.0008720 0.0012437
mu0[7] -9.8787 5.6200 0.0187334 0.0220341
mu0[8] -2.6677 0.8093 0.0026977 0.0031049
beta[1] 0.6573 0.2289 0.0007630 0.0014638
beta[2] -0.2070 0.1295 0.0004315 0.0008562
eta1 1.0738 0.2413 0.0008043 0.0079004
eta2 1.7969 0.2876 0.0009585 0.0114148
sigma_p 0.3952 0.1452 0.0004840 0.0102081
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu0[1] -0.99140 -0.5565 -0.3506 -0.1583 0.17633
mu0[2] -1.37444 -0.9813 -0.7896 -0.6073 -0.28961
mu0[3] -22.86689 -12.2703 -7.7293 -4.3924 -1.11586
mu0[4] -23.23184 -12.6819 -8.2386 -4.9687 -1.78125
mu0[5] -3.34135 -2.4865 -2.1031 -1.7598 -1.18702
mu0[6] -1.50144 -1.1368 -0.9572 -0.7861 -0.47632
mu0[7] -23.49693 -13.0395 -8.6939 -5.5558 -2.46794
mu0[8] -4.49800 -3.1426 -2.5798 -2.0991 -1.33107
beta[1] 0.23342 0.4987 0.6492 0.8057 1.13136
beta[2] -0.47441 -0.2914 -0.2016 -0.1172 0.03232
eta1 0.60702 0.9144 1.0683 1.2308 1.57088
eta2 1.21043 1.6079 1.8001 1.9937 2.34348
sigma_p 0.09682 0.2861 0.4268 0.5201 0.58382
| Sample area and month | Mean (deer/km\(^2\)) | Median (deer/km\(^2\)) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 46.63 | 45.40 | 25.30 | 74.86 |
| Haugen (April) | 72.45 | 70.43 | 37.11 | 119.28 |
| Raa (March) | 1.77 | 0.03 | 0.00 | 16.84 |
| Raa (April) | 3.39 | 0.04 | 0.00 | 32.76 |
| Søre Bjørkum (March) | 39.38 | 38.39 | 22.28 | 62.11 |
| Søre Bjørkum (April) | 13.45 | 12.21 | 3.54 | 30.51 |
| Sprakehaug (March) | 9.17 | 7.58 | 1.11 | 26.42 |
| Sprakehaug (April) | 0.90 | 0.02 | 0.00 | 8.48 |
4.3.9.3 Naive model without overdispersion and assuming perfect detection
We also fitted a model without extra-Poisson variation, and also assuming perfect detection (ignoring \(y\) and fitting the model to only the total count data \(Y\) while fixing \(p^*\) to 1). Density estimates from this model are presented in Figure 11 and Table 10. As expected, posterior predictive checks clearly show overdispersion of the data relative to this model (Figure 12), which leads to too narrow (overconfident) credible intervals.
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu0[1] -0.4040 0.2986 0.0009952 0.0010890
mu0[2] -0.8417 0.2778 0.0009259 0.0010206
mu0[3] -8.9473 5.8481 0.0194937 0.0225227
mu0[4] -9.4048 5.7377 0.0191257 0.0217053
mu0[5] -2.1928 0.5513 0.0018377 0.0020845
mu0[6] -1.0110 0.2602 0.0008672 0.0011816
mu0[7] -9.8938 5.5725 0.0185750 0.0215506
mu0[8] -2.7176 0.8074 0.0026912 0.0030743
beta[1] 0.6572 0.2308 0.0007694 0.0014556
beta[2] -0.2055 0.1295 0.0004318 0.0008412
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu0[1] -1.0284 -0.5948 -0.3906 -0.1975 0.14143
mu0[2] -1.4160 -1.0221 -0.8299 -0.6479 -0.33098
mu0[3] -22.9616 -12.2898 -7.7944 -4.4286 -1.14594
mu0[4] -23.3359 -12.6396 -8.2257 -4.9870 -1.78270
mu0[5] -3.3889 -2.5345 -2.1531 -1.8061 -1.22516
mu0[6] -1.5434 -1.1803 -1.0024 -0.8323 -0.52323
mu0[7] -23.2947 -13.0701 -8.7508 -5.5895 -2.50882
mu0[8] -4.5469 -3.1956 -2.6335 -2.1470 -1.38314
beta[1] 0.2287 0.4991 0.6490 0.8078 1.13414
beta[2] -0.4735 -0.2899 -0.2002 -0.1161 0.03404
| Sample area and month | Mean (deer/km\(^2\)) | Median (deer/km\(^2\)) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 44.75 | 43.61 | 24.27 | 71.82 |
| Haugen (April) | 69.72 | 67.67 | 35.76 | 115.19 |
| Raa (March) | 1.73 | 0.03 | 0.00 | 16.82 |
| Raa (April) | 3.29 | 0.04 | 0.00 | 31.79 |
| Søre Bjørkum (March) | 37.62 | 36.70 | 21.36 | 59.26 |
| Søre Bjørkum (April) | 12.85 | 11.61 | 3.37 | 29.37 |
| Sprakehaug (March) | 8.73 | 7.18 | 1.06 | 25.08 |
| Sprakehaug (April) | 0.86 | 0.02 | 0.00 | 8.14 |
4.3.9.4 Sensitivity to prior for \(\sigma_p\)
We used a somewhat informative prior for \(\sigma_p\) (sigma_p) (Section 2.2 and Section 4.3.8). To assess the sensitivity of the model results to this prior, we also fitted a model with a \(\text{Uniform}(0, 10)\) prior for \(\sigma_p\). The summary of this model is shown below and posterior summaries from this model are compared to the summaries from the original model in Figure 13.
Despite the substantially higher posterior values for \(\sigma_p\) with the wider prior, the posterior distributions for other parameters are little affected (Figure 13). There are no strong correlations between posterior samples of the the parameters for neither the original model (Table 11) nor the model with the wider prior for \(\sigma_p\) (Table 12). We conclude that our results are not sensitive to the choice of prior for \(\sigma_p\).
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu0[1] -1.0752 0.7193 0.0023975 0.003165
mu0[2] -1.6304 0.5546 0.0018486 0.003355
mu0[3] -9.3395 5.7063 0.0190211 0.174803
mu0[4] -9.4139 5.2895 0.0176318 0.189126
mu0[5] -3.1941 0.8153 0.0027176 0.005469
mu0[6] -1.8422 0.5553 0.0018510 0.004445
mu0[7] -10.2834 5.5520 0.0185066 0.141999
mu0[8] -3.6850 1.0454 0.0034848 0.007307
beta[1] 0.9313 0.4406 0.0014685 0.004492
beta[2] -0.2988 0.2738 0.0009127 0.002858
sigma 1.4603 0.1083 0.0003611 0.002109
eta1 1.3498 0.3235 0.0010783 0.003534
eta2 1.7470 0.3334 0.0011112 0.003486
sigma_p 2.8149 1.4725 0.0049083 0.113930
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu0[1] -2.3994 -1.5552 -1.1077 -0.6339 0.4513
mu0[2] -2.7124 -1.9991 -1.6345 -1.2706 -0.5166
mu0[3] -22.4682 -12.7631 -8.3043 -4.9496 -1.3547
mu0[4] -22.2840 -12.3709 -8.3578 -5.4134 -2.1778
mu0[5] -4.8709 -3.7195 -3.1704 -2.6449 -1.6685
mu0[6] -2.9164 -2.2108 -1.8455 -1.4797 -0.7379
mu0[7] -23.8318 -13.3143 -9.1482 -6.1587 -2.8319
mu0[8] -5.9109 -4.3375 -3.6297 -2.9737 -1.7799
beta[1] 0.1226 0.6309 0.9093 1.2109 1.8554
beta[2] -0.8511 -0.4758 -0.2944 -0.1170 0.2296
sigma 1.2386 1.3897 1.4632 1.5342 1.6641
eta1 0.7303 1.1286 1.3454 1.5671 1.9944
eta2 1.0858 1.5244 1.7497 1.9717 2.3952
sigma_p 0.9534 1.7604 2.4686 3.5025 6.7068
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu0[1] -1.1067 0.7291 0.0024305 0.002946
mu0[2] -1.6812 0.5543 0.0018478 0.002838
mu0[3] -8.9378 5.3741 0.0179137 0.157768
mu0[4] -9.6603 5.5136 0.0183785 0.208530
mu0[5] -3.2823 0.8189 0.0027296 0.005188
mu0[6] -1.9036 0.5562 0.0018538 0.004154
mu0[7] -10.2412 5.3528 0.0178427 0.128269
mu0[8] -3.7547 1.0429 0.0034762 0.007175
beta[1] 0.9337 0.4436 0.0014788 0.004359
beta[2] -0.3039 0.2745 0.0009152 0.002705
sigma 1.4751 0.1045 0.0003482 0.001906
eta1 1.0757 0.2493 0.0008310 0.008396
eta2 1.7942 0.2830 0.0009432 0.010332
sigma_p 0.4192 0.1358 0.0004525 0.007925
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu0[1] -2.4300 -1.5932 -1.1431 -0.6668 0.4428
mu0[2] -2.7541 -2.0489 -1.6909 -1.3232 -0.5628
mu0[3] -21.4400 -12.2023 -7.9519 -4.8066 -1.3833
mu0[4] -22.5264 -12.9091 -8.5306 -5.4130 -2.1854
mu0[5] -4.9540 -3.8125 -3.2621 -2.7311 -1.7371
mu0[6] -2.9876 -2.2724 -1.9030 -1.5407 -0.8046
mu0[7] -23.4464 -13.2376 -9.1630 -6.1994 -2.9120
mu0[8] -5.9628 -4.4033 -3.7076 -3.0438 -1.8489
beta[1] 0.1254 0.6309 0.9123 1.2137 1.8693
beta[2] -0.8517 -0.4833 -0.3013 -0.1233 0.2322
sigma 1.2633 1.4065 1.4765 1.5458 1.6742
eta1 0.5925 0.9076 1.0724 1.2402 1.5750
eta2 1.2474 1.5995 1.7905 1.9880 2.3457
sigma_p 0.1146 0.3292 0.4558 0.5321 0.5848
gamma2 model; orignal model in blue and model with wider \(\text{Uniform}(0, 10)\) prior for \(\sigma_p\) (sigma_p) in red. The lower credible bounds for \(\mu_0\) where no deer has been seen (mu0[3], mu0[4], and mu0[7]) are naturally sensitive to priors as it impossible to tell how close to zero the densities may be - on the arithmetic scale, the density estimates for these area-months are essentially zero.
| mu0[1] | mu0[2] | mu0[3] | mu0[4] | mu0[5] | mu0[6] | mu0[7] | mu0[8] | beta[1] | beta[2] | sigma | eta1 | eta2 | sigma_p | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mu0[1] | 1.000 | 0.021 | -0.008 | 0.004 | 0.027 | 0.048 | 0.000 | 0.025 | 0.001 | -0.043 | -0.090 | -0.018 | -0.010 | -0.005 |
| mu0[2] | 0.021 | 1.000 | 0.000 | 0.001 | 0.041 | 0.068 | 0.009 | 0.028 | 0.022 | -0.053 | -0.191 | -0.013 | -0.017 | -0.011 |
| mu0[3] | -0.008 | 0.000 | 1.000 | 0.044 | -0.003 | 0.010 | -0.013 | -0.001 | 0.047 | -0.041 | -0.002 | -0.021 | -0.025 | -0.004 |
| mu0[4] | 0.004 | 0.001 | 0.044 | 1.000 | 0.003 | 0.021 | -0.016 | 0.008 | 0.040 | -0.041 | -0.015 | -0.030 | -0.020 | -0.023 |
| mu0[5] | 0.027 | 0.041 | -0.003 | 0.003 | 1.000 | 0.153 | 0.024 | 0.128 | -0.237 | -0.068 | -0.175 | -0.008 | -0.014 | 0.006 |
| mu0[6] | 0.048 | 0.068 | 0.010 | 0.021 | 0.153 | 1.000 | 0.016 | 0.119 | 0.215 | -0.549 | -0.185 | -0.009 | -0.015 | -0.001 |
| mu0[7] | 0.000 | 0.009 | -0.013 | -0.016 | 0.024 | 0.016 | 1.000 | 0.009 | -0.031 | 0.002 | -0.038 | 0.035 | -0.019 | 0.007 |
| mu0[8] | 0.025 | 0.028 | -0.001 | 0.008 | 0.128 | 0.119 | 0.009 | 1.000 | -0.149 | -0.072 | -0.132 | -0.005 | -0.020 | 0.008 |
| beta[1] | 0.001 | 0.022 | 0.047 | 0.040 | -0.237 | 0.215 | -0.031 | -0.149 | 1.000 | -0.682 | 0.064 | 0.006 | 0.002 | -0.006 |
| beta[2] | -0.043 | -0.053 | -0.041 | -0.041 | -0.068 | -0.549 | 0.002 | -0.072 | -0.682 | 1.000 | -0.010 | -0.008 | 0.000 | 0.007 |
| sigma | -0.090 | -0.191 | -0.002 | -0.015 | -0.175 | -0.185 | -0.038 | -0.132 | 0.064 | -0.010 | 1.000 | -0.024 | -0.003 | -0.022 |
| eta1 | -0.018 | -0.013 | -0.021 | -0.030 | -0.008 | -0.009 | 0.035 | -0.005 | 0.006 | -0.008 | -0.024 | 1.000 | 0.118 | 0.198 |
| eta2 | -0.010 | -0.017 | -0.025 | -0.020 | -0.014 | -0.015 | -0.019 | -0.020 | 0.002 | 0.000 | -0.003 | 0.118 | 1.000 | 0.038 |
| sigma_p | -0.005 | -0.011 | -0.004 | -0.023 | 0.006 | -0.001 | 0.007 | 0.008 | -0.006 | 0.007 | -0.022 | 0.198 | 0.038 | 1.000 |
sigma_p)
| mu0[1] | mu0[2] | mu0[3] | mu0[4] | mu0[5] | mu0[6] | mu0[7] | mu0[8] | beta[1] | beta[2] | sigma | eta1 | eta2 | sigma_p | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mu0[1] | 1.000 | 0.031 | -0.005 | 0.004 | 0.035 | 0.052 | 0.001 | 0.031 | -0.012 | -0.033 | -0.113 | -0.004 | -0.027 | 0.047 |
| mu0[2] | 0.031 | 1.000 | 0.003 | 0.001 | 0.047 | 0.079 | 0.001 | 0.035 | 0.020 | -0.052 | -0.208 | -0.006 | -0.037 | 0.083 |
| mu0[3] | -0.005 | 0.003 | 1.000 | -0.034 | -0.007 | 0.019 | 0.006 | 0.001 | 0.038 | -0.033 | -0.015 | 0.001 | 0.005 | 0.012 |
| mu0[4] | 0.004 | 0.001 | -0.034 | 1.000 | -0.015 | 0.027 | 0.013 | -0.018 | 0.062 | -0.059 | -0.003 | 0.025 | -0.002 | 0.026 |
| mu0[5] | 0.035 | 0.047 | -0.007 | -0.015 | 1.000 | 0.149 | 0.010 | 0.124 | -0.226 | -0.069 | -0.172 | -0.015 | -0.028 | 0.053 |
| mu0[6] | 0.052 | 0.079 | 0.019 | 0.027 | 0.149 | 1.000 | 0.018 | 0.117 | 0.213 | -0.543 | -0.201 | -0.014 | -0.035 | 0.091 |
| mu0[7] | 0.001 | 0.001 | 0.006 | 0.013 | 0.010 | 0.018 | 1.000 | 0.018 | -0.023 | -0.008 | -0.004 | -0.005 | -0.003 | -0.011 |
| mu0[8] | 0.031 | 0.035 | 0.001 | -0.018 | 0.124 | 0.117 | 0.018 | 1.000 | -0.148 | -0.069 | -0.145 | -0.007 | -0.034 | 0.043 |
| beta[1] | -0.012 | 0.020 | 0.038 | 0.062 | -0.226 | 0.213 | -0.023 | -0.148 | 1.000 | -0.685 | 0.075 | 0.000 | 0.007 | 0.004 |
| beta[2] | -0.033 | -0.052 | -0.033 | -0.059 | -0.069 | -0.543 | -0.008 | -0.069 | -0.685 | 1.000 | -0.021 | 0.007 | -0.007 | -0.003 |
| sigma | -0.113 | -0.208 | -0.015 | -0.003 | -0.172 | -0.201 | -0.004 | -0.145 | 0.075 | -0.021 | 1.000 | 0.014 | 0.050 | -0.116 |
| eta1 | -0.004 | -0.006 | 0.001 | 0.025 | -0.015 | -0.014 | -0.005 | -0.007 | 0.000 | 0.007 | 0.014 | 1.000 | 0.013 | 0.162 |
| eta2 | -0.027 | -0.037 | 0.005 | -0.002 | -0.028 | -0.035 | -0.003 | -0.034 | 0.007 | -0.007 | 0.050 | 0.013 | 1.000 | -0.190 |
| sigma_p | 0.047 | 0.083 | 0.012 | 0.026 | 0.053 | 0.091 | -0.011 | 0.043 | 0.004 | -0.003 | -0.116 | 0.162 | -0.190 | 1.000 |
4.3.10 Proof of concept - Identifiability assessment by simulations
To check that there are no inherent identifiability problems with the model, we fitted the model to two different sets of simulated data. In the first simulation, we set parameters to the posterior mean of the original fit and just boosted sample size to 500 sites within each of the area-month combinations. Then, in the second simulation, we also reduced the parameter \(\sigma\) (sigma) describing the among site variation in density to 5% of the estimated value (posterior median), in addition to boosting sample size. The same prior distributions as in the original fit was used in both simulations.
4.3.10.1 Boosting sample size
Posterior summaries from the model fitted to the boosted data are shown below, overlap between prior and posterior distributions are shown in Table 13, predicted densities are plotted in Figure 14, and predicted densities at mean distance from the field are shown in Table 14.
With this much data (8 areas with 500 sites each), we obtain substantially higher precision of the estimates, as well as lower overlap between posterior and prior distributions. This shows that that there are no inherent identifiability problems in the model. One would perhaps expect higher precision (narrower credible intervals) when having 8 \(\times\) 500 sites, but the high variation in number of deer per sites contributes to much uncertainty, and since the average area of sites is just under 0.6 ha (0.6% of a km\(^2\)), the expected counts per site is quite small. In the next simulation, we also reduced the random variation among sites within areas.
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
beta[1] 0.8995 0.10428 0.0003476 0.0009246
beta[2] -0.2924 0.07542 0.0002514 0.0007335
eta1 0.9857 0.10396 0.0003465 0.0071195
eta2 1.7378 0.12180 0.0004060 0.0096278
mu0[1] -1.1607 0.16759 0.0005586 0.0012902
mu0[2] -1.6801 0.17255 0.0005752 0.0013398
mu0[3] -7.1290 1.20691 0.0040230 0.0385446
mu0[4] -12.8709 5.00279 0.0166760 0.6756462
mu0[5] -3.2925 0.24855 0.0008285 0.0017874
mu0[6] -1.6387 0.17670 0.0005890 0.0012403
mu0[7] -13.0874 5.21458 0.0173819 0.7356025
mu0[8] -4.6047 0.39887 0.0013296 0.0041903
sigma 1.4497 0.03246 0.0001082 0.0006676
sigma_p 0.3354 0.15154 0.0005051 0.0188382
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
beta[1] 0.69761 0.8285 0.8983 0.9705 1.1054
beta[2] -0.43831 -0.3435 -0.2930 -0.2418 -0.1431
eta1 0.79104 0.9115 0.9848 1.0558 1.1953
eta2 1.51099 1.6499 1.7371 1.8246 1.9718
mu0[1] -1.48827 -1.2731 -1.1612 -1.0476 -0.8329
mu0[2] -2.01966 -1.7967 -1.6791 -1.5639 -1.3436
mu0[3] -9.99584 -7.7731 -6.9418 -6.2783 -5.3060
mu0[4] -25.56872 -15.4319 -11.7917 -9.1060 -6.4719
mu0[5] -3.78775 -3.4578 -3.2898 -3.1240 -2.8168
mu0[6] -1.98488 -1.7576 -1.6392 -1.5192 -1.2908
mu0[7] -26.77264 -15.6869 -11.8302 -9.1332 -6.5592
mu0[8] -5.43604 -4.8623 -4.5871 -4.3278 -3.8701
sigma 1.38511 1.4282 1.4499 1.4714 1.5132
sigma_p 0.07493 0.2129 0.3384 0.4721 0.5771
| Parameter | Study area and month | Prior distribution | Overlap |
|---|---|---|---|
| beta[1] | N(0, sd = 2) | 0.106 | |
| beta[2] | N(0, sd = 2) | 0.087 | |
| eta1 | N(1.4, sd = 0.36) | 0.275 | |
| eta2 | N(1.4, sd = 0.36) | 0.363 | |
| mu0[1] | Haugen (April) | N(0, sd = 10) | 0.043 |
| mu0[2] | Haugen (March) | N(0, sd = 10) | 0.043 |
| mu0[3] | Raa (April) | N(0, sd = 10) | 0.182 |
| mu0[4] | Raa (March) | N(0, sd = 10) | 0.286 |
| mu0[5] | Søre Bjørkum (April) | N(0, sd = 10) | 0.058 |
| mu0[6] | Søre Bjørkum (March) | N(0, sd = 10) | 0.044 |
| mu0[7] | Sprakehaug (April) | N(0, sd = 10) | 0.282 |
| mu0[8] | Sprakehaug (March) | N(0, sd = 10) | 0.084 |
| sigma | Gamma(shape = 0.01, rate = 0.01) | 0.002 | |
| sigma_p | T[Gamma(shape = 1, rate = 0.05), 0, 0.59] | 0.867 |
| Sample area and month | Mean (deer/km\(^2\)) | Median (deer/km\(^2\)) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 54.06 | 53.31 | 38.67 | 73.97 |
| Haugen (April) | 90.82 | 89.50 | 65.75 | 123.53 |
| Raa (March) | 0.05 | 0.00 | 0.00 | 0.45 |
| Raa (April) | 0.39 | 0.28 | 0.01 | 1.42 |
| Søre Bjørkum (March) | 56.41 | 55.45 | 39.89 | 78.22 |
| Søre Bjørkum (April) | 10.95 | 10.66 | 6.54 | 17.00 |
| Sprakehaug (March) | 3.09 | 2.91 | 1.25 | 5.95 |
| Sprakehaug (April) | 0.05 | 0.00 | 0.00 | 0.41 |
4.3.10.2 Boosting sample size and reducing \(\sigma\)
To investigate the impact of variation in density within study areas, we reduced \(\sigma\) to 5% of the estimated value while keeping the mean densities the same. This entails changing \(\sigma\) to \(\sigma'=\alpha\sigma\) and \(\mu_{0,j}\) to \(\mu_{0,j}'= \mu_{0,j} + \sigma^2(1-\alpha^2)/2\) where \(\alpha = 0.05\). The resulting variation in density among sites is shown in Figure 15.
Posterior summaries from the model fitted the simulated data are shown below, overlap between prior and posterior distributions are shown in Table 15, predicted densities are plotted in Figure 16, and predicted densities at mean distance from the field are shown in Table 16.
Iterations = 1:30000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 30000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
beta[1] 1.08400 0.07952 0.0002651 0.012521
beta[2] -0.41645 0.05512 0.0001837 0.012074
eta1 1.06999 0.09161 0.0003054 0.004536
eta2 1.84020 0.11817 0.0003939 0.007088
mu0[1] -0.03231 0.06629 0.0002210 0.008476
mu0[2] -0.61594 0.07560 0.0002520 0.009180
mu0[3] -7.04693 1.84428 0.0061476 0.072798
mu0[4] -7.89821 1.00514 0.0033505 0.092203
mu0[5] -2.18414 0.16818 0.0005606 0.018155
mu0[6] -0.65359 0.07894 0.0002631 0.009018
mu0[7] -7.20157 1.27744 0.0042581 0.120449
mu0[8] -2.68511 0.21937 0.0007312 0.015500
sigma 0.04003 0.04843 0.0001614 0.008764
sigma_p 0.37152 0.13475 0.0004492 0.023799
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
beta[1] 0.913492 1.008271 1.11365 1.14087 1.21291
beta[2] -0.514581 -0.444248 -0.43159 -0.37142 -0.31456
eta1 0.890593 1.007851 1.06926 1.13148 1.25187
eta2 1.615310 1.758540 1.83732 1.91603 2.08174
mu0[1] -0.167719 -0.078850 -0.03162 0.01811 0.08219
mu0[2] -0.769082 -0.657829 -0.61693 -0.58667 -0.43182
mu0[3] -10.196259 -9.329650 -6.33193 -5.30688 -4.92005
mu0[4] -9.230805 -8.992867 -7.86135 -7.05486 -6.28998
mu0[5] -2.431293 -2.355133 -2.14915 -2.07305 -1.92029
mu0[6] -0.799318 -0.701028 -0.66041 -0.61560 -0.48165
mu0[7] -9.238339 -8.388753 -7.35461 -5.96905 -4.63586
mu0[8] -3.094282 -2.890581 -2.60810 -2.55403 -2.26910
sigma 0.002172 0.006771 0.01813 0.06216 0.17384
sigma_p 0.114691 0.259565 0.38750 0.48131 0.57872
| Parameter | Study area and month | Prior distribution | Overlap |
|---|---|---|---|
| beta[1] | N(0, sd = 2) | 0.067 | |
| beta[2] | N(0, sd = 2) | 0.046 | |
| eta1 | N(1.4, sd = 0.36) | 0.302 | |
| eta2 | N(1.4, sd = 0.36) | 0.279 | |
| mu0[1] | Haugen (April) | N(0, sd = 10) | 0.017 |
| mu0[2] | Haugen (March) | N(0, sd = 10) | 0.021 |
| mu0[3] | Raa (April) | N(0, sd = 10) | 0.119 |
| mu0[4] | Raa (March) | N(0, sd = 10) | 0.089 |
| mu0[5] | Søre Bjørkum (April) | N(0, sd = 10) | 0.049 |
| mu0[6] | Søre Bjørkum (March) | N(0, sd = 10) | 0.020 |
| mu0[7] | Sprakehaug (April) | N(0, sd = 10) | 0.148 |
| mu0[8] | Sprakehaug (March) | N(0, sd = 10) | 0.045 |
| sigma | Gamma(shape = 0.01, rate = 0.01) | 0.050 | |
| sigma_p | T[Gamma(shape = 1, rate = 0.05), 0, 0.59] | 0.786 |
| Sample area and month | Mean (deer/km\(^2\)) | Median (deer/km\(^2\)) | Lower 95% C.I. | Upper 95% C.I. |
|---|---|---|---|---|
| Haugen (March) | 54.27 | 54.05 | 46.51 | 64.98 |
| Haugen (April) | 97.22 | 97.14 | 84.80 | 108.78 |
| Raa (March) | 0.06 | 0.04 | 0.01 | 0.19 |
| Raa (April) | 0.24 | 0.18 | 0.00 | 0.73 |
| Søre Bjørkum (March) | 52.29 | 51.71 | 45.06 | 61.88 |
| Søre Bjørkum (April) | 11.44 | 11.68 | 8.80 | 14.74 |
| Sprakehaug (March) | 7.00 | 7.38 | 4.53 | 10.39 |
| Sprakehaug (April) | 0.17 | 0.06 | 0.01 | 0.97 |
References
Banerjee, Sudipto. 2024. “Finite Population Survey Sampling: An Unapologetic Bayesian Perspective.” Sankhya A 86 (1): 95–124. https://doi.org/10.1007/s13171-024-00348-8.
Brack, Ismael V., Andreas Kindel, and Luiz Flamarion B. Oliveira. 2018. “Detection Errors in Wildlife Abundance Estimates from Unmanned Aerial Systems ( UAS ) Surveys: Synthesis, Solutions, and Challenges.” Edited by Kylie Scales. Methods in Ecology and Evolution 9 (8): 1864–73. https://doi.org/10.1111/2041-210X.13026.
Conn, Paul B., Devin S. Johnson, Perry J. Williams, Sharon R. Melin, and Mevin B. Hooten. 2018. “A Guide to Bayesian Model Checking for Ecologists.” Ecological Monographs 88 (4): 526–42. https://doi.org/10.1002/ecm.1314.
de Valpine, Perry, C Paciorek, D Turek, N Michaud, C Anderson-Bergman, F Obermeyer, C Wehrhahn Cortes, A Rodrìguez, D Temple Lang, and S Paganin. 2024. NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling. Zenodo. https://doi.org/10.5281/zenodo.10602714.
de Valpine, Perry, Daniel Turek, Christopher Paciorek, Cliff Anderson-Bergman, Duncan Temple Lang, and Ras Bodik. 2017. “Programming with Models: Writing Statistical Algorithms for General Model Structures with NIMBLE.” Journal of Computational and Graphical Statistics 26 (2): 403–13. https://doi.org/10.1080/10618600.2016.1172487.
de Valpine, P, C Paciorek, D Turek, N Michaud, C Anderson-Bergman, F Obermeyer, C Wehrhahn Cortes, A Rodrìguez, D Temple Lang, and S Paganin. 2024. NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling. Zenodo. https://doi.org/10.5281/zenodo.10602714.
Gelman, Andrew, Carlin, John B., Stern, Hal S., and Rubin, Donald B. 2004. Bayesian Data Analysis. 2nd ed. Texts in Statistical Science. Chapman & Hall/CRC.
Gelman, Andrew, and Donald B. Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Stat. Sci. 7 (4): 457–511. https://doi.org/doi:10.1214/ss/1177011136.
Nichols, James D., James E. Hines, John R. Sauer, Frederick W. Fallon, Jane E. Fallon, and Patricia J. Heglund. 2000. “A Double-Observer Approach for Estimating Detection Probability and Abundance from Point Counts.” The Auk 117 (2): 393–408. https://doi.org/10.1093/auk/117.2.393.
Plummer, M, N Best, K Cowles, and K Vines. 2006. “CODA: Convergence Diagnosis and Output Analysis for MCMC.” R News 6: 7–11.
R Core Team. 2024. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Watanabe, Sumio. 2010. “Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory.”
Williams, B. K., J. D. Nichols, and M. J. Conroy. 2002. Analysis and Management of Animal Populations. Academic Press.