Motivation

• Climate change is intensifying extreme environmental conditions: heatwaves, cold spells are becoming more frequent and severe

• Environmental exposures rarely act instantaneously, as their health impact unfolds over time:

  • A heatwave today may elevate mortality risk several days later

  • Cold temperatures can trigger cardiovascular events long after exposure

  • Air pollution accumulates biological damage across days or weeks

• Standard regression approaches struggle to capture these delayed and usually non-linear effects simultaneously

Distributed Lag Non-Linear Models (DLNMs)

DLNMs (Gasparrini et al. 2010) are the standard framework for capturing these dynamics, simultaneously modelling two dimensions of risk:

• Usually splines (natural cubic, B-splines or P-splines) are used to define these usually non-linear relationships in each dimension.

Distributed Lag Non-Linear Models (DLNMs)

• The tensor product of the two spline bases defines the cross-basis: a smooth bivariate surface over exposure and lag

DLNMs: Model formulation

\[\log(E(Y_t)) = \alpha + \underbrace{cb(x_t, \ldots, x_{t-L}) \cdot \boldsymbol{\beta}}_{\text{cross-basis}} + \sum_{k} \gamma_k u_{kt}\]

Bayesian DLNMs

Same model, with prior distributions placed on all parameters:

\[\log(E(Y_t)) = \alpha + \underbrace{cb(x_t, \ldots, x_{t-L}) \cdot \boldsymbol{\beta}}_{\text{cross-basis}} + \sum_{k} \gamma_k u_{kt}\]

\[\alpha, \boldsymbol{\beta}, \gamma_k \sim \text{(prior distributions)}\]

Every parameter now has its posterior distribution

This enables:

• Full uncertainty quantification of all estimates via posterior distributions

• Uncertainty propagation in derived measures (optimal exposure, attributable burden) without additional approximations

• Hierarchical extensions to add layers of complexity in the model (e.g. Spatial Bayesian DLNMs, Quijal-Zamorano et al. (2024))

{bdlnm}

Published on CRAN since March 2026:

        

Overview of the package

Temperature-mortality study

• We assess the immediate and delayed non-linear effect of daily mean temperature on mortality (age 75+) in London, 2000–2011

• We will use the {bdlnm} built-in london dataset: one observation per day, containing mortality counts, mean temperature, and calendar variables:

(Vicedo-Cabrera et al. 2019)

Model specification

\[Y_t \sim \text{Poisson}(\mu_t), \quad \quad \log(\mu_t) = \alpha + cb(x_t, \ldots, x_{t-L}) \cdot \beta + \sum_{k} \gamma_k u_{kt}\]

Building the cross-basis

• We have to use first the {dlnm} package to build the cross-basis:

library(dlnm)
argvar <- list(
  fun = "ns",
  knots = quantile(london$tmean, c(10, 75, 90) / 100, na.rm = TRUE),
  Bound = range(london$tmean, na.rm = TRUE)
)
arglag <- list(
  fun = "ns",
  knots = logknots(21, nk = 3)
)
cb <- crossbasis(london$tmean, lag = 21, argvar, arglag)

• And build the seasonality and temperature prediction grid:

library(splines)
seas <- ns(london$date, df = round(8 * length(london$date) / 365.25))
temp <- round(seq(min(london$tmean), max(london$tmean), by = 0.1), 1)

Fitting the Bayesian DLNM

mod <- bdlnm(
  mort_75plus ~ cb + factor(dow) + seas,
  data = london,
  family = "poisson",
  sample.arg = list(n = 1000, seed = 5243)
)

bdlnm() handles two things internally:

1. Fits the model using R-INLA
2. Returns posterior samples and summaries for all parameters

List of 4 (class "bdlnm")
 $ model               : list [49]      — fitted R-INLA model object
 $ basis               : list [1]       — cross-basis
 $ coefficients        : num [123, 1000] — posterior samples
 $ coefficients.summary: num [123, 6]    — posterior summary

Fitting the Bayesian DLNM

Now we have a matrix with the estimated posterior samples of the model coefficients, together with their summary across samples:

Posterior samples mod$coefficients

                 sample1      sample2      sample3      sample4      sample5
(Intercept)  5.177126982  5.078820055  5.087683979  5.240468411  5.259946480
cbv1.l1     -0.139254820 -0.154092463 -0.112435103 -0.110085718 -0.153345807
cbv1.l2     -0.007785162  0.008057831  0.001246025 -0.003957057  0.006234502
cbv1.l3     -0.057349056 -0.056023923 -0.045405679 -0.044017954 -0.055955602
cbv1.l4      0.048904544  0.057243423  0.025429828  0.028536089  0.049960368

Posterior summary mod$coefficients.summary

                    mean         sd  0.025quant     0.5quant  0.975quant
(Intercept)  5.204972340 0.10244169  5.00329053  5.205428027  5.39781786
cbv1.l1     -0.134951449 0.02148864 -0.17689662 -0.134049849 -0.09601856
cbv1.l2     -0.006084204 0.01031606 -0.02556688 -0.006093614  0.01455691
cbv1.l3     -0.049708682 0.01141099 -0.07279473 -0.049392945 -0.02789958
cbv1.l4      0.049567985 0.01806032  0.01484373  0.050353168  0.08412770

* Visualizing only the first 5 coefficients and the first 5 posterior samples

Minimum mortality temperature (MMT)

We will use the optimal_exposure() function to calculate it:

mmt <- optimal_exposure(mod, exp_at = temp)

List of 2 (class "optimal_exposure")
 $ est    : num [1:1000] — posterior samples of the MMT
 $ summary: num [1:6]    — posterior summary (mean, sd, quantiles)

MMT posterior distribution

Predicting exposure-lag-response effects

From the fitted model, bcrosspred() computes the full bivariate exposure-lag-response association. Relative risks will be centered at the MMT:

cen <- mmt$summary[["0.5quant"]]
cpred <- bcrosspred(mod, exp_at = temp, cen = cen)

List (class "bcrosspred")
 $ coefficients        : num [20, 1000]      — posterior samples of cross-basis coefficients
 $ coefficients.summary: num [20, 6]         — posterior summary of cross-basis coefficients
 $ matRRfit            : num [326, 22, 1000] — posterior samples of RR (exposure x lag x samples)
 $ matRRfit.summary    : num [326, 22, 6]    — posterior summary of RR (exposure x lag x summaries)
 $ allRRfit            : num [326, 1000]     — posterior samples of overall RR (exposure x samples)
 $ allRRfit.summary    : num [326, 6]        — posterior summary of overall RR (exposure x summaries)

Predicting exposure-lag-response effects

Now we have an array of estimated posterior RRs for each temperature–lag combination, together with their summary across samples:

, , sample1

          lag0     lag1     lag2     lag3
-3.4 0.8489699 1.078465 1.144323 1.100592
-3.3 0.8496304 1.077665 1.143102 1.099767
-3.2 0.8502913 1.076867 1.141883 1.098942
-3.1 0.8509526 1.076069 1.140665 1.098119

, , mean

          lag0     lag1     lag2     lag3
-3.4 0.8526728 1.071554 1.134939 1.094427
-3.3 0.8533699 1.070944 1.133932 1.093743
-3.2 0.8540676 1.070335 1.132926 1.093060
-3.1 0.8547656 1.069726 1.131921 1.092377

plot.bcrosspred(): visualise predictions

plot(cpred, ptype = "3d")

plot.bcrosspred(): visualise predictions

plot(cpred, ptype = "3d")
plot(cpred, ptype = "slices", exp_at = round(quantile(london$tmean, 0.99), 1))

plot.bcrosspred(): visualise predictions

plot(cpred, ptype = "3d")
plot(cpred, ptype = "slices", lag_at = 0)

plot.bcrosspred(): visualise predictions

plot(cpred, ptype = "3d")
plot(cpred, ptype = "overall")

Attributable burden

• The attributable() function estimates the mortality burden attributable to non-optimal temperature exposures:

  • Attributable fraction (AF): proportion of all mortality events attributable to temperature
  • Attributable number (AN): absolute number of deaths attributable to temperature

• If we have estimated \(\beta_x\) associated with a certain exposure (centered to an optimal reference one) we can calculate the attribution to non-optimal temperature exposure as:

\[\text{AF}_{x,t} = 1 - \exp(-\beta_x) \qquad \text{AN}_{x,t} = n_t \cdot \text{AF}_{x,t}\]

• In a DLNM setting, the overall cumulative effect \(\beta_{x,\,\text{overall}}\) can be used as \(\beta_x\) (forward perspective)

Computing attributable burden

attributable() returns daily and aggregated AF and AN using either the forward or backward perspective:

attr <- attributable(mod, london, name_date = "date", name_exposure = "tmean", name_cases = "mort_75plus", cen = cen, dir = "forw")

List of 8
 $ af             : num [4383, 1000] — posterior samples of daily AF
 $ an             : num [4383, 1000] — posterior samples of daily AN
 $ aftotal        : num [1:1000]     — posterior samples of total AF
 $ antotal        : num [1:1000]     — posterior samples of total AN
 $ af.summary     : num [4383, 6]    — posterior summary of daily AF
 $ an.summary     : num [4383, 6]    — posterior summary of daily AN
 $ aftotal.summary: num [1:6]        — posterior summary of total AF
 $ antotal.summary: num [1:6]        — posterior summary of total AN

Daily attributable fraction

Time series of daily attributable fractions (AFs), together with the temperature exposition in that given day:

Total attributable burden

Total mortality burden attributable to non-optimal temperatures in the London 75+ population across 2000–2011:

Conclusions

• The {bdlnm} package provides a complete and accessible implementation of Bayesian Distributed Lag Non-Linear Models in R

• Combining the flexibility of DLNMs with full Bayesian inference enables researchers to fit more realistic and complex models

• The posterior distribution of all estimates allows direct uncertainty quantification of derived measures (e.g., MMT, attributable fractions, attributable numbers) without additional approximations

• Posterior summaries are fully compatible with frequentist workflows, making adoption straightforward for practitioners familiar with {dlnm}

• {bdlnm} constitutes a valuable tool for environmental risk research and, more generally, for any time series analysis with delayed non-linear effects

Present and future work

Coming soon:

• Multi-location analysis: model different exposure-lag-response curves across cities or regions within a single model (e.g., SB-DLNM)

• Parallelization: faster inference via {futurize} for large datasets

• Compatibility beyond INLA: use MCMC algorithm (e.g. NIMBLE, Stan) for Bayesian model estimation

Further resources

📦 Full tutorials and documentation are available on the package website
pasahe.github.io/bdlnm

🐛 Bug reports and contributions are welcome via the GitHub repository
github.com/pasahe/bdlnm

Spatial Bayesian DLNMs (SB-DLNMs)

SB-DLNMs (Quijal-Zamorano et al. 2024) adds a spatial structure to the cross-basis coefficients for reliable small-area estimation:

\[Y_{it} \sim \text{Poisson}(\mu_{it}) \qquad \log(\mu_{it}) = \alpha_i + \underbrace{cb(x_{it}, \ldots, x_{it-L}) \cdot \beta_i}_{\text{spatial cross-basis}} + \sum_{k} \gamma_k u_{kit}\]

\[\beta_i = \mu_\beta + \phi_i + \epsilon_i, \qquad \phi_i \mid \phi_{-i} \sim \mathcal{N}\!\left(\bar{\phi}_{\delta_i},\, \frac{\sigma^2_\phi}{|\delta_i|}\right)\]

Exposure-lag-response surface

The full bivariate association can be represented as a 3-D surface over temperature and lag, giving the RR at each specific temperature-lag combination:

Overall cumulative effect