Heat stress reduces growth rate of red deer calf: Climate warming implications

are flexible on the shape of the response. GAM models were implemented using the R software package mgcv with multiple smoothing parameter estimation by restricted marginal likelihood [42]. Both analyses used two types of explanatory terms, animal and physical environment vari- ables. The animal terms were the calf age (days) and mother age (years), calf sex and their per- tinent interactions, and were fitted as fixed effects in GAM and non-linear mixed effects models. In addition, calf and mother identity, and calf cohort were fitted as random effects. The second type of variables fitted in the above models were those that represent heat environ- ment. More concretely, the two indices of heat stress defined above (THI, THIWS), and the aforementioned four physical environment variables (T, Hr, W, SR). Although heat stress indi- ces are widely used in the literature [43], as any other composite index, they have the problem that when used in statistical analyses the particular relationships and error structures between the response variable and the variables that constitute the index are lost [44]. For the analysis of calf weight we undertook three non-linear mixed effects models fitting (i) animal variables and THI; (ii) animal variables and THIWS; and (iii) animal and physical environment vari- ables. Similarly, for mother weight but using GAM models. Although GAM models are additive by definition, we fitted interactions between explana- tory variables using the ti function within GAM model formulae. This produces a tensor prod-

uct interaction which is appropriate when the main effects and any lower interactions are also present [45, 46]. To deal with independent random effects, these were treated as smooths (using “re” basis in model formulae), the terms produce a parametric interaction of the predic- tors and penalize the corresponding coefficients with a multiple of the identity matrix, corre- sponding to an assumption of normality [45]. We assumed that the response could be non- linear but relatively simple, and because penalized regression smoothers gain computational efficiency using a basis of modest size [45] we used small k values (k = 3, 4). The non-linear mixed models evaluated an asymptotic regression function and its gradient of the type, f ðxÞ ¼ Asym þ ðR0 AsymÞ � e elrc�t Eq 3 where t is calf age in days; Asym represents the horizontal asymptote; R0 is the response at t = 0; and lrc is the natural logarithm of the rate constant. The starting values of the parameters of this function were initially estimated using calves body weights against age through lactation using the self-starting model implemented in the nonlinear least-squares regression package (nlm) of R software [47]. The starting parameters of the other explanatory variables and their pertinent interactions were guessed and fitted in nlme, the new estimated parameters by the model were fitted as starting parameters and the model was re-run to corroborate consistency in the parametrisation of the model, which was achieved in all cases. Model selection was per-

formed using Akaike (AIC) weights aided by the normalised probability of the Kullback–Lei- bler discrepancy ratio; in which model A is to be preferred over competing model B [48]. The coefficients of the linear mixed model were calculated using REML, as the estimates are more accurate than using maximum likelihood [49]. To improve the readability of the models out- put we plotted the response of interest by fixing the other explanatory variables to their mean values. Graphics were constructed in R [47] using the ggplot2 package based on The Grammar of Graphics [50]. Results Physical environment The mean date of birth was 23th of May (sd = 10.2 d, n = 583) and through the next 143 days till weaning and across the 19 years of the study, mothers and offspring experienced the PLOS ONE Heat stress and calf growth in red deer PLOS ONE | https://doi.org/10.1371/journal.pone.0233809 June 1, 2020 5 / 26