Factors influencing the success of aerial rabies vaccination of foxes

Sylvatic rabies has been eradicated from most of Central Europe, but cases still occur in the Balkans. Oral rabies vaccination of foxes is an effective method for controlling the disease. The aim of this study was to evaluate the success of aerial vaccination campaigns conducted in Montenegro by identifying ecological, environmental and climatic factors that influenced the prevalence of antibodies to the rabies vaccine. To monitor the bait uptake and the serological responses to vaccination, foxes were shot by hunters. Of 175 shot foxes, 142 foxes (81.1%) had consumed baits. Of these only a total of 81 (57.0%) tested positive for rabies vaccine antibodies, possibly, due to the delayed uptake of bait in which the rabies vaccine was already inactivated. We found that low vaccination responses were associated with high fox density and bait delivery in open areas. In high fox density habitat, bait uptake might be delayed as other food and prey options for foxes are abundant. Similarly, delayed bait uptake probably occurred in open areas as such areas are less frequently used by foxes. The findings of this study suggest that efficacy of oral rabies vaccination by aerial delivery is associated with landscape features.


Supporting information 1 Testing residual spatial correlation
Let y i denote the number of foxes with rabies antibodies out of n i at location x i . We define the Person's residuals asẐ wherep i is the predicted probability that a fox has rabies antibodies.
We define the empirical semivariogram of the residualsẐ i aŝ where N (u) = {(i, j) : ||x i − x j || = u}, i.e. the set of all pairs of data-points such that their Euclidean distance is u, and |N (u)| is the number of pairs within the set. Intuitively, in the presence of residual spatial correlation in the data, we would expect the squared difference (Ẑ h −Ẑ k ) 2 to be smaller, on average, at smaller distances u, as a result of the stronger correlation betweenẐ h andẐ k . Conversely, in the absence of spatial correlation, the empirical semivariogram would only show random fluctuations around a constant value, sinceẐ h andẐ k would vary independently of each other at any given distance u.
To reliably distinguish the first case from the latter, we use the following Monte Carlo strategy to simulate empirical semivariograms under the assumption of spatial independence.
1. Permute the locations x i while holding fix the order of theẐ i .
2. Compute the empirical semivariogram (1) using the permuted set of locations.
3. Repeat 1 and 2 a large enough number of times, say B.
4. Based on the resulting B empirical semivariograms, compute the 95% confidence intervals for each of the pre-specified distance bins.
If the empirical semivariogram lies outside outside the 95% tolerance bandwidth of step 4, we then interpret this as evidence of residual spatial correlation. If, instead, the empirical semivariogram lies within the 95% bandwidth, we then conclude that the data do not show evidence of residual spatial correlation.
2 Assessing the impact of home ranges above 1 km 2 on parameter estimation and predictions In this section, we shall usex to denote the location where the fox consumed the bait.
Sincex may not coincide with the location x, where the fox was shot, we modelx as a inhomogeneous Poisson process (Diggle, 2014) with intensity λ(x) given by the fox density at locationx. Let h denote the home range of a red fox, measured in km 2 ; since each grid cell has an area of 1 km 2 , we then consider h ∈ {1, 2, 3, 4} (Cavallini & Lovari, 1994). For a given home range h, we then define all possible combinations with adjacent grid cells to x, such that the total area is h. Figure 1 shows an example for h = 2 km 2 . In general, we have that the total number of possible habitat configurations, for a given a home range h, is The log-likelihood function for the regression coefficients β, based on the binary outcome y i , for i = 1, . . . , n, indicating the presence/absence of rabies antibodies is then given by where C j is the j-th habitat configuration out of N h , and l i (β|x) = y i η(x) − log 1 + e η(x) , i = 1, . . . , n with linear predictor η(x) = d(x) β and explanatory variables d(x). To estimate β, we then maximize (2) using numerical optimization. Standard errors are computed by taking the square root of the diagonal elements of the inverse of the negative Hessian matrix at the maximum likelihood estimate. Figure 2 and Figure 3 show the point estimates with associated 95% confidence intervals for 2011 and 2012, respectively. Table 1 and Table 2 report the percentage relative change in the estimates of the regression coefficients with respect to those at a home range of 1 km 2 , for 2011 and 2012, respectively. For example, standardized percentage of grass has a relative change of 10.68% at a home range of 4 km 2 with respect to a home range of 1 km 2 .