Occupancy modeling
We implemented a multispecies occupancy model of two or more interacting
species (Rota et al., 2016) in program R 3.5.1 (R Core Team, 2019) via
package unmarked (Fiske & Chandler, 2011) to explore how
environmental and anthropogenic variables affect the marginal occupancy
(occupancy without accounting for interactions with other species),
co-occupancy (overlap in marginal occupancy between species), and
conditional occupancy (effects of each species presence on other species
detection and occupancy) of lynx, wildcat, and wolf in the Romanian
Carpathians. Unlike traditional co-occupancy models, multispecies
occupancy models do not requires a priori assumptions of
asymmetric interactions, therefore species were not considered dominant
or subordinate to one another (Rota et al., 2016). Data from the two
seasons were analyzed separately, and sessions were divided into 14-day
sampling occasions, with the winter and autumn seasons having eight and
seven sampling occasions respectively. Camera trap photos were cataloged
by FCC staff and volunteers, and the date, time, location, and species
identification were recorded for each animal detection (Iosif et al.,
2022). Covariates were checked for correlation using Pearson’s
correlation tests and Pearson’s Chi-squared test (for numerical and
factors respectively), those with high correlations r >0.7
were not included in the same models for the same parameter. We first
explored combinations of five detection covariates for species-specific
detection probabilities (Table 1) by comparing models with the same
marginal occupancy parameterization for each species. Detection
covariates were kept the same for all three species as we did not have a
biological reason to vary them between species. We also included the
latent presence/absence of every other species as species-specific
detection covariates (e.g., lynx detection predicted by the
presence/absence of wildcat and wolf). Although multispecies occupancy
models do not assume asymmetric interactions between species, we wanted
to explore the possibility that dominant species could exist in our
system and affect the presence of other species. Therefore, we also
included species-specific detections of lynx as a function of the latent
presence/absence of potentially dominant wolf, and wildcat as a function
of lynx and wolf.
From these models, we determined a best model for each season based on
Akaike Information Criterion (AIC), using R package MuMIn(Bartoń, 2020). We included the top detection covariates in the models
exploring marginal occupancy and co-occupancy. We then ran a series of
models to assess the marginal occupancy of our three species using
environmental and anthropogenic variables (Table 1) that were determineda priori and we hypothesized would affect the marginal occupancy
of each species. The candidate set of marginal occupancy models was
similar for both seasons, models were only removed if variation in
covariates was not great enough to allow estimation (i.e. models
produced NAs or unreasonable estimates and standard errors). We compared
the marginal occupancy models for each season using AIC to identify the
best covariates explaining occupancy of each individual species. Using
the top covariates from the marginal occupancy analysis, we ran a series
of additional candidate models that reflected a priori hypotheses
regarding pairwise co-occupancy between lynx and wildcat, lynx and wolf,
and wildcat and wolf, and compared the models using AIC and biological
relevance (Table 2). Due to data limitations (small sample size), we did
not implement a three-species co-occupancy parameterization.