Animal counts.

Animals were counted using attempted total counts. The attempted total count method is widely used for counting animals in African savannas. For example, attempted total counts from a vehicle are used to count animals every two months in the Lake Nakuru National Park adjoining the Nakuru Wildlife Conservancy since 1970 [25–27,52] and in Nairobi National Park of Kenya since 1960 [6,30,53–57]. Attempted aerial total counts have also been used to count animals in the Kruger National Park in South Africa from 1977 to 1997 [44,58,59], Masai Mara National Reserve of Kenya since 1962 [60] and the Serengeti National Park since 1958 [61]

Even though it cannot provide estimates of the accuracy of population size estimates or detectability of animals, the attempted total count method is able to provide estimates for even rare animal species, which are of great conservation concern but which are hard to count using sampling techniques. For example, after the total counting procedure was replaced with aerial distance sampling in the Kruger National Park in 1998, it became difficult to obtain reliable estimates of trends for several species, such as warthog, wildebeest and waterbuck, even with a coverage intensity exceeding 22% of the park area [62]. Similarly, in the Nakuru Wildlife Conservancy, [11] used distance sampling to count potential cheetah prey species along 13 foot transects, covering a total of 55.3 km of the Conservancy. As expected, sightings of rare species were insufficient to model detection functions. Population size estimates for several species could only be obtained by pooling together several species of “similar size” and that “provide similar visual cues”. Even so, population estimates could not be calculated at all for some species, such as warthog, which had insufficient sightings to model a detection function and are hard to group with others [11].

The Nakuru Wildlife Forum has conducted biannual (twice a year) game counts throughout the conservancy area since March 1996, with the exceptions of 2001 and 2013 when only one count was conducted and 2015 for which the dry season count was not yet available. Counts on all properties are conducted on the same day and at the same time, as much as possible, to reduce the possibility of double counting animals. The different properties are subdivided into blocks using road networks, vegetation and area to be covered. Teams of counters consisting of at least one driver, one spotter and one recorder drive through each block and record total counts of all medium to large animal species encountered. The counts are undertaken jointly by volunteers from the Kenya Wildlife Service (KWS) and the individual conservancies who know the landscape and animal distribution in the conservancy well. The counts start at 0600 am and end at 1100 am when animals are most active and likely to be encountered in the open landscape. During the counts the conservancy is partitioned into the southern and northern sections, separated by Lake Naivasha and dense human settlements that limit interchange of wildlife between the two sections within the counting period. The counts are done using 4-wheel drive vehicles, on foot where the terrain is not accessible by vehicle due to steep or rugged topography, or vegetation cover and using boats on the shores of Lake Naivasha. The maximum vehicle speed during the surveys is 20 km/h. The counts are done with the aid of binoculars. On average, the entire counting exercise takes about three days to complete. The first day is used to train the new counters on the counting method, and for volunteers to meet the landowners and coordinate the survey logistics. The landowners also use the first day to inform volunteers about any changes in the conservancies or ranches since the last census. Actual counting is then carried out during the next two consecutive days. The expansive northern ranches, including Delamere, Kikopey and Marula ranches and the surrounding areas are covered in one day by about 50 counters and 15 vehicles. The Kedong ranches and other southern properties, including Hell’s Gate National Park, Longonot National Park and Oserian, are also covered in one day by about 50 counters and 20 vehicles. The foot counts are only carried out in the densely vegetated areas that are inaccessible by car, located predominantly in Oserian and Delamere properties and along the Lake Naivasha shoreline in Marula; all together totaling less than 20% of the entire Conservancy area. A boat is used for about one hour only to count buffaloes occurring inside the reeds as well as hippos along the shoreline. The total counts cover virtually all of the NWC as shown in Fig 1b.

Volunteers from a conservancy, park or ranch always count in the same counting blocks but volunteers from outside, including those from KWS, change occasionally as the individual volunteers vary over time. The same tracks are followed as much as possible during the surveys. The counts are done mainly in September-October and March-May of each year. These periods were selected to represent the dry and the wet season conditions, respectively. Availability of green plant biomass for wildlife is highest in March-May and lowest in September-October. The block counts are later summed to obtain total counts for each property and for the whole conservancy to track changes in animal numbers over time. We used all the counts carried out during 1996–2015 for this study. Because the same methods are applied in each count, any potential biases in the counts should be comparable, so that changes in density over time should reliably reflect overall changes. Except for the national parks and the Oserian sanctuary, all the properties constituting the conservancy support both wildlife and livestock, with Marula and Soysambu supporting the largest livestock herds. However, the livestock counts were not made available to us.

Rainfall.

We analysed temporal variation in rainfall from 1964 to 2011 to provide a broad historical context for evaluating rainfall variation during our study period spanning 1996–2015. Monthly rainfall was averaged across all the seven recording stations to account for spatial variation. The resulting average monthly rainfall was used to calculate the total annual (January-December), wet (March-August) and dry (September-January) season rainfall components. It was also used to compute the total quarterly rainfall in terms of the early (March-May) and late (June-August) wet season and the early dry (September-November) and late dry (December-February) season components.

We used an unobserved components model (UCM), a special case of the linear Gaussian state space model [63], to decompose the time series of the annual, wet season and dry season rainfall series (r_t) into the trend (μ_t), cyclical (ψ_t), seasonal (γ_t) and irregular (ϵ_t) components as follows (1) where β_j are regression coefficients, x_jt are regression variables with fixed effects and (ϵ_t) are independent and identically (i.i.d) normally distributed errors or disturbances with zero mean and variance . In other words, we assume ϵ_t to be a Gaussian white noise process. The trend, seasonal, cyclical and irregular components are assumed to be statistically independent of each other.

We assume a random walk (RW) model for the time trend, which is equivalent to assuming that the trend (μ_t) remains roughly constant through time. The RW trend model is (2) where . Note that if then μ_t = constant.

We further assume a stochastic cycle (ψ_t) with a fixed period (p > 2) and damping factor (ρ) but a time-varying amplitude and phase specified by (3) where 0 ≤ ρ ≤ 1, λ = 2 × π/p is the angular frequency of the cycle, ϑ_t and are independent Gaussian disturbances with zero mean and variance and 0 < λ < π. Values of ρ, p and are estimated from the data alongside the other model parameters. The damping factor ρ governs the stationarity properties of the random sequence ψ_t such that ψ_t has a stationary distribution with mean zero and variance if ρ < 1 but is nonstationary if ρ = 1. We allowed and tested for up to three cycles in the annual, wet season and dry season rainfall components.

In addition to the random walk model (2) we modelled the trend component using a locally linear time trend incorporating both a level and slope and described as (4) where the disturbance variances and are assumed to be independent. The UCM models (1) and (4), without the seasonal and regression components, were fitted by the diffuse Kalman filtering and smoothing algorithm in the SAS UCM procedure [64].

Modeling temporal trends in animal density.

Because the number of properties participating in the counts, and thus the total area surveyed varied over time, the total number of animals of each species counted in each census was divided by the total area covered by the census to obtain density (number / km²) estimates for each species. We therefore analyzed time trend in the density of each species and the biomass density of all herbivores, primates and carnivores. Both the wet and dry season counts were analyzed together because there are no large scale migratory or seasonal dispersal movements, nor substantial seasonal variation in visibility, and hence detectability of the species, due to a preponderance of short to medium grasses.

We modelled trends in the densities of all the wildlife species simultaneously using a flexible semi-parametric model that accommodates irregularly spaced, missing and non-normally distributed counts with many zeroes. The model assumed that the animal counts follow a negative binomial distribution in which the variance is a quadratic function of the mean.

The negative binomial distribution (NB) of animal counts (Y) can be given by (5)

The mean of Y is given by μ = E(Y) and its variance by Var(Y) = μ(1 + μ/k). k is a shape or dispersion parameter and quantifies the amount of overdispersion [65].

Let z ≡ (x, u) be a covariate vector with x = (x₁, …, x_p, x_p+1, …, x_q), a 1 × (p + q) covariate vector and u a continuous independent variable. Further, denote the mean of Y with u(z) = E(Y|z). Then the negative binomial probability distribution in Eq (5) can be recast in exponential format by (6)

This formulation shows that the NB model belongs to the exponential family of distributions for known k.

Next, we model the temporal trends in the animal counts by letting u be an unspecified smooth function of time and the fixed effects in x to be linear. Assuming a log link function, because the canonical link log[μ(z)/{μ(z) + k}] implied by Eq (6) can be problematic to work with as it is always negative, the expected counts given the covariates is then specified as (7) where β = (β₁, …, β_p, β_p+1, …, β_p+q)^T is a 1 × (p + q)-parameter vector for the p + q -covariate vector x = (x₁, …, x_p, x_p+1, …, x_p+q), s₁(·) and s₂(·) are unspecified smooth functions that capture the effects of u₁ and u₂. The logarithm of the total area counted in each census, log(area), is an offset used to calculate animal density. The offset has a coefficient equal to unity by construction. The unspecified smooth functions s₁(·) and s₂(·) are approximated by penalized B-spline basis functions as (8) where b_l and c_j are the penalized cubic B-spline coefficients to be estimated. The computational details and mathematical properties of B-splines can be found in [66].

We now use the time series of the 44 species monitored in NWC from 1996 to 2015 to clarify how the smooth functions s₁(⋅) and s₂(·) in Eq (8) are approximated by a mixed model. The number of all observations is , where n_i is the number of all observations for the i-th of the total of 44 species counted. Furthermore, p = 44 is the number of coefficients to be estimated for the species main effect and q = 44 is the number of coefficients representing the species × time interaction effect. Accordingly, the full design matrix of fixed effects x in Eq (7) has dimension n × (p + q) = [1400 × (44 + 44)] whereas the vector of fixed effect parameters β has dimension 1 × (p + q) = 1 × 88.

If we let represent the (n × K) matrix of B-splines of degree d and D_r the (K − r × K) matrix of r-th order difference penalty, then the (n × K − r) matrix used to fit the mixed model specified by 7 and 8 equals (9)

The total number of B-spline knots used to specify equals the number m of equally spaced interior knots plus d knots placed at the starting date and max{1, d} knots placed at the ending date of the censuses. The total number of columns in the B-spline basis is thus K = m +d + 1.

The number of variables (columns) representing the penalized B-spline (called P-spline) random effect of time trend common to all the species is given by . Here, the number of interior knots is m = 20, the degree of the B-spline basis is d = 3 and the order of differences of the spline coefficients is r = 3. The number of variables (columns) representing the random P-spline effect for the time trend specific to each species (species interaction × time) is then given by . It follows that for these data the random effect design matrix in Eq (8) has dimension whereas in Eq (8) has dimension . The full design matrix of random effects therefore has dimension . The vector of parameters of random effects in Eq (8) has dimension whereas has dimension . The full vector of parameters of random effects has dimension . The random effects whereas .

The full model specified by Eqs (7), (8) and (9) can therefore be described as a multivariate, semiparametric generalized linear mixed model. The model is semiparametric because it consists of parametric and non-parametric components. The fixed part of the model is parametric and consists of the species main effect, representing species-specific overall densities, and species-by-time trend interaction, representing species-specific trends. The random part of the model, specifying the non-parametric smooth functions, consists of two continuous random spline effects, each specified by a penalized spline variance-covariance structure. The first random spline effect fits a penalized cubic (d = 3) B-spline [67] with a third-order difference penalty (r = 3) to the random spline coefficients common to all the species and thus models the time trend common to all the species. The second random spline effect similarly fits a penalized cubic B-spline with random coefficients specific to each species and hence models the time trend unique to each species. Each of the two random spline effects had m = 20 equally spaced interior knots placed on the running date of censuses (March 1996, …, May 2015) plus 3 evenly spaced exterior knots placed both at the start date and end date of the censuses, for a total of 26 knots. The model yields estimates of three variance components, corresponding to the random spline time trend common to all species , random spline effect for the time trend specific to each species and the scale parameter for the negative binomial distribution k.

Fitting the model exploits the idea that spline smoothing and mixed modeling address equivalent minimization problems and produce the same solutions. However, the solutions for the spline coefficients within the mixed modeling framework are solutions of random effects, not fixed effects, as is the case for solutions of spline coefficients in the classical framework. Standard errors of the predicted counts thus account for the variation in spline coefficients associated with treating the coefficients as random effects in mixed models [68]. A distinct advantage of formulating spline smoothing as a mixed model is that the smoothing parameter is selected ‘automatically’ because it is a function of the covariance parameter estimates produced by the mixed model. The trend model was fitted via maximization of residual penalized quasi-likelihood, called pseudo-likelihood by [69], in the SAS GLIMMIX procedure. The SAS codes (SAS Version 9.4, SAS /STAT version 14.1) used to fit the rainfall and wildlife trend models are provided in Supplementary Materials SM1.

To establish if population density for each species changed significantly between 1996 and 2015 we used constructed spline effects [64] to compare the expected population size for 1996 to that for 2015. For species whose numbers initially increased to a peak and then declined, the expected peak population size was additionally compared with the expected population sizes in March 1996 and May 2015.

Relationship between rainfall and animal biomass and density.

We calculated the total biomass (kg/km²) of all the 32 herbivore species (S1 Table) for each year and related the biomass to 1- to 10-year moving averages of the annual, seasonal and quarterly rainfall components. The 10-year cut-off window was selected to coincide with the period of the longest cycle established for the annual and wet season rainfall components. The cumulative past rainfall component most strongly correlated with the aggregate herbivore biomass or the densities of each of the 30 most abundant species, chosen to minimize potential influences of stochastic noise due to small population sizes of some species, was selected using Pearson correlation coefficients and the corrected Akaike Information Criterion [70].

The role of conservancies in community-based wildlife conservation and management.

The conservation efforts in the NWC spearheaded by various agencies under the auspices of the Nakuru Wildlife Forum had evidently positive effects on the 32 wildlife species whose numbers were either stable or increasing but not on the 12 species whose numbers were declining. It is unlikely that the population increases were due solely to the two small national parks within the NWC, which made relatively minor contributions to the total wildlife population in the NWC throughout 1996–2015. In fact, the number of individuals of each species counted outside the protected National Parks as a percent of the total number counted both inside and outside the parks from 1996 to 2015 averaged 92.8% ±13.2%, emphasizing the importance of the private and communal lands to the success of conservation efforts. This success reinforces the view that private and communal lands play a pivotal role in wildlife conservation in the NWC. The private and communal lands are not only important to conservation but are becoming increasingly more so, not only in NWC but throughout most of Kenya’s pastoral lands that cover about 88% of the terrestrial land surface of Kenya and support between 65 and 70% of all wildlife in Kenya [1,2]. By comparison, officially protected wildlife reserves and parks cover a mere 10–12% of Kenya [77,79].

Conservation in the NWC is due to the collective action of the landholders working in partnerships with governmental, international and non-governmental organizations under the umbrella of the Nakuru Wildlife Forum [21,24]. Similar collaborative arrangements are also promoting wildlife conservation on private and communal lands in Laikipia, Masai Mara, Amboseli and other parts of Kenya. Given the increase in numbers of many wildlife species in the NWC, it is tempting to conclude that collaborative conservation had a generally positive impact on wildlife population performance in the NWC. However, the decline in some wildlife species and other growing rafts of challenges call for carefully re-examining and improving the conservation and management strategies used by members of the NWF to minimize further wildlife losses.

A fundamental question is the extent to which such collaborative management protects and promotes healthy wildlife populations and ecosystems. The involvement of local landowner associations bringing together private, communal and public landowners in the conservation and management of biodiversity outside protected areas is critical because most biodiversity lies outside the protected areas in East Africa. The community associations are gaining prominence as government resources available for conservation become increasingly inadequate in the wake of a growing raft of challenges. The land owner associations are well placed to drive collaborative and locally-based conservation initiatives given their detailed knowledge, skills and strong vested interests in the success of the conservancies. The associations enable individual landowners to integrate their land parcels into broader landscape and regional biodiversity conservation frameworks. By choosing to conserve and protect wildlife, landowners greatly expand the area available for conservation and critical ecosystem services in rural landscapes, buffer protected areas from surrounding human impacts and complement the limited capacity and skills of government agencies [79].

Conservancies and landowner associations that conserve and manage them are emerging as the centerpiece of conservation practice, policies and strategies in pastoral lands of Kenya and should be encouraged and supported as much as possible. The varied circumstances of the different conservancies and the landowners allow for pluralistic and locally-adaptive solutions to conservation challenges prevailing in specific localities. The future success of the conservancies will most strongly hinge on whether and how well wildlife conservation evolves to become a major component of sustainable livelihoods in the rangelands. This is likely to happen if landowners are able to derive meaningful income from wildlife conservation. Currently, land owners in the NWC, like elsewhere in Kenya, are faced with limited variety in sources of income from wildlife after KWS withdrew, in 2002, the wildlife use rights granted to the Nakuru-Naivasha region in 1992. Yet, landowners or landholders incur considerable costs by allowing wildlife to compete with their livestock without expecting any compensation from the government that owns wildlife, thereby enabling populations of many wildlife species to flourish outside the small protected areas. Moreover, most properties in NWC lack natural water sources in the dry season and pump water from rivers, boreholes and the lake for both wildlife and livestock. This is another substantial cost borne entirely by the landholders. Additionally, the landowners provide security for wildlife on their properties. Currently, a number of properties in NWC derive some benefits from wildlife by operating various ecotourism facilities.

Kenya’s Parliament passed a new Wildlife Conservation and Management Act in 2013 that provides new legislative, administrative and policy frameworks and principles for devolving rights and responsibilities for wildlife conservation, management, utilization, habitat recovery and restoration efforts and compensation for losses caused by wildlife in Kenya [80]. The Act advocates principles, policies and practices that promote co-existence between people and wildlife and expand the range of benefits of wildlife conservation and offset wildlife-related losses to those living with wildlife. By aiming to turn wildlife from a liability to an asset for individuals or communities living with wildlife and giving legal backing to organizations such as NWF, the new Act is stimulating more effective engagement in, and support for, collaborative wildlife conservation and management in human-dominated and modified landscapes of Kenya by private and communal land owners.