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.