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 ).