Relationships between environmental factors and various health outcomes are sometimes characterized by lag patterns,1–4 and epidemiologists have proposed a number of modelling approaches to account for this additional complexity in defining exposure–response associations.5–8 In particular, these methods are required in the presence of long lag periods, which have been previously reported in associations with cold temperatures3,9,10 and, to a lesser extent, air pollutants.11–14 However, there is little evidence on the comparative performance of alternative models, or on their sensitivity to assumptions about the shape of the lag structure or the length of the lag interval. In addition, recent published articles, specifically evaluating the health burden of cold temperature, have contended that part of the estimated associations can be due to the confounding effect of season, if long lags are assumed.15,16 On the other hand, no empirical or well-grounded theoretical evidence is provided in support of these statements.
Clearly, the problem of modelling such phenomena has important implications in environmental epidemiology, and more generally in biomedical research. This contribution aims to address the issue through a simulation study on the comparative performance of two alternative methods for investigating lagged associations.
MODELLING LAGGED ASSOCIATIONS
In environmental time series studies, the outcome, commonly daily mortality or morbidity counts, is usually compared with levels of environmental exposures in a Poisson regression allowing for overdispersion.17 This model represents the counts of outcome events Yt at day t as
;)
The function g(t), usually a spline whose flexibility is defined by the number of degrees of freedom per year, controls for seasonal and long-term trends. Additional potential confounders Zp showing day-to-day variability are controlled for through functions Vp. The function S describes the association with the exposure of interest x measured over the lag interval
, with L as maximum lag.
Two main methods have been proposed to define
.
Moving Average Models
A simple approach relies on averaging the exposure over the lag interval. For linear exposure–response relationships, the parameterization reduces to the moving average of lagged exposures. Adopting some liberty in the algebraic notation, here and in equations below, with the exclusion of the model coefficients, the moving average can be computed as
However, the definition is not straightforward when nonlinear dependencies are assumed, and modeled through an exposure–response function f(x). In this case, the moving average can be defined in two alternative forms:
The form in Equation (3.1) applies the function f to the moving average of x, while the form in Equation (3.2) is the moving average of the values transformed by f. The latter is defined over the whole range of the observed data. The former, although being defined in a narrower range of x, as the extreme values are averaged out, is simpler to compute and frequently applied in environmental studies.9,15,18,19
Distributed Lag Linear and Nonlinear Models
In an alternative and more advanced approach, the lag structure can be explicitly modeled through a lag-response function
expressed in the lag dimension
, with the function S defined as
This parameterization, applicable to describe linear lagged associations, is known as a distributed lag model, and simplifies to the moving average in Equation 2 when
is a constant function. In its original development for econometric time series,20 and the following application in environmental epidemiology,5 a polynomial function was used to specify
. A simpler alternative is represented by an unconstrained distributed lag model, where lag-specific effects are modeled by indicators for each lag. Further developments of the framework are detailed in recent publications. Specifically, these illustrate the use of alternative functions for modelling the lag response,21,22 the extension to nonlinear exposure–response relationship,21,22 and the generalization beyond time series data.23 In particular, the nonlinear extension to distributed lag nonlinear models is obtained by expressing the association through the combination of the exposure–response and lag–response functions. This combination produces a bidimensional exposure–lag–response function
. In this case S, termed cross-basis function, is represented as
DLNMs simplifies to the linear counterpart in Equation 4 when f(x) is assumed linear, and to the second moving average form defined in Equation (3.2) when the lag-response function
is a constant term. In this contribution, I will refer generally to distributed lag models, pointing out if linear or nonlinear when required. These models have been recently applied to study health effects of environmental factors.12–14,24,25
SIMULATION SETTINGS
These alternative methods were comparatively assessed by simulating associations between air pollution and temperature with all-cause mortality, using real daily time series data from Chicago in the period 1987–2000.26 The data represent aggregated series and are publicly available, and therefore no ethical review was required. The simulated exposures were derived as the ozone series, standardized over the range 0 to 50 part per billion (ppb), and the temperature series, standardized over the range −20 to 35°C. The baseline mortality was obtained as the predicted curve from the regression model fitting a natural cubic spline of time with 10 df/year to the real mortality series, thus simulating the complex seasonal trend represented in Figure 1.
FIGURE 1: Baseline mortality trend as death counts showing the complex seasonal and long-term patterns simulated using real daily time series from Chicago 1987–2000.
The excess risk associated with the exposures was defined by scenarios of linear ozone-mortality and nonlinear temperature-mortality dependencies, assuming either short or long-lagged associations. Specifically, the following scenarios were generated:
- Linear short-lag scenario (scenario 1): ozone is associated linearly with mortality over a short lag period, limited to 0–3 days.
- Linear long-lag scenario (scenario 2): ozone is associated linearly with mortality at longer lags, with the excess risk smoothly decreasing over the following 20 days.
- Nonlinear short-lag scenario (scenario 3): both hot and cold temperatures are associated nonlinearly with mortality over a short lag period, limited to 0–3 days.
- Nonlinear long-lag scenario (scenario 4): temperature shows a nonlinear relationship, with hot temperatures associated only with short lags as in scenario 3, and cold temperatures associated at longer lags, extended for up to 20 days with a peak at lag 4.
The four simulated exposure–lag–response surfaces, compatible with estimates previously reported in the literature,12,14,24,25,27 are illustrated in Figure 2. The mortality series in each replicate was simulated assuming a Poisson distribution with expectation equal to the baseline counts multiplied by the excess due to the risk of the exposure cumulated over the lag period. Technical details are provided in eAppendix 1 (https://links.lww.com/EDE/B80).
FIGURE 2: Simulated exposure–lag–response surfaces as relative risk in four scenarios describing the linear association between ozone and mortality and the nonlinear association between temperature and mortality, each characterized by short and long-lag patterns. The bold lines represent lag–response relationships at 10 ppb of ozone and at −15 and 30°C of temperature.
Performance and inferential properties of the two approaches were compared by fitting regressions with the same splines of time used when simulating the data, and the following models for the associations with the exposures:
- Three moving average models, computed using the forms in Equation 2 (in scenarios 1 and 2) and in Equation (3.2) (in scenarios 3 and 4), with lag periods of 0–3, 0–7, and 0–20, respectively.
- A distributed lag linear (in scenarios 1 and 2) and nonlinear (in scenarios 3 and 4) model with the lag-response function
- specified by a natural cubic spline with three knots equally spaced in the log scale, plus intercept, over lag 0–20.
The nonlinear exposure–response function f(x) for both the moving average models and the distributed lag nonlinear model in scenarios 3 and 4 was specified by a quadratic B-spline with three knots at the 10th, 75th, and 90th percentiles of temperature distribution. These choices follow models applied in previous work.24
A set of sensitivity analyses was carried out to test the impact of alternative specifications of the function
in distributed lag models (with details provided below), and differences between moving average models adopting the two forms of moving average in Equation (3). In addition, the impact of overdispersion was assessed by simulating data from a negative binomial distribution with overdispersion parameter
.
The simulation results are based on 5,000 replicates, and are summarized in terms of bias (difference between the true simulated log-relative risk [log-RR] and average of its estimates across the replicates), coverage (percentage of times the confidence interval of the estimates includes the true log-RR), and root mean square error (RMSE, average across replicates of the squared difference between the true simulated and estimated log-RR). The latter can be interpreted as the sum of the bias and the imprecision of the estimator. These statistics are reported as the average across the curve representing the exposure–response association cumulated over the lag period. An algebraic definition is provided in eAppendix 2 (https://links.lww.com/EDE/B80).
The R scripts illustrating a simple example of simulation of the data and fit of the models, and then fully reproducing the results, are available at www.ag-myresearch.com.
RESULTS
Figures 3 and 4 demonstrate the good performance of the distributed lag models in recovering the lag–response curves corresponding to 10 ppb of ozone and −15 and 30°C of temperature, respectively, simulated in the four scenarios and also represented as bold black lines in Figure 2. These findings are confirmed by Figure 5, which reports the average estimates of the overall cumulative exposure–response relationships. The graphs suggest that distributed lag models provide nearly unbiased estimates of the associations.
FIGURE 3: Lag–response relationships as relative risk at 10 ppb of ozone in the linear short-lag (left column) and linear long-lag (right column) scenarios, estimated by the distributed lag model. The curves represent the real simulated relationship (continuous black lines), the average of the estimated relationship across 5,000 replicates (dashed red lines), and a sample of the estimated relationship in 20 replicates (continuous grey lines), respectively.
FIGURE 4: Lag–response relationships as RR at 30°C (top row) and −15°C (bottom row) in the nonlinear short-lag (left column) and nonlinear long-lag (right column) scenarios, estimated by the distributed lag non-linear model (DLNM). The curves represent the real simulated relationship (continuous black lines), the average of the estimated relationship across 5,000 replicates (dashed red lines), and a sample of the estimated relationship in 20 replicates (continuous grey lines), respectively.
FIGURE 5: Overall cumulative exposure–response relationships as RR estimated by the distributed lag models and two moving average models with lag 0–3 and 0–20 (MA0–3 and MA0–20). The graphs represent the linear association with ozone (top panels) and the nonlinear association with temperature (bottom panels) in the short-lag (left panels) and long-lag (right panels) scenarios. The slopes represent the real simulated relationship (continuous black lines) and the average of the estimated relationship across 5,000 replicates (dashed lines, with different colours and patterns depending on the model). MA indicates moving average.
The panels of Figure 5 report also comparable estimates for moving average models with lag periods 0–3 and 0–20. The results suggest that the lag 0–3 model shows little bias in compatible short-lag scenarios (left panels), but, as expected, it considerably under estimates the excess risk when long lagged associations occur, namely with ozone in the top-right panel and with cold temperature in the bottom-right panel of Figure 5. The reverse occurs to the lag 0–20 model, with the left panels of Figures 5 indicating a strong downward bias in moving average models when the lag interval is extended further than the simulated lag period. Interestingly, in the nonlinear long-lag scenario, moving average models show biases even when the lag period is correctly specified. The moving average model with lag 0–7 shows an intermediate behavior, as illustrated in eFigure 1 (https://links.lww.com/EDE/B80). Results in eFigure 2 (https://links.lww.com/EDE/B80) indicate that the two forms of moving average defined in Equations (3.1)–(3.2) provide almost identical estimates, although the latter is defined over the whole exposure range. The bias generated by the strong assumptions of moving average models about the shape of the exposure–lag–response risk surface, compared with the flexibility of distributed lag models, are illustrated in eFigure 3 (https://links.lww.com/EDE/B80).
Table 1 complements the illustration of the findings, comparing the inferential performance of the four models. Distributed lag models show the lowest bias and nominal coverage in both scenarios. In terms of RMSE, the higher precision favors the moving average model with lag 0–3 in the short-lag scenario, but this is counterbalanced in the long-lag scenario by the bias previously discussed. The extension of the lag interval of the moving average to 0–20 shows no improvement, as this model provides less precise estimates without eliminating the bias. All moving average models are also affected to some extent by under coverage, which is very pronounced in scenarios with nonlinear relationships.
TABLE 1: Bias, Coverage of 95% Confidence Intervals, and RMSE, Reported for the DLM and DLNM and for MA Models with Different Lag Intervals (0–3, 0–7, and 0–20), in Four Alternative Simulated Scenarios
Results of the sensitivity analysis are summarized in Table 2, and can be compared with the main results reported in Table 1. The simplification of the function
of distributed lag models in short-lag scenarios, either by using an unconstrained parameterization over lag 0–3 or by reducing the number of knots to 2 over the lag period 0–7, substantially increase the precision of the estimates, with RMSE comparable with the moving average lag 0–3 models (analyses 1–4). The different specification of
in long-lag scenarios, using two or four knots, does not seem to affect the estimates (analyses 5–6). The presence of overdispersion, if modeled through a quasi-Poisson family, slightly decrease precision but does not introduce biases in distributed lag models (analysis 7). Results of the sensitivity analyses are reported graphically in eFigure 4 (https://links.lww.com/EDE/B80).
TABLE 2: Bias, Coverage of 95% Confidence Intervals, and RMSE, Reported for distributed lag model and distributed lag non-linear model (DLNM), in Sensitivity Analyses Using Different Model Specifications in Various Scenarios
DISCUSSION
This simulation study compares alternative approaches to model epidemiologic associations characterized by lag patterns. Results indicate that standard methods based on moving average models can introduce important biases, in particular in the presence of extended lag periods. In contrast, more flexible methods based on distributed lag models appropriately account for lagged dependencies, with no or minimal bias in point estimates and confidence intervals.
The framework of distributed lag models has the advantage of offering a direct representation of the lag structure, through lag–response relationships estimated by the data. This extra dimension provides additional information on the phenomenon under study. Results in Figures 3–5 and Table 1 demonstrate that these models can correctly retrieve the underlying association, summarized as overall cumulative (net) risk or as lag-specific contributions, even when the lag interval is extended well beyond the true lag period. Simpler methods based on moving average models represent a viable alternative only in the presence of relatively short lag periods and when the lag interval is correctly specified. Conversely, the application of moving average models to roughly approximate long and complex lag patterns, or the extension of the moving average beyond the actual lag interval, can result in substantial biases, in particular when modelling nonlinear relationships.
The inflexible specification of the lag structure is the main reason why moving average models are outperformed by distributed lag models. As mentioned above, the former can be interpreted as special cases of the latter in which the lag–response function
is defined as a constant. This definition implies the strong assumption that all the exposures experienced within the lag interval equally contribute to the overall cumulative excess risk, as demonstrated in eFigure 3 (https://links.lww.com/EDE/B80). In contrast, flexible specifications within the distributed lag modelling framework allow more realistic lag structures, when for instance more recent exposures are associated with higher lag-specific excess risks if compared with more lagged exposures.
In the simple scenario describing a situation in which the lag interval is short and correctly specified, moving average models benefit from a higher precision in the estimates, which can avoid artifactual interpretation of the results. However, it should be noted that alternative and more efficient specifications of the lag–response function
, obtained for instance by reducing the lag interval and/or by selecting simpler functions, are available within the distributed lag modelling framework. These alternative choices for the lag–response function, illustrated in the sensitivity analysis, improve the performance of the distributed models in simpler scenarios by increasing the precision, without introducing biases in point estimates and confidence intervals. It should be acknowledged, however, that the flexibility of the distributed framework presents the additional complexity of selecting bidimensional exposure–lag–response functions, if compared with unidimensional moving average models. The choice can be based on evidence from previous studies, as in the examples above, or on existing selection methods (such as the Akaike information criterion) proposed in methodologic publications.22,23 Furthermore, model selection in distributed lag models is an issue of current research, and future developments are likely to provide improved selection criteria.
An important result of this simulation study is the provision of empirical evidence that the flexible modelling of the lag dimension is not affected by the presence of complex seasonal trends. These findings contradict claims, recently suggested in the literature,15,16 that increased mortality risks attributed to lagged associations with cold temperature are partly due to confounding effects by season. This simulation study demonstrates that complex exposure–lag–response associations can be effectively disentangled from strong seasonal trends.
In conclusion, this study provides some guidance on the use of alternative methods to model lagged associations, identifying limitations of standard approaches and advantages of recent and more flexible methods.
REFERENCES
1. Zanobetti A, Schwartz J, Samoli E, et al. The temporal pattern of mortality responses to air pollution: a multicity assessment of mortality displacement. Epidemiology. 2002;13:87–93.
2. Goodman PG, Dockery DW, Clancy L. Cause-specific mortality and the extended effects of particulate pollution and temperature exposure. Environ Health Perspect. 2004;112:179–185.
3. Carder M, McNamee R, Beverland I, et al. The lagged effect of cold temperature and wind chill on cardiorespiratory mortality in Scotland. Occup Environ Med. 2005;62:702–710.
4. Pope CA III. Mortality effects of longer term exposures to fine particulate air pollution: review of recent epidemiological evidence. Inhal Toxicol. 2007;19(Suppl 1):33–38.
5. Schwartz J. The distributed lag between air pollution and daily deaths. Epidemiology. 2000;11:320–326.
6. Dominici F, McDermott A, Zeger SL, Samet JM. Airborne particulate matter and mortality: timescale effects in four US cities. Am J Epidemiol. 2003;157:1055–1065.
7. Roberts S. Using moving total mortality counts to obtain improved estimates for the effect of air pollution on mortality. Environ Health Perspect. 2005;113:1148–1152.
8. Muggeo VM, Hajat S. Modelling the non-linear multiple-lag effects of ambient temperature on mortality in Santiago and Palermo: a constrained segmented distributed lag approach. Occup Environ Med. 2009;66:584–591.
9. Anderson BG, Bell ML. Weather-related mortality: how heat, cold, and heat waves affect mortality in the United States. Epidemiology. 2009;20:205–213.
10. Analitis A, Katsouyanni K, Biggeri A, et al. Effects of cold weather on mortality: results from 15 European cities within the PHEWE Project. Am J Epidemiol. 2008;168:1397.
11. Bell ML, McDermott A, Zeger SL, Samet JM, Dominici F. Ozone and short-term mortality in 95 US urban communities, 1987-2000. JAMA. 2004;292:2372–2378.
12. Samoli E, Zanobetti A, Schwartz J, et al. The temporal pattern of mortality responses to ambient ozone in the APHEA project. J Epidemiol Community Health. 2009;63:960–966.
13. Zanobetti A, Schwartz J, Samoli E, et al. The temporal pattern of respiratory and heart disease mortality in response to air pollution. Environ Health Perspect. 2003;111:1188–1193.
14. Fraga J, Botelho A, Sá A, Costa M, Quaresma M. The lag structure and the general effect of ozone exposure on pediatric respiratory morbidity. Int J Environ Res Public Health. 2011;8:4013–4024.
15. Kinney PL, Schwartz J, Pascal M, et al. Winter season mortality: will climate warming bring benefits? Environ Res Lett. 2015;10:064016.
16. Woodward A. Heat, cold and climate change. J Epidemiol Community Health. 2014;68:595–596.
17. Bhaskaran K, Gasparrini A, Hajat S, Smeeth L, Armstrong B. Time series regression studies in environmental epidemiology. Int J Epidemiol. 2013;42:1187–1195.
18. Chung Y, Lim YH, Honda Y, et al. Mortality related to extreme temperature for 15 cities in northeast Asia. Epidemiology. 2015;26:255–262.
19. Todd N, Valleron AJ. Space-time covariation of mortality with temperature: a systematic study of deaths in France, 1968-2009. Environ Health Perspect. 2015;123:659–664.
20. Almon S. The distributed lag between capital appropriations and expenditures. Econometrica. 1965;33:178–196.
21. Armstrong B. Models for the relationship between ambient temperature and daily mortality. Epidemiology. 2006;17:624–631.
22. Gasparrini A, Armstrong B, Kenward MG. Distributed lag non-linear models. Stat Med. 2010;29:2224–2234.
23. Gasparrini A. Modelling exposure-lag-response associations with distributed lag non-linear models. Stat Med. 2014;33:881–899.
24. Gasparrini A, Guo Y, Hashizume M, et al. Mortality risk attributable to high and low ambient temperature: a multicountry observational study. Lancet. 2015;386:369–375.
25. Huang Z, Lin H, Liu Y, et al. Individual-level and community-level effect modifiers of the temperature-mortality relationship in 66 Chinese communities. BMJ Open. 2015;5:e009172.
26. Samet JM, Zeger SL, Dominici F, Curriero F, Coursac I, Dockery DW. The National Morbidity, Mortality, and Air Pollution Study (NMMAPS). Part 2. Morbidity and Mortality from Air Pollution in the United States. Health Effects Institute. 2000.
27. Guo Y, Barnett AG, Pan X, Yu W, Tong S. The impact of temperature on mortality in Tianjin, China: a case-crossover design with a distributed lag nonlinear model. Environ Health Perspect. 2011;119:1719–1725.
