Attributable risk from distributed lag models

Attributable risk measures

with n as the total number of cases. The parameter β _x used in Eq. (1a) represents the risk associated with the exposure, and it usually corresponds to the logarithm of a ratio measure such as relative risk, relative rate or odds ratio. It is generally obtained from regression models while adjusting for potential confounders. The general definition of β _x used here refers to the association with a specific exposure intensity x compared to a reference value x ₀. For linear exposure-response relationships, the association can also be reported as β·x, where in this case β refer to a unit increase in x. For binary variables reporting presence/absence of the exposure, Eq. (1a) simplifies to AF=(RR−1)/RR, with RR as relative risk, as reported by Steenland and Armstrong [2]. We keep the more general definition of β _x , which is easily applicable to non-linear exposure-response relationships, throughout the manuscript.

The theoretical nature of these effect measures is based on a counterfactual, where the observed condition is compared with a reference state which never occurred. This state postulates that the same population is followed in an identical situation where only the exposure level changes to the reference value x ₀. Typically, such a reference is represented by the absence of association, meaning x ₀=0 and ${\beta }_{x_0}=0$. However, different counterfactual conditions can be used, for example a lower exposure which can be determined by an intervention. In this case the quantity β _x can be simply re-parameterized as ${\beta }_x^{\ast }={\beta }_x-{\beta }_{x_0}$, and Eq. (1) still applies.

with ${\text{AN}}_{x_1,\dots ,x_p}$ obtained by substituting Eq. (2) in Eq. (1b) [2]. For the specific form of Eq. (2), it should be noted that ${\text{AF}}_{x_1,\dots ,x_p}\le {\text{AF}}_{x1}+\dots +{\text{AF}}_{\mathit{\text{xp}}}$, i.e. the sum of the attributable risk measured for individual exposures is usually higher than their concurrent attributable risk.

A review of the DLNM modelling framework

The function s(x,t) is computed as the approximate integral of the exposure-lag-response function over the lag dimension, representing the cumulated risk over the lag period. The parameterization in the final step of Eq. (3) is obtained through a cross-basis, involving a tensor product between the basis chosen for f(x) and w(ℓ), generating the transformed variables w _x,t linearly combined with the parameters η. Simpler DLMs are defined by Eq. (3) by assuming f(x) as linear. Algebraic details and additional information are provided elsewhere [11]. The cross-basis is specified with a reference value x ₀ used later as a centering point for the function f(x), which is used to define the counterfactual condition.

This overall cumulative association is composed of the sum of contributions β _x,ℓ from exposures $x_{t-{\ell }_0},\dots ,x_{t-L}$ experienced within the lag period. Algebraic definitions have been previously provided [11].

Forward and backward perspectives

The term β _x,ℓ for each intensity x can be interpreted using two complementary perspectives, illustrated graphically in Figure 1. From a forward standpoint, looking from current exposure to future risks, the terms β _x,ℓ are the contributions from the exposure x _t occurring at time t to the risk at times t+ℓ ₀,…,t+L, identified by green circles. From a backward standpoint, looking from current risk to past exposures, the terms β _x,ℓ are the contributions to the risk at time t from exposures $x_{t-{\ell }_0},\dots ,x_{t-L}$ experienced at t−ℓ ₀,…,t−L, identified by yellow squares. The underlying curve in Figure 1 depicting such associations is called the lag-response curve related to a given exposure intensity x. The sum of these contributions over the whole lag period can be interpreted as the overall cumulative risk.

Attributable risk from DLNMs

with n _t as the number of cases at time t. This structure is consistent with the configuration of the regression model usually applied to fit the data, where the risk at time t is associated with lagged exposures at times t−ℓ. The definition of backward attributable risk requires an extended version of the counterfactual condition accounting for the additional lag dimension: b-AN _x,t and b-AF _x,t are interpreted as the number of cases and the related fraction at time t attributable to past exposures to x in the period t−ℓ ₀,…,t−L, compared to a constant exposure x ₀ throughout the same period.

This alternative version has some advantages if compared to the backward definition. First, the counterfactual condition is simpler: f-AF _x,t and f-AN _x,t are interpreted as the fraction and number of future cases in the period t+ℓ ₀,…,t+L attributable to the single exposure x occurring at time t, compared to x ₀. Moreover, the overall cumulative risk $\sum {\beta }_{x_t,\ell }$ for a given exposure x _t in (6a) is available also when the bi-dimensional exposure-lag-response is reduced to uni-dimensional exposure-response relationship, a step often needed in multi-site studies [14]. In contrast, all the lag-specific contributions are needed to compute $\sum {\beta }_{x_{t-\ell },\ell }$ in (5a) for the backward counterpart.

However, the forward version also has an important limitation, related to the fact that the contributions are associated to risks measured at different times. The attributable number f-AN _x,t in (6b) is computed by averaging the total counts experienced in the next ℓ ₀,…,L times, thus only approximating the lag structure of risks. This approximation is likely to produce some bias, which is expected as an underestimation of the attributable number if compared to the backward version.

Separating attributable components

simply selecting the risk contributions from past exposures included in the range r. The related attributable number ${\mathrm{b}-\mathit{\text{AN}}}_{x,t}^r$ is computed by substituting Eq. (7) into Eq. (5b). Attributable components referring to different ranges can be summed up, as all are defined using the same counterfactual condition of a constant exposure x _ℓ =x ₀ for the whole lag period ℓ=ℓ ₀,…,L.

The forward version has the additional advantage that for two non-overlapping ranges r ₁ and r ₂ the sum of the components is equal to the overall attributable risk, namely ${\mathrm{f}-\mathit{\text{AF}}}_{x,t}^{r_1+r_2}={\mathrm{f}-\mathit{\text{AF}}}_{x,t}^{r_1}+{\mathrm{f}-\mathit{\text{AF}}}_{x,t}^{r_2}$. In contrast, adopting a backward perspective ${\mathrm{b}-\mathit{\text{AF}}}_{x,t}^{r_1+r_2}\le {\mathrm{b}-\mathit{\text{AF}}}_{x,t}^{r_1}+{\mathrm{b}-\mathit{\text{AF}}}_{x,t}^{r_2}$, as the risks are simultaneously computed for the same time t in the like of Eq. (2).

Total attributable risk

The equations above can be applied either to forward or backward attributable risk and to separate components, simply substituting the related attributable numbers in Eq. (8a).

Computing uncertainty intervals

Analytical formulae for confidence intervals of attributable risk measures are not easily produced [15], and this also applies to the extended versions developed here. Although approximated estimators have been proposed [15, 16], in this context the most straightforward approach is to rely on interval estimation obtained empirically through Monte Carlo simulations [17, 18]. Basically, we take random samples η ^(j) of the original parameters η of the cross-basis in Eq. (3) from the assumed multivariate normal distribution with point estimate $\widehat{\mathit{\eta }}$ and (co)variance matrix $V\left(\widehat{\mathit{\eta }}\right)$ derived from the regression model. These samples η ^(j) are used to compute ${\beta }_{x,\ell }^{\left(j\right)}$ for ℓ=ℓ ₀,…,L and each intensity x, empirically reconstructing the distributions of the attributable measures defined in Eq. (5)–(8). The related 2.5^th and 97.5^th percentiles of such distributions are interpreted as 95% empirical confidence intervals (eCI).

Modelling strategy

We fitted a standard time series Poisson model allowing for overdispersion, controlling for seasonal and long term trends and day of the week, using a 10 df/year spline and indicator variables, respectively. Model selection is still an issue of current research within the DLNM framework, although simulation studies indicate a good performance of methods based on the Akaike Information criterion (AIC) [11]. Considering the illustrative purpose of the example, we selected a-priori the cross-basis function in Eq. (3) for representing the association between mean daily temperature and mortality, basing our choice on previous analyses. Specifically, the cross-basis is composed of a quadratic B-spline with two equally-spaced knots as the exposure-response function f(x), and a natural cubic B-spline with three equally-spaced knots in the log-scale as the lag-response function w(ℓ) over lags 0–25.

In the specific case of temperature where a null exposure condition cannot be defined, a reasonable choice is to center the cross-basis in Eq. (3) to the temperature of minimum risk, as suggested in previous publications [13]. This optimal temperature corresponds to 20°C and 21°C for London and Rome respectively, and it represents the reference point x ₀ for the computation of the attributable risk measures. These are obtained for the whole temperature range, and then for cold and heat contributions by separating the associations with temperatures lower or higher than x ₀. In addition, the attributable components are separated further in mild and extreme cold and heat by selecting as cut-off values the 1^st and 99^th percentiles of city-specific distributions, corresponding to 0.4°C and 23.7°C in London and 2.6°C and 28.6°C in Rome.

We derived empirical confidence intervals for backward total attributable numbers and fractions, computed overall and for separated components, by simulating 5,000 samples from the assumed distribution of $\widehat{\mathit{\eta }}$.

Risk attributable to temperature

The estimated associations between temperature and all-cause mortality in the two cities are illustrated in Figure 2. The left panels show the bi-dimensional exposure-lag-response surfaces, while the right panels display the overall cumulative exposure-response curves, interpreted as the risk cumulated over the entire lag period of 0–25 days. The associations are represented in the RR scale, with the centering point and the cut-off values for defining extreme cold and heat displayed as dotted and dashed vertical lines respectively.

Cold, heat and extreme components

The total backward attributable risk is then separated into components due to cold and hot temperatures, defined as those below and above the optimal temperature, respectively. The estimates, computed using Eq. (7), are reported in Table 1. The comparison of the two contributions clearly indicates that cold is responsible for most of the mortality attributable to temperature, with b-AF _{t o t} equal to 12.95% and 10.84%, compared to 0.66% and 1.74% for heat, in the two cities. Estimates of forward attributable risk are very similar, and as expected their sum is equal to the overall burden, differently than for the backward version.

In contrast, the comparison between the two cities is rather different for the components attributable to mild and extreme hot temperatures. In spite of the stronger risk in London, the attributable fraction is similar for extreme heat and even higher in Rome for mild heat (1.32%–1.45% versus 0.25%–0.33%). This apparent contradiction is explained by the different temperature distribution, and in particular the percentile corresponding to the optimal temperature, corresponding to 93.6^th and 72.5^th in London and Rome. This result suggests the hypothesis that the population in Rome is more adapted to the range of temperatures corresponding to extreme hot if compared to London, where the population experienced only a few days of unusually high temperatures.

The harvesting paradox

Accounting for the additional lag dimension in exposure-lag-response associations involves further complexities in the interpretation of attributable risk measures. We now focus our attention to the association with hot temperature in Rome. The left panel of Figure 3 shows the estimated lag-response curves at various temperatures, computed versus the reference optimal temperatures of 21°C. The graph indicates a strong risk in the first lags, then followed by a protective association at longer lags. This pattern is consistent with the harvesting hypothesis, which assumes that the initial risk is partly discounted by the depletion of the pool of susceptible after an extreme heat event [19]. Also, the extent of harvesting seems more pronounced for extreme temperatures.

This phenomenon has interesting implications. An example is offered by the right panel of Figure 3, illustrating the estimated daily deaths b-AN _x,t and f-AN _x,t attributable to heat, computed backward and forward for the first summer in the time series for Rome, with related temperature trend. As expected, a substantial number of deaths are attributable to temperatures above the optimal value, represented by the horizontal dotted line, in the period mid-July to mid-August. The trend of forward attributable deaths f-AN _x,t closely follows the daily temperatures, consistently with the definition of number of deaths attributable to the temperature in day t cumulated in the next L days. In contrast, the backward attributable number b-AN _x,t decreases to zero and even becomes negative in late summer days, although the overall cumulative exposure-response in Figure 2 (bottom-right panel) does not show a RR below 1 for any temperature.

This paradox is explained by the counterfactual condition associated with the backward perspective. Specifically, each b-AN _x,t compares the association with the observed temperatures in the past L days to a constant exposure x ₀. In the presence of harvesting, the observed population becomes ‘healthier’ than the counterfactual population after a series of heat days, due to the depletion of the susceptible pool. This explains the negative attributable numbers for specific combinations of lagged exposures. This fact emphasises that harvesting should not be interpreted as a true protective association at longer lags, but rather as an artefact due to a change in the underlying population following a stress, which affects the counterfactual condition. This issue is relevant when using backward attributable risk measures b-AN _x,t and b-AF _x,t to assess the contribution of specific days. However, similarly to the net overall cumulative risk, the total attributable number b-AN _{t o t} and fraction b-AN _{t o t} , produced by Eq. (8) and reported in Tables 1 and 2, account for the discount by summing the contributions over the whole series.