Statistical analysis
We used linear mixed models (LMM) to examine the effects of climate,
urbanization, and life history traits on the emergence, termination, and
duration of adult insect activity across North America. Estimates of
emergence, termination, and duration were the response variables, and we
included human population density, mean annual temperature, annual
precipitation, temperature seasonality, precipitation seasonality as
predictor variables along with two, two-way interactions based on
previous biological knowledge. The two-way interaction terms tested
were: 1) whether the effects of annual temperature changed along a
precipitation gradient and 2) whether the impacts of human population
density changed along a temperature gradient. Areas with higher human
population densities may have more incidental data records, biasing
phenological metrics and we controlled for this by including the number
of observations as a fixed effect in our initial model. We scaled
variables to have a mean of zero and standard deviation (s.d.) of one to
ensure comparable model effect sizes across variables. Cell identity and
insect species were included as random terms for intercepts, and insect
species were included as a random term for the slope of each predictor
variable. Next, we used backward model selection using the step function
from the R package lmerTest(Kuznetsova et al.,
2017). If any variance inflation factor (VIF) was greater than five for
a variable in the best model, we removed the correlated variable with
lower coefficient and reran our backward selection process.
After selecting the best model with only climate, human population
density, and number of observation variables, we added insect traits to
our models. We added two-way interaction terms between each trait
variable with the climate and human population density variables to
examine if adult insect phenology along climate and population density
gradients changed for each trait. Backward model selection was performed
as described above to reduce model variables and select a best model.
Again, if any VIF was greater than five for the variables in the best
model, we removed the correlated variable with the lowest coefficient
and reran our backward selection process.
LMM can lead to false conclusions and inflated type I error rates if
phylogenetic relationships are ignored
(Li and Ives, 2017).
We therefore generated a subtree from the Open Tree of Life for the 101
insects in our analysis using the R package rotl(Michonneau et al.,
2016). Branch lengths were generated by searching the TimeTree of Life
database (Kumar et
al., 2017) to get the estimated divergence time of each internal node.
The branch lengths were then scaled from these ages using the
ph_bladj() function from the R package phylocomr(Ooms and Chamberlain,
2019). We used this synthesis phylogeny to fit phylogenetic linear
mixed models (PLMM) using a Bayesian framework with the default
uninformative INLA priors
(Rue et al., 2009).
We used the R package phyr(Li et al., 2020) to
fit our top LMM as PLMM. Our results differed slightly between the PLMM
and LMM, and therefore, we present the results based on PLMM in the main
text. Results of linear mixed models are included in the supporting
information (Supporting Information Tables 1-3). We checked residuals of
the top models to ensure models did not have any obvious deviation from
model assumptions, including spatial autocorrelations (See Supporting
Information for more details and Figure S1 for residuals of the top
models). We measured the goodness of fit of our PLMM using the packagerr2 (Ives and
Li, 2018) to generate partial R2 values (see
supporting information for additional details).