Trends in abundance and phenology: We estimated relative abundance using ‘list length’ analysis (Breed et al., 2013; Szabo et al., 2010; Franklin, 1999). List length analysis uses the total number of species observed in an outing, i.e. the list’s length, as a proxy for search effort. List length may be more reliable than observer hours, observer miles, or number of observers because the probability of seeing a species does not scale directly with these metrics for observation (Szabo et al., 2010; Franklin, 1991). Relative abundance of each species was estimated using a binomial generalized linear mixed model (GLMM) with fixed effects of year and list length, and a presence/absence response (see Appendix S2 ). We also included year as a random effect to account for interannual variation in relative abundance of species. The yearly slope coefficients of these models represent changes in relative abundance, measured as the change in species presence on lists over time. Analyses were implemented in R using the lme4 package (R Core Team, 2018; Bates et al., 2014). To supplement inference from list length analysis, we analyzed trends in the number of butterflies of each species seen per trip. These counts are a very coarse metric of abundance because they were collected idiosyncratically (not with a uniform survey protocol). However, congruent results from analyzing these counts and list length analysis would suggest that both are reasonable metrics of abundance. Count models used Poisson GLMMs with fixed effects of year and random effects of year. Models were fit using the glmer() function from the package ‘lme4’ (Bates et al., 2014).
We used quantile regression (Cade & Noon, 2003) to estimate trends in the onset, end, and total length of observations of each species across years. Changes in the onset of flight were estimated as the slope of the 0.1 quantile of observation dates for each species as a function of year. Therefore, negative values indicate advances in phenology and positive values indicate delays in phenology. Changes in the end of flight were estimated as the slope of the 0.9 quantile of observation dates of each species as a function of year. Therefore, negative values indicate advances and positive values indicate delays. The change in the total flight period was estimated as the difference between the 0.9 and 0.1 quantile slopes, where negative values indicate shortening flight periods and positive values indicate lengthening flight periods. Quantile regression was carried out using the rq() function from the package ‘quantreg’ (Koenker, 2019). To make an explicit comparison between the changes in phenology of Massachusetts and British butterflies from Macgregor et al. (2019), we also use linear regression to measure the annual change in mean observation date as a function of time with year included as both a fixed and random effect. Negative values again indicate advances in the mean observation date and positive values indicate delays. We used the glmer() function from the package ‘lme4’ (Bates et al., 2014).
As a simple metric of the association between changes in phenology and changes in abundance, we correlated metrics of phenological change (slope for each species with respect to time of the 0.1 quantile, 0.9 quantile, flight period, and mean date) with changes in abundance (slope for each species from list length analysis). Because all slopes were measured with error, we verified the relationship with parametric bootstrapping, i.e., resampling estimates for each species from a random normal distribution using the true mean and standard deviation for each species (10,000 bootstrap samples). We evaluated taxonomic effects as predictors of all response variables. Trends in abundance and various phenology metrics were not associated with taxonomic groups (at the family or genus level), and these effects are not discussed further (see Appendix S3 ).