10.5061/DRYAD.S1RN8PK4W
Terraube, Julien
0000-0002-7060-3988
University of Helsinki
Van doninck, Jasper
University of Turku
Helle, Pekka
Natural Resources Institute Finland
Cabeza, Mar
University of Helsinki
Data from: Assessing the effectiveness of a national protected area
network in maintaining carnivore populations
Dryad
dataset
2020
protected areas
biodiversity loss
Maj ja Tor Nesslingin Säätiö
http://dx.doi.org/10.13039/501100004157
201700201
2020-05-01T00:00:00Z
en
5329573 bytes
4
CC0 1.0 Universal (CC0 1.0) Public Domain Dedication
Protected areas (PAs) are essential to prevent further biodiversity loss
yet their effectiveness varies largely with governance and external
threats. Although methodological advances have permitted assessments of PA
effectiveness in mitigating deforestation, we still lack similar studies
for the impact of PAs on wildlife populations. Here we demonstrate the
application ofuse an innovative combination of matching methods and
hurdle-mixed models with a large-scale and long-term dataset of
unprecedented coverage for Finland’s large carnivore species. We show that
the national PA network , at the national level, PAs does not support
higher densities than non-protected habitat for 3 of the 4 species
investigated. For the brown bear, PAs appear to have lower densities than
non-protected areas. For some species, PA effects interact with region or
time, i.e. wolverine densities decreased inside PAs over the study period
and lynx densities increased inside eastern PAs. Although we show that
matching approaches could and should be applied to wildlife population
data, Wwe support their application of matching methods in combination of
additional analytical frameworks for deeper understanding of conservation
impacts on wildlife populations. These methodological advances are crucial
for improving PA targets and extremely timely for preparing ambitious PA
targets a post-2020 global framework for biodiversity.
Wildlife population time series data Track count data collected within the
Finnish Wildlife Triangle Scheme were used as indices of lynx, wolf,
wolverine and bear densities. The scheme is a long term, large-scale
survey which provides yearly estimates of the distribution and relative
abundance of game species. The methods are described in detail by Linden
et al.55. A wildlife unit is a permanent line transect route of 12 km (4 x
3 km). Unit locations represent different forest habitats in proportion to
local occurrence. The monitoring operation is carried out twice a year.
The minimum number of monitoring experts is three for the late-summer
count and one for the winter count. Each winter, 800–900 units distributed
throughout Finland are surveyed by local hunters who count the number of
fresh snow tracks crossing the transect. The count is either done 1–2 days
after a snow fall or 24 h after a pre-check where all old tracks are
marked. For each species, the data are compiled by the Finnish Game and
Fisheries Research Institute as the number of crossings per 24 h per 10
km. Thus, the unit used in statistical modelling was the average number of
crossings per 24 h per 10 km for lynx wolf and wolverine, per triangle and
for each year. As brown bears hibernate during winter counts, density
indices for this species were collected during late summer counts and also
focused on counting tracks crossing the transect. From the initial 2171
sampling units, a total of 1805 units were included in subsequent analyses
(we discarded transects in mixed agricultural-forest landscapes), 1195
were located in non-protected sites while 610 were located inside
protected areas. Explanatory variable selection and extraction The set of
matching covariates extracted for each wildlife unit represented
ecological characteristics considered most likely to be important
determinants of large carnivore occurrence in European ecosystems based on
empirical observation56. These covariates can be regarded as two groups of
possible influences: biophysical context (e.g. latitude) and human impacts
(e.g. distance to closest settlements). Several variables suggested by the
literature to be important in determining PA effectiveness, e.g. PA
management, were available for a restricted amount of sites and therefore
were not included in the matching analysis. Description of the matching
covariates and the rationale for their inclusion are found in the
Supplementary Table S8. Spatial data were analyzed using the R packages
”sp”, “raster”, “rgdal” and “rgeos”57. PA boundaries were calculated using
spatial information from the World Database of Protected Areas58. A
sampling unit was considered as protected if the 1km circular buffer
centered on the unit’s centroid intersected a protected area polygon.
Other explanatory variables corresponding to each wildlife unit were
extracted for the unit’s centroid –for distance variables or gridded
variables with a spatial resolution larger than the unit’s dimensions– or
over a circular buffers centered on the unit’s centroid and passing
through the three vertices –for gridded variables with a spatial
resolution smaller than the unit’s dimensions (Supplementary Table S8).
Estimating protected area effectiveness with matching methods Matching
methods We used the Matchit package59 in R environment, which fits a
logistic generalized linear model where the treatment assignment (land
protection) is the response variable and the matching variables are the
predictors. We chose one-to-one, nearest-neighbor covariate matching with
replacement using a generalized version of the Mahalanobis distance
metric. We used a caliper of 0.2 standard deviations of the propensity
scores as our minimum matching criterion. To assess the quality of the
matches, we checked the resulting covariate balance testing the normalized
differences in covariate means and their distribution. The normalized
difference in means is the difference in the average covariate value
divided by its standard deviation60. We tested for differences in the
distribution using eQQ plots that graph the covariate values in the same
quantile of the treated (protected sites) against those in the control
(non-protected sites), allowing us to observe if characteristics are
distributed similarly across both treatment and control groups (see
Supplementary Figures, S1a and S1b). We were able to match 100% of the
original sample to controls that suited the criteria. Therefore, our
unit-to-unit matching yielded a final dataset of 1220 sampling units, i.e.
610 protected wildlife units that were matched to 610 similar unprotected
units across Finland. These matched pairs of PA and non-PA sampling units
covering, a period of ~30 years, were included in subsequent analyses.
These sampling units were not homogeneously distributed among the four
Game Management clusters with Central Finland being the most represented
cluster and Lapland the least (see Supplementary Table S2). Measuring the
absolute effect of PAs Previous studies have quantified the effect of PAs
in reducing threats in different ways (Absolute PA effect, Relative PA
effect, Pooled relative effect17). However, due to limitations imposed by
the structure of wildlife time series and the natural low densities of
large carnivores (zero inflation of all time-series), an effectiveness
metric based on the ‘absolute PA effect’ was the most relevant approach
for this study. The absolute PA effect is the difference between densities
in a unit located inside a PA and its matched control unit located in a
non-protected area. Therefore a positive value means that sampling units
located inside PA show higher large carnivore densities than its control
unprotected units. We calculated this metric at the national PA network
level, for each pair of ‘protected-non protected’ units. We computed the
median absolute effect of the PA network and its associated 95% confidence
interval for each species of large carnivore, globally and at the regional
level. We performed iterative random sampling (1000 iterations) to control
for differences in the number of repeated observations (number of years)
among pairs of matched units. This was done by selecting only one value of
absolute PA effect for each pair of matched units before estimating the
median PA effect across units. To calculate the absolute PA effect per
region, we pooled the 15 Game Management Areas covering the whole country
in four main clusters: South= ’Satakunta’ + ‘Etelä-Häme’ +
‘Kaakkois-Suomi’ + ‘Uusimaa’ + ‘Varsinais-Suomi’ + ‘Pohjois-Häme’;
Central= ‘Etelä-Savo’ + ‘Rannikko-Pohjanmaa’ + ‘Keski-Suomi’ +
‘Pohjois-Savo’ + ‘Pohjois-Karjala’ + ‘Pohjanmaa’; North= ‘Oulu’ +
‘Kainuu’; Lapland = ‘Lapland’. We compared the annual absolute PA effect
per year with zero using one sample t-tests. Assessing Estimating PA
effectiveness with two-part mixed effects model for semi-continuous data.
We extracted the matched dataset obtained through the matching process
described previously (1220 unpaired wildlife units, see Supplementary
Table S7) to test for the effect of land protection status on large
carnivore densities. We implemented two-part zero-inflated mixed effects
models in the novel GLMMadaptive package that uses adaptive Gaussian
quadrature61. This approach allowed us to account for the data structure
(zero-inflated, right-skewed continuous data, GLMMzi) and assess the
relationship between density indices of the four large carnivore species
and the set of explanatory variables described in Table S8. All models
were implemented in R using the packages GLMMadaptive and parallel. To
account for potential problems of pseudoreplication, unit identity number
was kept consistently as a random effect in all models for each of the 4
large carnivore species. Zero-inflated structures were added on all the
fixed effects included in the models. We reduced the full list of
variables based on co-linearity and biological relevance to produce a set
of 8 variables (e.g. collinearity between distance to closest roads and
distance to closest settlements was too high and we chose to remove
distance to closest roads from all the models). These 8 covariates
included the 6 matching covariates described earlier (percentage of forest
cover, terrain ruggedness, distance to closest settlements, human
population density, latitude and longitude), to which we added year and
protection status of the wildlife unit (protected or not, coded 0/1).
Three interaction terms between fixed effects were also added (between PA
and Latitude, PA and Longitude and PA and years) to assess spatial and
temporal variations in PA effectiveness. A set of 12 models was built for
each species. All 12 models included the 6 matching covariates and
protection status at the unit scale (PA), therefore the difference in
model structure resided in the addition of year, the three interactive
terms and their different possible combinations (structure of the 12
models is described in the Supplementary Table S9). All first-order model
fits were ranked using the Akaike Information Criterion, the best model
having the lowest AIC values from the set of 12 models built for each
species. The fixed effects estimates in mixed models with nonlinear link
functions have an interpretation conditional on the random effects. To
take this into account, we extracted the marginalized coefficients and
their standard errors from the two part mixed models following the
approach described by Hedeker et al.62 and implemented by Rizopoulos61 in
the GLMMadaptive package. Goodness-of-fit of the models was assessed using
residual diagnostics following the procedures described in the DHARMa
package. All statistical analyses were performed using the software R 3.
5. 157. In order to test if PA size could affect our results, we built
four additional models following the same procedure highlighted above.
Data were restricted to protected wildlife units, density of the four
species of carnivores was the response variables and covariates included:
percentage of forest cover, human population density, distance to the
closest settlements, terrain ruggedness, latitude, longitude, year and PA
size. We extracted the marginalized coefficients and checked the residuals
of the different models as highlighted above.