Focal results and statistical analyses
Throughout the manuscript, we focus on the median values of the 31-year time series of OAB, TDM and PLD in both climate scenarios and compare them between lake types. We do so because the objective of our work is to identify general, continental-scale and lake type-specific patterns of phenology and their responses to warming, and not to describe phenological responses to interannual variation in the weather. Similarly, we define the dominant controlling process of OAB in a given lake type at a given geographic location as the process (seasonal change in incident radiation, ice-off, or thermal stratification) which controls OAB in most of the 31 simulated years.
We analyzed the influence of environmental drivers on predicted phenology metrics with generalized additive models (GAMs) (R package mgcv (Wood 2017)). Environmental drivers included geographic location (described by latitude, longitude, and elevation), OD, and the dominant process controlling OAB. Predicted phenology metrics included OAB, TDM, and PLD in both climate scenarios as well as OABdiff, TDMdiff, and PLDdiff. Since OAB and TDM could not be determined for all lake types at all geographical locations and scenarios (see above), GAMs were run on a slightly reduced data set including those 24501 (our of 30512 possible) combinations of lake types x geographic locations for which OAB and TDM were defined in both the reference and warming scenarios (i.e. PLDdiff could be calculated).
We used GAMs to allow for non-linear relationships between environmental drivers and phenology metrics. As we were interested in main effects, we did not model any interactions between environmental drivers. Hence, for each phenology metric, we compared three different models: 1 - a model including only smooth functions of latitude, longitude, and elevation, 2 - a model including only a smooth function of log-transformed OD, and 3 - a model including smooth functions of all four independent variables. For OABdiff, TDMdiff, and PLDdiff, we ran additional models using the dominant process controlling OAB as the sole independent variable. For all models, we report R2 and plot the component smooth functions ± 1 standard error. We do not report p-values as metrics were calculated from the results of deterministic models. Furthermore, due to the large number of observations (grid points x lake types), p values are always highly significant and standard errors are usually too small to show on the plots.