Accounting for spatial habitat and management boundaries when estimating forest bird population distribution and density: inferences from a soap film smoother

Study area

As previously described in Camp et al. (2020), Hakalau (19°51′N, 155°18′W) was established in 1985 and is actively managed to preserve native forest birds, rainforest plants and their habitats. The 15,390-ha montane forest is dominated by native ‘ōhi’a lehua Metrosideros polymorpha and koa Acacia koa with a mixture of native and non-native understory plants. Temperature averages 15 °C with annual variation <5 °C, and annual precipitation averages 2,500 mm with a maximum of 6,100 mm (Juvik & Juvik, 1998). Hakalau is roughly divided into three strata: afforested-pasture, open- and closed-forest (Fig. 1; see Camp et al., 2016 for a detailed map). The afforested-pasture stratum extends from the western edge of the open-forest stratum to the western edge of the Hakalau boundary, covering an area of 1,271 ha, with elevation between 1,650 and 2,000 m above sea level (asl). The afforested-pasture was intensely grazed through the mid-1980s and is being gradually reforested through outplanting koa. The open-forest stratum occurs at an elevation between 1,400 and 1,920 m asl and encompasses an area of 3,061 ha. Previously heavily grazed, the open-forest stratum has undergone natural regeneration since the removal of cattle in 1988 (Maxfield, 1998). The closed-forest stratum extends down slope from the eastern edge of the open-canopy stratum to the eastern edge of the Hakalau boundary, with elevation between 1,400 and 1,700 m asl. The closed-forest stratum was least modified by grazing and surveys commenced in 1999 covering an area of 1,541 ha.

Bird surveys

Following the creation of Hakalau, in 1987 refuge staff initiated an annual abundance monitoring program by establishing 350 point samplers (hereafter, points) on 14 transect lines following a systematic random design spanning the upper elevations of Hakalau. U.S. Fish and Wildlife, Big Island National Wildlife Refuge Complex, coordinated and conducted forest bird monitoring on the Hakalau Forest Unit. The distance among points is approximately 150 m with transects being either roughly 500 or 1,000 m apart (Fig. 3). For this analysis we selected a single survey from the ‘Ākepa time series based on broad sampling of the study area and sufficient numbers of detections. In 2002, 289 points were sampled using PTDS methods within the forest stratum where 276 ‘Ākepa were detected at 121 points. Observers recorded the species, detection type (heard, seen, or both), and horizontal distance from the station center point to individual birds detected during an 8-min count. Observers also recorded cloud cover, rain intensity, wind strength on the Beaufort scale, and gust strength on the Beaufort scale.

Study boundary

To exemplify extrapolating to management unit or political boundaries we set the western edge of the boundary to the coordinates of the pasture-forest boundary, while the north, east and south boundaries were squared off (hereafter, referred to as the “soap boundary”) using the function locator (Becker, Chambers & Wilks, 1988) (Fig. 3). In addition to defining the soap boundary, soap film smoothing also required a priori specification of the number and location of knots inside the boundary. We defined a set of knots on a regular grid with locations every 730 m east and 670 m north across the study area (Fig. 3).

Stage 1 modelling detectability

Using the R package Distance (Miller, 2017, Miller et al., 2019) we fitted half-normal and hazard-rate key detection functions without either adjustment terms or covariates (Fig. 2). Preliminary analysis revealed key detection functions with adjustment terms were not strictly monotonic, and there were insufficient detections per factor to model covariates. Data were truncated at a distance w where the estimated detection probability (using a preliminary detection function model) was about 0.1. Model selection for the detection function used AIC (Burnham & Anderson, 2002) and model fit was evaluated with a Cramér-von Mises test (Buckland et al., 2015).

Stage 2 design-based density estimation

For comparative purposes we estimated the ‘Ākepa abundance using conventional distance sampling methods (Fig. 2; Buckland et al., 2015). We combined the estimated detection probability with counts to estimate densities as birds ha⁻¹ using package Distance. Abundance was estimated as density times the area inside the boundary (5,671.8 ha). Variance and 95% CI were calculated using analytic methods (Buckland et al., 2015).

Stage 2 model-based density estimation

Comparison between estimates

Estimated coefficients cannot be directly compared between most spline models. Instead, we predicted density and uncertainty estimates and mapped the distribution of ‘Ākepa densities within the soap boundary on a grid with cells 200 m × 200 m. We also evaluated the model predictions at even distances along the boundary. Estimates of abundance and 95% confidence intervals (CI) were computed from the posterior distribution of the fitted soap model following methods described by Wood (2017) and detailed in Camp et al. (2020). A total of 10,000 replicate parameter value sets were drawn from the posterior distribution of the model coefficients, $\hat{\mathit{\beta }}$, assuming a multivariate normal distribution. We estimated abundance for each of these parameter samples and calculated the mean of the replicate sets, SE as the standard deviation of the replicates, and 95% CIs were computed from the 2.5th and 97.5th quantiles. The detection function uncertainty was not propagated with the uncertainty from either model-based approach as the detection probability does not contain any spatial variability (modelled without covariates since there were no detection covariates). All analyses were conducted using the R statistical language and code is provided in Appendix B. We computed the change in uncertainty as the percentage change in interval widths (CIW) between soap and TPRS, and soap and Distance abundance estimates expressed as $\left(\frac{\mathrm{C}\mathrm{I}{\mathrm{W}}_{TPRS}}{\mathrm{C}\mathrm{I}{\mathrm{W}}_{\mathrm{s}\mathrm{o}\mathrm{a}\mathrm{p}}}-1\right)\times 100\mathrm{\%}$ and $\left(\frac{\mathrm{C}\mathrm{I}{\mathrm{W}}_{Distance}}{\mathrm{C}\mathrm{I}{\mathrm{W}}_{\mathrm{s}\mathrm{o}\mathrm{a}\mathrm{p}}}-1\right)\times 100\mathrm{\%}$, respectively.

Design-based density estimation

Using a preliminary detection function model, truncation was set at 58 m. The AIC for the hazard-rate detection function model was 11 units smaller than that of the half-normal model (Table 1). The Cramér-von Mises test was non-significant at the $\alpha =0.05$ level indicating that the detection function provided a satisfactory fit to the distance histogram (Table 1). Inspection of diagnostic plots also indicated that the model adequately fit the data (Appendix C). The shoulder of the detection function extends out to 30 m before decaying rapidly. The estimated detection probability of a bird that was within 58 m of a point was 0.631 (SE = 0.035) and the effective area surveyed per point was 6,668.6 m².

Soap film smoother

Our soap film smoother explained 54.6% of the deviance in the data using 2.98 of 18 boundary degrees of freedom and 14.45 of 108 interior degrees of freedom, and the index to ensure adequate complexity was 1.1 (desired value should be near one for each term). The negative binomial dispersion parameter was 10.55 indicating that the counts were over-dispersed.

The main advantage of using the soap smooth is that there was no leakage across the pasture-forest boundary (Fig. 4, left panel). The density surface maps showed an ‘Ākepa density hotspot in the southern portion of the domain that extended north-east to a second hotspot in the central portion of the domain (Fig. 4, left panel). Both soap smooth hotspots occurred interior to the boundary, even though the parametrization allowed for fitted values to persist adjacent to and including the boundary. Densities throughout the northern portion of the domain were zero or near zero. The uncertainty estimates portrayed the same pattern with large SEs predicted in the southern portion that extended north-east to the central portion of the domain (Fig. 4, right panel). SEs in the northern portion of the domain were near zero while SEs were moderate adjacent to the south boundary and throughout a swath along the south-eastern boundary.

On the boundary the largest densities predicted by the soap smoother occurred in the south-western corner of Hakalau (labelled zero in Fig. 5) with densities of about 1.25 birds ha⁻¹. Densities declined smoothly to near zero by and along the northern boundary (occurring between 300 and 600 in Fig. 5) and then increased progressively along the eastern boundary until again reaching maximum densities at the south-western corner.

TPRS smoother

The TPRS smoother had an estimated negative binomial dispersion parameter of 12.40 and explained 55.1% of the deviance. There was a clear pattern of lower densities in the north than in the south part of Hakalau. The EDF on the Easting and Northing smoothing term was greater than zero, significant and the function was nonlinear (Supporting Online Information Appendix A Table 2; Supporting Online Information Appendix A Fig. 7).

There were two density hotspots and two SE hotspots in the density surface maps predicted from the TPRS model (Fig. 4, bottom left and right panels respectively). The hotspot in the south portion was relatively large spanning north-east and extending to the south and west boundaries. The second hotspot occurred in the central portion of the domain and extended to the east boundary. The contours of the TPRS model continued beyond the soap boundary, which resulted in very convoluted densities along the boundary (Fig. 5).

Comparison of smooths predictions

Both smoothers had visually similar cold and hotspots in the domain with similar predicted densities and SEs (Fig. 4). Both smoothers predicted zero or near zero densities in the northern part of the domain. The extent of the coldspot was much larger for the soap smooth than for the TPRS smooth extending from the central hotspot to the north boundary, while the TPRS smooth predicted approximately zero densities to about half the area of that predicted by the soap smooth. The large hotspot in the southern part of the study area was similar in shape and extent between the two smoothers. The local hotspot in the central part of the study area was more symmetrical for the soap smooth than for the TPRS smooth which was more triangularly-shaped. For the TPRS, both hotspots extended to the boundary. In contrast, both soap smooth hotspots occurred interior to the boundary, even though the parametrization allowed for fitted values to persist right to and including the boundary.

Greater differences between the two smoothers were observed in the fitted SE than the fitted density estimates (Fig. 4). The TPRS and soap both had two hotspots, but the SE in the central part of the domain was adjacent to and extended along the boundary. The equivalent soap SE hotspot was more centrally located between the east and west boundaries and the SE estimates were generally smaller. The TPRS and soap smooth global maximum SEs were in the southern part of the domain and the soap hotspot portrayed a patchwork of varying SEs while the TPRS SEs were relatively uniform across the entirety of the hotspot with less cell-to-cell variability.

The soap-film basis argument controlling the boundary smooth eliminated the extreme rough densities of the TPRS smooth (Fig. 5). On the boundary the largest densities predicted by the soap smooth occurred in the south-west corner of Hakalau (labelled zero in Fig. 5). Densities declined smoothly to near zero by and along the northern boundary (occurring between 300 and 600 in Fig. 5) and then increased progressively along the eastern boundary until again reaching maximum densities of nearly 1.25 birds ha⁻¹ at the south-west corner. Predicted densities from the TPRS smooth were more rough with peaks at points 80, near 700 and at 950 along the boundary (Fig. 5). There were minor peaks at about points 100, 200 and 800. Each of these peaks occurred where contours intersected the boundary (compare Fig. 1, right panel, and Fig. 5). In between each peak the TPRS smooth predicted small to very small densities before increasing rapidly to the next peak. Similar to the soap smooth, the TPRS smooth predicted near zero densities along the northern boundary, between points 400 and 600. The peaks along the boundary indicates leakage by the TPRS smooth.

Abundance estimates

The soap film smooth estimate was 6,518–9,630 ‘Ākepa with a median point estimate of 7,743 ± 794 (Table 2). The TPRS smooth produced similar estimates to the soap estimates and the abundances of the two approaches differed by more than 100 birds with a TPRS median point estimate of 7,619 ± 673 birds. Estimates of the 2002 ‘Ākepa abundance using design-based, conventional distance sampling analyses were similar to the soap smooth estimates. The 95% CIs of each approach bracketed the mean abundance estimates of the other two approaches. The soap CIW was 16.5% wider than the TPRS CIW, and was 0.8% wider than the Distance CIW.

Statistical methods

We used a soap film smoother to model spatial densities of ‘Ākepa from PTDS count data. The two dimensional soap film smoother comprises two separate but linked bases; one for the boundary and one for the film itself. In this case the soap basis arguments provide a good approach for estimating the boundary and interior surface splines. This approach of estimating densities along the boundary allowed us to answer the question posed in the Introduction: do ‘Ākepa occur at relatively high densities adjacent to the pasture-forest boundary? Densities at the boundary varied along the pasture-forest boundary as well as within the domain (see Figs. 4 and 5). As seen in Fig. 5, the predicted values on the boundary were smallest in the north-western corner and increased along the western boundary toward the south-western corner. Thus, in the southern portion of Hakalau ‘Ākepa occur at relatively high densities near to the edge of the pasture-forest boundary while in the northern portion ‘Ākepa densities were near zero inside the domain and on the boundary, as well as outside the domain in the pasture stratum.

The soap smooth model estimating the boundary reflects that ‘Ākepa densities may extend beyond the soap boundary. As seen in the ‘Ākepa densities, estimating density along the boundary is more realistic since it is not certain that the boundary value will be zero. Along the eastern and southern boundaries the coefficient indicates that the population extends beyond the boundary into the adjacent forest habitat and the boundary that poses no real physical barrier in the ‘Ākepa distribution. Fixing the boundary to zero instead of estimating it would force the smooth of the interior surface to decrease to zero at the boundary. This approach may be appropriate for island populations where birds are restricted to suitable habitat that is located within an inhospitable matrix. Examples of this include bird populations on Pacific islands such as Aguiguan in the Mariana Islands (Amidon et al., 2014), Tau in American Samoa (Judge et al., 2013) and Nihoa in the Hawaiian Islands (Gorresen et al., 2016) where suitable forest bird habitat occurs across these islands and extends to the coastal non-forest habitat or high-tide water line. For more complex systems, the approach applied here can be used to address and overcome challenges of holes (i.e., unsuitable habitat such as lakes, urban areas and agricultural fields for forest obligate species) within the domain and to model the domain boundary.

The choice of whether to start with a TPRS or soap film smoother should take into consideration features of the study area and data that may affect predicted densities. Features such as a lake or a change in habitat along the study area may result in leakage. A general approach may be to start with a TPRS smoother and if leakage is identified switch to a soap film smoother. The peaks in predicted densities along the boundary indicate leakage by the TPRS smooth, and thus the soap film smoother provides more biologically realistic ‘Ākepa densities. The most common ecological applications of soap film smoothers are in marine ecosystems, with species such as fish (Augustin et al., 2013, Polansky et al., 2018), cetaceans (Dellabianca et al., 2016) and seabirds (Grecian et al., 2016). Alternative methods that respect complex boundaries include finite element L-splines (Ramsay, 2002), geodesic low rank thin plate spline methods (Wang & Ranalli, 2007), multi-dimensionally-scaled Duchon splines (Miller & Wood, 2014), complex region spatial smoother (CReSS; Scott-Hayward et al., 2014), and barrier models (Bakka et al., 2019). Each method has several advantages and disadvantages, and comparisons conducted by Miller & Wood (2014) and Scott-Hayward et al. (2014) indicate that the CReSS and soap film methods perform better than the other methods. Moreover, soap film smoothers fit easily within the GAM framework and are readily modelled in R (R Core Team, 2017) package mgcv (Wood, 2021). Limitations of GAMs are that they are restricted to additive models, and while interactions can be captured it can be difficult to construct interactions involving many variables. The main issues with using the soap film are that it can be hard to parameterize, requiring specifying the boundary, checking model convergence, and interpreting the results. The main advantages of using GAMs are that they are easy to interpret, the flexible predictor can identify complex data patterns, and regularisation of smoothing functions avoids overfitting (Wood, 2017). Moreover, the software package mgcv (Wood, 2021) facilitates modelling smoothers including model selection and model checking, making the analyses accessible to a broad user community.

Biological findings

The hotspot in the southern part of the domain coincides with the locally abundant density estimates identified by Scott et al. (1986) at 200+ birds km^--2. Similarly to the Scott et al. (1986) predictions, our results indicate that ‘Ākepa densities decrease outside the hotspot. ‘Ākepa extend outside the domain on three sides of the soap boundary but are restricted to forest habitats above 1,500 m elevation. On the northern and eastern sides of the domain this is within several 100 m of the boundaries. The ‘Ākepa distribution to the south of Hakalau continues along the 1,500 m elevation contour before terminating several kilometres south of the study area (Judge et al., 2018). The current sampling frame is centred to track the core, high density proportion of the population, and although the extent of the survey is limited, the soap smoother based density estimates appear to be a good approximation of the ‘Ākepa distribution and abundance in the region (Scott et al., 1986; Lepson & Freed, 1997).

Abundances of between 6,300 to 9,700 ‘Ākepa in the extended forest stratum boundary seems realistic for the 2002 survey in this high density ‘Ākepa population (Gorresen et al., 2009; Camp et al., 2010, 2016). Simulation studies of design-based distance sampling have shown that method to be unbiased when the critical assumptions are met or possess low bias when assumptions fail (Buckland et al., 2001, 2015). This is not necessarily the case with model-based approaches and model mis-specification can lead to bias. The standard procedure in this case is to check the residual structure of the selected model. The residual plots were qualitatively similar and the desired refitted model EDF was approximately zero for each term. While we are not able to assess bias in the soap film smoother estimates, it is insightful to compare these estimates to those produced using the design-based approach assuming that the latter approach is unbiased. In this case it appears that both the soap and TPRS smoothers are unbiased. The 95% CIs among the three estimators overlapped the abundance estimate produced by the other models (see Table 2). The coefficient of variation of abundance was very small at about 10% for the Distance, soap, and TPRS models.

Along the southern portion of the pasture-forest boundary ‘Ākepa occur in relatively large densities (see Fig. 5). The combination of suitable habitat and high densities provides the potential for ‘Ākepa to colonize the pasture-stratum naturally. It was not until 2011 that ‘Ākepa had been detected in the afforested pasture (Paxton et al., 2018). Surveys subsequent to 2012 indicate that ‘Ākepa continue to use the koa in the pasture stratum (Steve J. Kendall, 2020, personal observations). These occurrences are likely diel migrations to locate foraging sources (Lepson & Freed, 1997); however, there is no evidence of ‘Ākepa breeding outside of established and mature forested habitats. As afforestation of the pasture stratum progresses ‘Ākepa appear poised to colonize this restored forest.