Statistical analysis
All analyses were run in the R (R Core Team, 2020) with packagesMuMIn (Bartoń, 2020) and lme4 (Bates et al., 2018). We analysed each of the response variables separately with generalized linear mixed-effects models. We tested the effects of urbanization, local canopy cover and their interaction on leaf damage with Gaussian error distribution and identity link (the results were the same with a beta-distribution and log-link), and on the incidence of leaf-miners and gall-inducers with binomial error distribution and logit-link in separate models. The data was not overdispersed, visual inspection of raw-data did not call for zero-inflated models, the distribution of residuals met model assumptions.
In each model, Urbanization (%), Local canopy cover (%), Urbanization × Local canopy cover, Year (as a factor), Spring temperature (°C) and Spring precipitation (mm) were included as fixed effects; and Partner ID as a random factor to account for the fact that some partners surveyed multiple trees and/or several years.
We analysed the data in the framework of information theory (Burnham and Anderson, 2002). We first built three models, one for each response variable separately (leaf damage, gall-inducer and leaf-miner incidences). We scaled and centred all continuous predictor variables prior to modelling to make their coefficients comparable, and verified that uncontrolled correlations among explanatory variables were unlikely to bias model coefficient parameter estimates (all variance inflation factors [VIF] lower than 2) (Zuur et al., 2009). We then applied a procedure of parsimonious model selection based on the Akaike’s Information Criterion corrected for small sample sizes (AICc) and considered every model in a range of 2 units of AICc to the best model as equally likely (Arnold, 2010). We calculated the AIC weight (wi ) – i.e., the probability that a given model is the best model within the set of candidate models – and also the relative variable importance (RVI) as the sum ofwi of every model including this variable. When multiple models were competing with the best model (i.e., when several models with ΔAICc < 2), we implemented a multi-model inference approach, constructing a consensus model that comprised the selected variables from the set of best models. We subsequently averaged their effect sizes over all models in the set of best models, utilizingwi as the weighting parameter (i.e., model averaging). A certain predictor was deemed to have a statistically significant effect on the response variable if its 95% confidence interval did not bracket zero.