What this study adds
This study shows that the interannual variability of heat-related mortality might be partly driven by large-scale climatic teleconnections, specifically the Atlantic Multidecadal Oscillation in the case of the province of Quebec. To do so, we provide a rigorous methodology that uses recent relevant statistical development, in particular from functional data analysis that allows the consideration of several climate indices and lags. As far as the authors are aware, this study is the first of its kind, and we think it is valuable both for its informative results and its methodology.

Introduction
Heat-related mortality is an important public health burden in summer^{1} and is expected to increase with climate change.^{2} ^{,} ^{3} To inform adaptation policies, studies have focused on modeling the evolution of mortality attributable to heat^{4} or relating the vulnerabilities across the world to socio-economic drivers.^{5} ^{,} ^{6} However, in addition to the long-term trend and to spatial variations, there is important interannual variability in summer heat and thus in heat-related mortality.^{7} It is therefore of interest to model and predict these inter-annual changes to inform short-term action plans.

Temperature is widely driven by climatic teleconnections which contribute to define large-scale climatic variations. For instance, long-term cycles in sea surface temperature in the Atlantic Ocean are strongly linked to the North American and European summer climates. Climatic teleconnections and large-scale patterns are summarized by different climate indices that represent diverse variations across the globe, such as the Southern Oscillation Index (SOI) measuring the so-called El-Niño phenomenon or the Atlantic Multidecadal Oscillation (AMO) to represent variations of sea surface temperatures in the Atlantic.^{8} These indices can then be used to predict local climate characteristics like the frequency and intensity of heatwaves.^{9–11} In addition to heat waves, these indices also predict variations of other weather variables such as precipitations,^{12–14} thus potentially affecting the vulnerability of populations to extreme meteorological events.

Given the strong links between climate indices and local weather, they are expected to be instrumental for the prediction of heat-related mortality. As accurate forecasts of climate indices several months in advance are now a prospect,^{15} they can be used to provide heat-related health forecasts months in advance. This strategy is often considered for infectious diseases that are particularly climate-sensitive,^{16} ^{,} ^{17} but less so for heat-related mortality.

One limiting factor of integrating large-scale climate information is the difficult statistical problem of predicting heat-related mortality by using potentially numerous interrelated and lagged climate indices. In the only contribution attempting it, Majeed et al.^{18} focused on two indices for the United-Stated: the ENSO phenomenon and the AMO. They used mortality causes defined a priori for heat-related mortality and performed a spectral analysis of time series. However, such an approach does not necessarily represent accurately the mortality actually caused by heat, nor allow disentangling the different climate indices that could impact summer heat.

We propose here a two-stage statistical analysis that consists in first estimating heat-related mortality for each year, and second to model annual heat-related mortality using values from climate indices. Similar two-stage methodologies have been employed to study adaptation to heat and cold^{4} as well as the differential vulnerability patterns over the world.^{5} The first stage uses state-of-the-art epidemiological models to estimate the part of mortality that can be attributed to heat. The second stage considers functional regression to efficiently integrate the whole information provided by several months of climate indices and predict the magnitude of upcoming heat-related mortality.^{19} Functional regression takes continuous functions as inputs and is thus particularly adapted to time-indexed processes. In particular, we use the functional linear array model recently proposed by Brockhaus et al.^{20} to manage the high-dimensionality of the model. Its boosting-based algorithm is efficient both for prediction and for variable selection.^{21} This allows to select the most predictive indices as well as understand how their variations can help predict heat-related mortality.

The objective of the present article is to relate mortality attributable to heat with large-scale climate indices using the two-stage methodology sketched above. As climate indices, we consider both the AMO and the SOI measuring the ENSO phenomenon, as well as five additional potentially predictive climate indices: the Arctic Oscillation (AO), the North Atlantic Oscillation (NAO), the Oceanic Niño Index (ONI), the Pacific Decadal Oscillation (PDO), and the Pacific North American Index (PNA). The obtained model can then be used to predict the magnitude of upcoming heat-related mortality and inform heat early warning systems several months in advance.

Methods
Data
We consider data from two metropolitan areas (MA) in the province of Québec, Canada: Montréal and Québec. These communities correspond to the extended urban areas that include neighbor municipalities such as Laval for Montréal and Lévis for Québec as shown in Figure 1 . For both MAs, we consider daily mean temperature and all-cause mortality counts as well as monthly climate indices over the period 1981–2018.

Figure 1.: Map of the considered metropolitan areas: Montréal (green) and Québec (blue).

Daily time series are herein restricted to the five warmest months of the year which are May to September. All-cause mortality counts are provided for each MA by the National Institute of Public Health of Québec (Institut national de santé publique du Québec ). Mean temperature is constructed as the average between minimum and maximum temperature, both extracted from the Daymet database (daymet.ornl.gov). Daymet is a gridded weather database with a resolution of 1 km over North America.^{22} It is bias-corrected and validated using ground observations and its high resolution allows for an accurate representation of exposure to temperature.^{23} Daily time series of mean temperature are obtained by averaging over all grid points that cover each MA.

Monthly values of the seven climate indices given above are obtained from the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (https://psl.noaa.gov/data/climateindices/list/ ). These seven indices are suspected to influence the climate of Southern Québec to different degrees.^{24} Each index is detailed in Appendix A; https://links.lww.com/EE/A184 .

Statistical analysis
To study the ability of climate indices to predict heat-related mortality, we consider a two-stage approach applied to each MA. In the first stage, we estimate annual heat-attributable fractions (AF) of mortality through a distributed lag nonlinear model. In the second, AF estimates are considered as the response in a functional regression model with climate index curves as predictors.

Estimation of heat attributable fractions
To estimate heat AF in the first stage, we apply a time-varying distributed lag nonlinear model. It consists in a quasi-Poisson regression with a cross-basis to represent the nonlinear and delayed effect of temperature on all-cause mortality.^{25} The cross-basis is parameterized following previous studies on heat-related mortality,^{26} i.e., through a quadratic B-spline with knots at percentiles 50% and 90% on the temperature dimension, and a cubic natural spline with two interior knots uniformly distributed on the log scale on the lag dimension, up to a maximum lag of 10 days. To allow the exposure-response function to evolve along the study period, and thus account for potential long-term adaptation, we add an interaction term between the cross-basis and the date as described by Gasparrini et al.^{26} The model also includes a day-of-week factor to control for weekend effects, a natural spline component of year with one degree of freedom per decade to control for a long-term trend and a natural spline of day-of-season with four degrees of freedom to control for seasonality.

The fitted model is then used to compute daily AFs for heat using classical formulas.^{27} The daily AFs are then aggregated by years including only heat days, defined as days for which mean temperature exceeds a specific threshold. Several such thresholds are considered to test whether overall heat or only the most extreme events can be predicted by climate indices (see examples in Appendix B; https://links.lww.com/EE/A184 ).^{28} First, we consider the minimum mortality temperature (MMT) for each year, i.e., we retrieve the overall cumulative exposure-response function^{29} at the first of July of each year and extract the temperature at which the risk is the lowest, with the constraint that it is between its 10^{th} and 90th percentile.^{30} We also consider percentiles 95%, 97.5%, and 99% of summer temperature, computed on the whole period. We focus on AF instead of attributable number (AN) as it is less prone to confounding by population.

Uncertainty of reported MMTs and AFs is assessed by 95% empirical confidence intervals (eCI) that simulate new first-stage coefficients from a multivariate normal distribution and computes the MMT and AFs for each set of simulated coefficients.^{27}

Prediction by climate indices
For each MA and each heat threshold, estimated AFs are predicted by climate index curves through a functional regression model in which the response is scalar and the predictors are functional.^{31} This model estimates the effect of each considered month of the climate indices on AF through the following model:

$$A{F}_{i}={\displaystyle \sum}_{j=1}^{7}\underset{0}{\overset{16}{{\displaystyle \int}}}{\beta}_{j}(l){x}_{ij}\left(l\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}dl\text{}+s\left(i\right)+{\epsilon}_{i},$$

where $A{F}_{i}$ is the heat AF of year i , ${x}_{ij}\left(l\right)$ is the curve of the ${j}^{th}$ climate index ($j=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\dots \text{\hspace{0.17em}}\text{\hspace{0.17em}},7$ ) for year i with the associated functional coefficient ${\beta}_{j}(l)$ . This coefficient indicates the change in AF induced by ${x}_{.j}\left(l\right)$ at each lag l . As the impact of climate indices on the local climate can be delayed by up to two years,^{12} in this study ${x}_{.j}\left(l\right)$ spans a period of 16 months from January of year $i-1$ to May of year i . Finally, $s\left(i\right)$ is a smooth component of the year expanded by P-splines, i.e., a dense cubic B-spline with 20 interior knots that is penalized to control its smoothness.^{32} The last component ${\epsilon}_{i}$ represents the residuals.

To estimate model (1), we consider the functional linear array model representation with associated boosting optimization developed by Brockhaus al.^{20} In the linear array representation, the functional variables ${\beta}_{j}(l)$ and ${x}_{ij}\left(l\right)$ are expanded by spline tensor products to map discrete data to the continuous domain of $l.$ Using this expansion, the model is then estimated through gradient boosting, which iteratively constructs the model by adding a very simple component of a single variable at each step to reduce the residual sum of squares.^{21} This algorithm served the double purpose of optimization as the residual sum of squares is iteratively reduced while the model gains in complexity, and variable selection as only the most strongly predictive explanatory variables are added in the model. The number of iterations is chosen by 10-fold cross-validation. More details on the estimation procedure and chosen number of iterations are given in Appendix C; https://links.lww.com/EE/A184 .

To evaluate the predictive power of each model, we compute the R ^{2} through 10-fold cross-validation, a procedure that does not favor more complex models, in contrast with its training sample counterpart.^{33} As an additional sensitivity analysis, we apply a model that adds the average summer temperature as an additional term in model (1) for each outcome and each location. Results of this sensitivity analysis are reported in Appendix D (https://links.lww.com/EE/A184 ) and show consistency with the main results reported below.

The analysis is performed using R 4.1.0^{34} with the addition of packages dlnm^{35} for first-stage heat-related mortality estimation and FDboost^{36} for the second-stage functional regression.

Results
Summary statistics of deaths and temperature are reported in Table 1 for each MA. The study covers more than 400,000 deaths, with five times more deaths reported in Montréal than Québec. Montréal also shows a higher average temperature than Québec, being on average 3°C warmer. MMTs are also much warmer in Montréal, ranging from 16.3°C (95%eCI: 13.8–18.8) in 1981 to 19.1°C (95%eCI: 11.2–23.2) in 2018, compared to MMTs of 8.0°C (95%eCI: 8.0–20.0) in 1981 that increased to 14.3°C (95%eCI: 10.3–16.3) in 2018 for Québec.

Table 1. -
Total number of deaths and mean temperature summary by area.

Area
Deaths
Daily mean temperature (mean) (range)
Annual minimum mortality temperature (MMT) range (95%eCI)
Montréal
354 764
17.6 (1.3–29.2)
16.3 (13.8–18.8) 19.1 (11.2–23.2)
Québec
73 199
14.4 (−1.4–26.6)
8.0 (8.0–20.0) 14.3 (10.3–16.3)

Figure 2 shows AF estimated for both MAs and different percentiles defining heat-related mortality, with detailed numbers and eCI in Tables S1 and S2; https://links.lww.com/EE/A184 . Related results such as dose-response curves and summary statistics are reported in Appendix C; https://links.lww.com/EE/A184 . On average 2.42% (95%eCI: 2.02–3.48) of mortality is attributed to heat above the MMT in Montréal and 2.84% (95%eCI: 1.32–5.88) in Québec, and the AF decreases with the percentile considered with an average of 0.40% (95%eCI: 0.36–0.45) in Montréal and 0.24% (95%eCI: 0.15–0.34) in Québec for heat defined with the 99th percentile. In both MAs, AF widely varies from year to year with a downward trend visible for MMT-related AFs. This downward trend is accompanied by an increase of MMT in recent years (see Figure S2; https://links.lww.com/EE/A184 in Appendix B) and suggests a possible diminution of the overall vulnerability to heat, consistently with previous studies.^{4} ^{,} ^{26} Figure 2 also highlights years with peaks for both mortality attributed to heat and AMO values such as 1983, 1988, 1999, 2005, and 2010.

Figure 2.: Estimated attributable fractions (AF) for each year, area and percentile (top), and the average of the same year March to May AMO (bottom). Grey areas indicate periods when both AF and AMO peak.

For both MAs and all heat definitions, results from the second-stage functional regression model are similar. Details on intermediate results are shown in Appendix D; https://links.lww.com/EE/A184 . In all cases, AMO was the single climate index selected, except for Québec and the MMT-related AF for which the model selected only the time component due to the important trend. Figure 3 shows the estimated functional coefficients $\beta {\widehat{}}_{j}\left(.\right)$ from Equation 1. In all cases, the coefficient is positive for lags 0 to 5, i.e., in the same year, winter and spring, and are close to zero afterward, with negative values during winter of the previous year. Therefore, a low AMO value in winter and spring is often followed by an important AMO value in winter and spring of the following year and then higher than usual heat-related mortality in summer. This is visually confirmed in Figure 2 which shows a similar pattern between AFs and the average of the same year January to May AMO. Figure 3 also indicates a higher amplitude of scaled $\beta {\widehat{}}_{j}\left(.\right)$ for when AFs are defined with the MMT in Montréal. This suggests a link between the AMO and the overall heat-related mortality, not only with the most extreme heat waves.

Figure 3.: Scaled functional coefficients

$\beta {\widehat{}}_{j}\left(.\right)$

with their 95% confidence intervals in Montréal and Québec for all heat definitions. The curves span 16 months lag from the beginning of summer (lag 0) on the left to January of the previous year on the right (lag 16).

Table 2 indicates the cross-validated R ^{2} statistic for each model. The prediction of AF is better in Québec than in Montréal with R ^{2} values decreasing from 21.8% for MMT-related AF to 2.0 % for the 99th percentile AF for Montréal and decreasing from 51.9% to 14.5% for Québec. These numbers are higher in Québec because of the important decreasing trends of AF that are mainly predicted by the time component while the AMO predicts the year-to-year variations.

Table 2. -
Cross-validated

R ^{2} (standard deviation) expressed in percentage for each second-stage model.

Percentile defining heat
Montréal
Québec
MMT
20.1 (1.9)
51.8 (0.8)
95th
1.8 (0.4)
21.5 (2.7)
97.5th
2.0 (0.5)
26.0 (4.7)
99th
2.0 (0.4)
13.2 (1.7)

Discussion
In the present work, we predicted summer heat-related mortality in two MAs of the province of Québec using climate indices. Results indicated that a non-negligible part of heat-related mortality can be well predicted by the AMO index, especially its values in the late winter and spring of the same year. Predictions are also better for the whole heat-related AF than for the most extreme cases. These results suggest that monitoring the AMO is relevant for public health authorities in the province of Québec to anticipate and adjust emergency preparedness for heat-related mortality on an annual basis.

The AMO index measures north-Atlantic sea surface temperatures, which shows an important impact on the North American weather.^{37} ^{,} ^{38} For instance, positive (warmer) AMO phases are associated with a decrease in rainfalls as well as an increase in droughts in the United States.^{39} In the province of Québec, AMO has been associated with generally hotter summers,^{24} as well as longer and warmer heat spells.^{11} Therefore, it is natural that AMO is partly predictive of heat-related mortality.

AMO is predictive of heat-related mortality in both Montréal and Québec, despite strong geographical and adaptation trends disparities between the two cities. Indeed, the MMT differs by approximately 6°C which confirms previous works indicating higher risks for similar temperature in Québec compared to Montréal.^{40} A variety of factors can explain this difference, including generally lower temperatures in Québec, a less pronounced urban heat island effect being a more spread and less dense city, as well as closer proximity to the Atlantic Ocean.^{41} First-stage results also indicate a decreasing heat-related risk during the study period, which is in line with previous results in Québec^{42} and across the world.^{26}

Results suggest that AMO is not only predictive of extreme heat waves-related mortality but of the whole heat-related mortality. Indeed, in Québec, the whole temperature distribution seems correlated to climate teleconnections.^{24} Besides, the number of extreme events is still low over less than forty years of data, meaningless data to properly train predictive models. Therefore, further research is needed to be able to predict the deadliest heatwave events.

The main strength of this work is the statistical methodology that uses state-of-the-art approaches in both stages. First, heat AFs are estimated by time-varying distributed lag nonlinear models , a methodology recently used to highlight long-term adaption phenomena across the world.^{4} ^{,} ^{26} Second, the associations between climate indices and heat AF are estimated through functional regression, an increasingly popular methodology in many fields. It has been considered a few times in epidemiology,^{43} including environmental epidemiology,^{44} ^{,} ^{45} as well as in climatic studies to forecast temperature.^{46} Its efficient use of intrinsically continuous variables allows the inclusion of monthly records of climate indices, without the need to remove the highest frequencies and select specific months beforehand as was done before.^{18} ^{,} ^{24}

The present work shows several shortcomings, due to the complexity of climate and health relationships. The number of years considered is still low to fully apprehend the impact of teleconnections on heat-related health, especially for the most extreme events. Indeed, many teleconnections represented by the considered climate indices present low-frequency patterns that cannot be fully captured by only 38 years of data. For instance, one of the main periodicities of AMO is several decades long and this mode is not captured by our dataset, although the smaller scale oscillations of AMO predict well heat-related mortality. In addition, the second stage model is fitted with only 38 individual records, preventing the fitting of very complex models with several indices and potential interactions. Therefore, studies relating large-scale climate patterns and heat may gain strength as data accumulates in the upcoming years.

The second shortcoming, as suggested by the prediction criteria computed here, is that many factors that drive the variation in heat-related mortality could not be accounted for in this study. These include the prevalence of air conditioning^{6} or the introduction of early heat-health warnings in the province of Québec.^{47} ^{,} ^{48} In addition, the already complex time-varying DLNM fitted in this work does not capture short-term variations in the heat-related risk such as seasonal adaptation, that could alter the true AF.^{49}

More research is needed to better understand how climatic teleconnections can predict heat-related mortality and how to efficiently integrate such information in heat-wave warning systems. This includes studies with much longer time series of mortality and temperature. As the present work focuses on the province of Québec, it is important to replicate it in locations with different climates in which other climatic indices will be better predictors of heat-related mortality. In addition, the present work focuses on all-cause mortality and specific mortality causes could be best predicted by teleconnections. Finally, the present work could also be extended to cold-related mortality, as winter extreme cold and precipitations are also widely linked to the different climatic patterns.^{50}

The main conclusion of the present work is the possibility to predict summer heat-related mortality using climate indices representing large-scale teleconnections. For the specific case of the Québec province, it was found that the best index for this task is the AMO. This information can help inform public health seasonal heat action plans.

ACKNOWLEDGMENTS
Conflicts of interest statement
The authors declare that they have no conflicts of interest with regard to the content of this report.

References
1. Weinberger KR, Harris D, Spangler KR, Zanobetti A, Wellenius GA. Estimating the number of excess deaths attributable to heat in 297 United States counties. Environ Epidemiol. 2020;4:e096.

2. Vicedo-Cabrera AM, Guo Y, Sera F, et al. Temperature-related mortality impacts under and beyond Paris Agreement climate change scenarios. Clim Change. 2018;150:391–402.

3. Wang Y, Wang A, Zhai J, et al. Tens of thousands additional deaths annually in cities of China between 1.5°C and 2.0°C warming. Nat Commun. 2019;10:1–7.

4. Vicedo-Cabrera AM, Sera F, Guo Y, et al. A multi-country analysis on potential adaptive mechanisms to cold and heat in a changing climate. Environ Int. 2018;111:239–246.

5. Sera F, Armstrong B, Tobias A, et al. How urban characteristics affect vulnerability to heat and cold: a multi-country analysis. Int J Epidemiol. 2019;48:1101–1112.

6. Sera F, Hashizume M, Honda Y, et al. Air conditioning and heat-related mortality: A Multi-country Longitudinal Study. Epidemiology. 2020;31:779–787.

7. Scortichini M, de’Donato F, De Sario M, et al. The inter-annual variability of heat-related mortality in nine European cities (1990-2010). Environ Health. 2018;17:66.

8. Yuan X, Kaplan MR, Cane MA. The interconnected global climate system-a review of tropical-polar teleconnections. J Clim. 2018;31:5765–5792.

9. Loikith PC, Broccoli AJ. The influence of recurrent modes of climate variability on the occurrence of winter and summer extreme temperatures over North America. J Clim. 2014;27:1600–1618.

10. Ouarda TBMJ, Charron C, Kumar KN, Phanikumar DV, Molini A, Basha G. Nonstationary warm spell frequency analysis integrating climate variability and change with application to the Middle East. Clim Dyn. 2019;53:5329–5347.

11. Ouarda TBMJ, Charron C. Nonstationary temperature-duration-frequency curves. Sci Rep. 2018;8:15493.

12. Ouachani R, Bargaoui Z, Ouarda T. Power of teleconnection patterns on precipitation and streamflow variability of upper Medjerda Basin. Int J Climatol. 2013;33:58–76.

13. Ouarda TBMJ, Yousef LA, Charron C. Non-stationary intensity-duration-frequency curves integrating information concerning teleconnections and climate change. Int J Climatol. 2019;39:2306–2323.

14. Rust W, Cuthbert M, Bloomfield J, Corstanje R, Howden N, Holman I. Exploring the role of hydrological pathways in modulating multi-annual climate teleconnection periodicities from UK rainfall to streamflow. Hydrol Earth Syst Sci. 2021;25:2223–2237.

15. Soares MB, Daly M, Dessai S. Assessing the value of seasonal climate forecasts for decision-making. Wiley Interdiscip Rev-Clim Change. 2018;9:e523.

16. Lowe R, Stewart-Ibarra AM, Petrova D, et al. Climate services for health: predicting the evolution of the 2016 dengue season in Machala, Ecuador. Lancet Planet Health. 2017;1:e142–e151.

17. Zaitchik BF, Hayden MH, Villela DAM, et al. Climate information for arbovirus risk monitoring: opportunities and challenges. Bull Am Meteorol Soc. 2016;97:ES107–ES111.

18. Majeed H, Coles JG, Moore GWK. Influence of Atlantic and Pacific sea surface temperatures on heat-related mortality in the United States. Geohealth. 2020;4:e2019GH000220.

19. Morris JS. Functional regression. Annu Rev Stat Its Appl. 2015;2:321-359.

20. Brockhaus S, Scheipl F, Hothorn T, Greven S. The functional linear array model. Stat Model. 2015;15:279–300.

21. Bühlmann P, Hothorn T. Boosting algorithms: regularization, prediction and model fitting. Stat Sci. 2007;22:477–505.

22. Thornton MM, Shrestha R, Wei Y, Thornton PE, Kao S-C, Wilson BE.

Daymet: Daily Surface Weather Data on a 1-km Grid for North America, Version 4 . ORNL DAAC[Internet] . December 15, 2020 [cited 2021 May 3]. Available at:

https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1840 . Accessed 3 May 2021.

23. Thornton PE, Shrestha R, Thornton M, Kao SC, Wei Y, Wilson BE. Gridded daily weather data for North America with comprehensive uncertainty quantification. Sci Data. 2021;8:190.

24. Assani AA, Maloney-Dumont V, Pothier-Champagne A, Kinnard C, Quéssy J-F. Comparison of the temporal variability of summer temperature and rainfall as it relates to climate indices in southern Quebec (Canada). Theor Appl Climatol. 2019;137:2425–2435.

25. Gasparrini A, Armstrong B, Kenward MG. Distributed lag non-linear models. Stat Med. 2010;29:2224–2234.

26. Gasparrini A, Guo Y, Hashizume M, et al. Temporal variation in heat-mortality associations: A Multicountry Study. Environ Health Perspect. 2015;123:1200–1207.

27. Gasparrini A, Leone M. Attributable risk from distributed lag models. BMC Med Res Methodol. 2014;14:55.

28. Guo Y, Gasparrini A, Armstrong BG, et al. Heat wave and mortality: A Multicountry, Multicommunity Study. Environ Health Perspect. 2017;125:087006.

29. Gasparrini A, Armstrong B. Reducing and meta-analysing estimates from distributed lag non-linear models. BMC Med Res Methodol. 2013;13:1.

30. Tobías A, Hashizume M, Honda Y, et al. Geographical variations of the minimum mortality temperature at a global scale: A Multicountry Study. Environ Epidemiol. 2021;5:e169.

31. Ramsay JO, Silverman BW. Functional Data Analysis. 2nd ed. Wiley Online Library; 2005.

32. Eilers PHC, Marx BD. Flexible smoothing with $B$-splines and Penalties. Stat Sci. 1996;11:89–102.

33. Hastie T, Tibshirani R, Friedman J. The Elements of Statistical Learning. 2nd ed. Springer with Series in Statistics; 2009.

34. R Core Team (2021). R:

A language and environment for statistical computing . R Foundation for Statistical Computing, Vienna, Austria. URL

https://www.R-project.org/ .

35. Gasparrini A. Distributed lag linear and non-linear models in R: the package dlnm. J Stat Softw. 2011;43:1–20.

36. Brockhaus S, Melcher M, Leisch F, Greven S. Boosting flexible functional regression models with a high number of functional historical effects. Stat Comput. 2017;27:913-926. 2016.

37. Schlesinger ME, Ramankutty N. An oscillation in the global climate system of period 65–70 years. Nature. 1994;367:723–726.

38. Thiombiano AN, El Adlouni S, St-Hilaire A, Ouarda TBMJ, El-Jabi N. Nonstationary frequency analysis of extreme daily precipitation amounts in Southeastern Canada using a peaks-over-threshold approach. Theor Appl Climatol. 2017;129:413–426.

39. Enfield DB, Mestas-Nuñez AM, Trimble PJ. The Atlantic Multidecadal Oscillation and its relation to rainfall and river flows in the continental U.S. Geophys Res Lett. 2001;28:2077–2080.

40. Doyon B, Bélanger D, Gosselin P. The potential impact of climate change on annual and seasonal mortality for three cities in Québec, Canada. Int J Health Geogr. 2008;7:23.

41. Curriero FC, Heiner KS, Samet JM, Zeger SL, Strug L, Patz JA. Temperature and mortality in 11 cities of the eastern United States. Am J Epidemiol. 2002;155:80–87.

42. Bayentin L, El Adlouni S, Ouarda TB, Gosselin P, Doyon B, Chebana F. Spatial variability of climate effects on ischemic heart disease hospitalization rates for the period 1989-2006 in Quebec, Canada. Int J Health Geogr. 2010;9:5.

43. Sera F, Griffiths LJ, Dezateux C, Geraci M, Cortina-Borja M. Using functional data analysis to understand daily activity levels and patterns in primary school-aged children: cross-sectional analysis of a UK-wide study. PLoS One. 2017;12:e0187677.

44. Arisido MW. Functional measure of ozone exposure to model short-term health effects. Environmetrics. 2016;27:306–317.

45. Masselot P, Chebana F, Ouarda TBMJ, Bélanger D, St-Hilaire A, Gosselin P. A new look at weather-related health impacts through functional regression. Sci Rep. 2018;8:15241.

46. Curceac S, Ternynck C, Ouarda TBMJ, Chebana F, Niang SD. Short-term air temperature forecasting using Nonparametric Functional Data Analysis and SARMA models. Environmental Modelling & Software. 2019;111:394-408.

47. Chebana F, Martel B, Gosselin P, Giroux JX, Ouarda TB. A general and flexible methodology to define thresholds for heat health watch and warning systems, applied to the province of Québec (Canada). Int J Biometeorol. 2013;57:631–644.

48. Toutant S, Gosselin P, Bélanger D, Bustinza R, Rivest S. An open source web application for the surveillance and prevention of the impacts on public health of extreme meteorological events: the SUPREME system. Int J Health Geogr. 2011;10:39.

49. Lee M, Nordio F, Zanobetti A, Kinney P, Vautard R, Schwartz J. Acclimatization across space and time in the effects of temperature on mortality: a time-series analysis. Environ Health. 2014;13:89.

50. Assani AA, Guerfi N. Analysis of the joint link between extreme temperatures, precipitation and climate indices in winter in the three hydroclimate regions of Southern Quebec. Atmosphere. 2017;8:75.