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