Statistical analysis
We used t-tests to compare the two nest box temperature measurements
Ta and Tbr of our control and heated
nests during the six hours of daily treatment (see above). For each nest
we compared growth rates (K ) for nestling mass using nonlinear
mixed models (package nlme Pinheiro et al. 2017) following the
methodology of Sofaer et al. 2013. This approach uses age and
treatment as fixed effects, and nest and nestling identity as random
effects to account for the lack of independence among siblings of the
same brood and for repeated measures of the same individual over time.
We initially tested the effects of simulated heat waves by fitting
linear mixed models using the packages lme4 and lmerTest (Bates et
al. 2014; Kuznetsova et al. 2017) with telomere length as the
dependent variable, body mass, growth rate, as covariates, and brood
size, hatching order, treatment, sex and age plus their interactions as
categorical fixed factors. Due to the strong correlation between
offspring mass and age during post-natal growth we created a variable
called mass index that encapsulated the standardized differences in mass
according to the age class of each subject, with a high mass index
meaning that individuals are relatively heavy for their age and vice
versa (see more details in S3.1). We included random slopes for
individuals to account for repeated measures of subjects belonging to
the same brood at different times. From this initial model we dropped
the non-significant variables and fitted our best model that included
only treatment, age, mass index and the interaction between treatment
and age, as well as individual identity as random factor (see Section
S3.2). To ensure that we were not over-fitting the data by including the
interaction between treatment and age, we also performed a model
comparison test using Akaike Information Criterion (AIC) (Akaike 1987)
and Bayesian Information Criterion (BIC) (Schwarz 1978). This test
indicated that our model was preferred by both approaches (AICweights =
99%; BICweights = 80%) to a model with an additive effect of treatment
and age and to another with no interaction at all (see Tables S3.2.3).
We conducted post-hoc pairwise comparisons between groups using the
emmeans package (Lenth et al. 2019), which applies Tukey
adjustment and Kenward-Roger degree of freedom adjustments. We also
built a General Linear Model with time spent brooding as a dependent
variable and treatment as a fixed factor to test the effects of heat
waves on time spent brooding by the parents. Since this latter model
indicated an effect of treatment on brooding effort, we attempted a test
of the effects that variation in brooding rate has on telomere length
and attrition by including in our initial model brooding and its
interaction with age and treatment as a covariate. This part of the
analysis however needs to be interpreted with extreme caution (see
sections S6). Indeed, AIC suggested the best model to be the one
including the three-way interaction while BIC gave opposite indications;
also results were highly sensitive to influential observations (S6.3.3).
All models assumptions and performances (homogeneity of variance,
normality of residuals, normality of random effects, collinearity, and
influential observations) were checked (see Supplemental Materials). All
analyses were performed with R version 4.1.2 for Mac (R Core Team
2021).