2.5 Statistical analysis
Changes in pollinator species richness were evaluated in response to
green area fragmentation (i.e. , the variable Edge Density) and
flower richness (i.e. , the number of flowering species per site).
To do this, a Generalised linear mixed model (GLMM) regression (glmmTMB
R package; Magnusson et al., 2017) with Poisson distribution was used,
with island included as a random effect. The flower richness was
included as a predictor along with the edge density, since it could
represent an important local driver of pollinator richness (Blüthgen &
Klein 2011).
The same variables (edge density, and flower richness) were used along
with the network size as predictors of change of the Connectance network
index. Network size, calculated as the product between the number of
insects and plants included in the networks for each site, was also
included as a predictor in the model to account for its effect on
Connectance variation (as in Biella et al., 2020). In this case, a GLMM
with beta distribution and island included as a random effect was
used.
Changes in individual pollinator degree were evaluated in response to
green area fragmentation and flower richness. The effects of these
covariates were evaluated in interaction with the pollinator species
identity to highlight differences among the considered pollinator
species. A GLMM with Poisson distribution was used, with sites nested in
the island as a random effect.
Variation in the pollination efficiency was evaluated using the
pollinator richness and the Connectance as covariates. Moreover, the
plant degree (mean number of pollinator species interacting with each
plant species considered in pollination efficiency analysis) was
calculated from DNA metabarcoding data to estimate the mean plant
generalism and included as model covariate. The role of these covariates
was evaluated in interaction with the plant species identity, to
highlight differences among the investigated plants. A GLMM with
negative binomial distribution was used to account for overdispersion.
Also in this case, the site nested in island was included as a random
effect.
All the analyses were performed with R (version 3.6.1; R CoreTeam 2019).
Predictor significance was evaluated through a log likelihood ratio test
(P < 0.05). The Vif function of the car package (i.e. ,
Variance Inflation Factor with an exclusion threshold of 3 (Zuur, Ieno,
& Smith, 2007)) was used to exclude collinearity among variables. In
all cases, the final models were obtained by removing the variables that
did not improve the model fit through backward stepwise regression based
on second-order Akaike Information Criterion (AIC) (Zuur et al., 2009)
calculated with the package MuMIn.