Skip to contents

Tutorial 9: Epidemiological Modeling (DLNM and Attributable Burden)

climasus4r Team

17/10/2026

Source: vignettes/09-modelagem.Rmd
09-modelagem.Rmd

Objective

By the end of this tutorial you will be able to turn a daily health series aligned with climate into epidemiological inference: estimate the non-linear and lagged relationship between temperature and mortality with a DLNM (sus_mod_dlnm()), quantify the attributable burden of heat and cold (sus_mod_af()), and get to know the other tools in climasus4r’s biostatistical toolkit — interrupted time series (sus_mod_its()), case-crossover (sus_mod_casecrossover()), and multi-city meta-analysis (sus_mod_pool()).

Modeling is the last stage of the pipeline. It consumes the output of the previous stages (the aggregated health series from Tutorial 4 and the aligned climate exposure from Tutorials 7 and 8) and produces relative risk estimates, attributable fractions, and publication-ready visualizations.

Scientific context

Modern environmental epidemiology treats the association between an environmental exposure (temperature, pollutants) and a daily count outcome (deaths, hospitalizations) as a problem that is both non-linear and lagged in time. Extreme heat kills on the same day or within 1–2 days (acute effect); cold kills over one to two weeks (delayed, prolonged effect). Modeling only the day’s temperature underestimates the effect of cold.

The Distributed Lag Non-linear Model (DLNM) of Gasparrini, Armstrong & Kenward (2010) solves this with a cross-basis: the outer product of two spline bases — one on the exposure dimension (captures the U/J-shaped curve: risk rising with both cold and heat) and one on the lag dimension (captures how the effect is distributed over the days). The R implementation is the dlnm package (Gasparrini, 2011), on top of which climasus4r builds sus_mod_dlnm().

The fit follows the good practices for time-series regression in environmental epidemiology systematized by Bhaskaran et al. (2013): a log-linear quasi-Poisson model (robust to overdispersion), control for long-term trend and seasonality via a natural spline of time, and adjustment for day of week. The attributable fraction decomposes the total burden of events into heat and cold components (Gasparrini & Armstrong, 2013).

Guiding case. We investigate temperature and health outcomes in older adults (60+) in the São Paulo Metropolitan Region, 2014–2019 (respiratory J00–J99 and cardiovascular I00–I99). Here we join the daily death series (caso_serie.rds, Tutorial 4) with the temperature from INMET station A701 (São Paulo – Mirante de Santana), already aligned (caso_clima_estacao.rds, Tutorial 7), build the distributed-lag structure, and fit the DLNM to estimate the exposure-response curve and the burden attributable to heat and cold.

Data used

Input Source Content Role
dados/caso_serie.rds Tutorial 4 Daily death series (older adults 60+) Outcome n_obitos
dados/caso_clima_estacao.rds Tutorial 7 Daily temperature from station A701 Exposure tair_dry_bulb_c

The two series are joined by date. The exposure is then expanded into distributed-lag columns (tair_dry_bulb_c_lag0, …, tair_dry_bulb_c_lag21), a structure that in a full workflow is produced by sus_climate_aggregate(temporal_strategy = "distributed_lag", lag_days = 21). The script guards this stage with requireNamespace()/tryCatch() and, if the inputs are missing, generates a minimal, realistic synthetic series (with a V-shaped temperature-mortality signal) purely for illustration.

Full pipeline

The blocks below reproduce what scripts/tutorial_09_modelagem.R executes.

library(climasus4r)

Building the DLNM input (distributed lags)

sus_mod_dlnm() expects a climasus_df at stage "climate" and type "distributed_lag", with the date column, the outcome, and the {variable}_lag0 … _lagL columns. In a real workflow this comes from sus_climate_aggregate():

df_dl <- sus_climate_aggregate(
  health_data       = caso_serie,          # health series (from Tutorial 4)
  climate_data      = caso_clima_estacao,  # aligned temperature (from Tutorial 7)
  climate_var       = "tair_dry_bulb_c",
  temporal_strategy = "distributed_lag",
  lag_days          = 21,
  lang              = "pt"
)

Fitting the DLNM (sus_mod_dlnm)

Cross-basis with exposure on a natural spline with 4 degrees of freedom (captures the U/J heat+cold shape) and lag on a natural spline with 3 degrees (smooth decay), quasi-Poisson family, and automatic seasonal control (dof_per_year = 6, recommended for subtropical Brazil):

fit <- sus_mod_dlnm(
  df           = df_dl,
  outcome_col  = "n_obitos",
  climate_col  = "tair_dry_bulb_c",
  lag_max      = 21,
  argvar       = list(fun = "ns", df = 4),  # exposure dimension
  arglag       = list(fun = "ns", df = 3),  # lag dimension
  family       = "quasipoisson",
  dof_per_year = 6L,                        # df/year for seasonality
  lang         = "pt"
)

print(fit)
fit$exposure_response   # cumulative RR on the percentile grid
fit$lag_response        # RR by lag at the 75th percentile
fit$diagnostics$disp_ratio

Attributable fraction from heat and cold (sus_mod_af)

From the fitted DLNM, we decompose the total burden of deaths into heat (above the comfort threshold) and cold (below it) components, with confidence intervals from Monte Carlo simulation:

af <- sus_mod_af(
  fit  = fit,
  by   = "month",   # optional monthly breakdown
  nsim = 200L,      # use 5000L for publication-quality results
  lang = "pt"
)

print(af)
af$total        # total AF/AN, heat and cold
af$by_period    # monthly breakdown

Visualizing (sus_mod_plot_dlnm)

The package offers ready-made scientific visualizations — exposure-response curve, lag-response, 3-D surface, contour, lag slices, distribution, and time series:

sus_mod_plot_dlnm(fit, type = "overall", lang = "pt")  # exposure-response
sus_mod_plot_dlnm(fit, type = "contour", lang = "pt")  # log(RR) surface
sus_mod_plot_dlnm(fit, type = "lag",     lang = "pt")  # lag-response

Function characteristics

sus_mod_dlnm() — Distributed Lag Non-linear Model

Argument Essential? Description
df yes climasus_df at stage = "climate", type = "distributed_lag"
outcome_col yes Daily count column (e.g., "n_obitos")
climate_col optional Exposure base without the _lagN suffix; NULL auto-detects
lag_max optional Maximum lag; NULL detects it from the largest _lagN present
argvar optional Exposure basis (e.g., list(fun = "ns", df = 4))
arglag optional Lag basis (e.g., list(fun = "ns", df = 3))
family optional "quasipoisson" (default), "poisson", "negbin" (→ quasi)
ns_df / dof_per_year optional Degrees of freedom of the time spline (seasonality)
ref_value optional Reference value (RR = 1); NULL uses the median
covariates optional Additional linear confounders (humidity, etc.)

Returns: a climasus_dlnm object with $model (GLM), $crossbasis, $pred (crosspred on a 100-point grid), $exposure_response, $lag_response, $models (one-line summary), $data_daily, and $diagnostics (dispersion, Ljung-Box autocorrelation, AIC).

sus_mod_af() — Attributable Fraction and Number

Argument Description
fit climasus_dlnm object from sus_mod_dlnm()
threshold Threshold separating heat (above) and cold (below); NULL uses ref_value
pred_at Percentiles of the ranges (default c(0.75, 0.90, 0.95, 0.99))
by Temporal breakdown: "month", "year", "season", or NULL
nsim Monte Carlo simulations for the CI (default 1000L)
alpha Significance level (default 0.05 → 95% CI)

Returns: climasus_af with $total (total/heat/cold), $by_quantile, $by_period, $daily, and $custom.

sus_mod_its() — Interrupted Time Series

Segmented Poisson regression for the effect of interventions (policies, extreme events, pandemics) on a count series (Bernal et al., 2017).

Argument Description
data climasus_df (aggregate/climate) or data.frame with date + count
outcome_col / date_col Outcome and date columns
interruption_dates One or more intervention dates (within the range)
harmonics Sinusoidal pairs for seasonality (default 2L)
family "quasipoisson" (default) or "poisson"
counterfactual Projects the pre-intervention trend (default TRUE)

Returns: climasus_its with $effects (level and trend change), $counterfactual, and $segments.

sus_mod_casecrossover() — Time-Stratified Case-Crossover

Each time stratum (month/week) acts as its own control, removing trend and seasonality by design (Maclure, 1991; Armstrong et al., 2014).

Argument Description
data climasus_df or data.frame with date, outcome, and exposure
outcome_col / exposure_col Outcome and exposure (no auto-detection)
stratum "month" (default), "week", or an existing column
lag Single integer lag or vector (moving average, e.g., 0:6)
method "conditional_poisson" (default) or "clogit"
family "quasipoisson" (default) or "poisson"

Returns: climasus_casecrossover with $or_table (OR, CI, p) and $diagnostics.

sus_mod_pool() — Multi-City Meta-Analysis

Combines DLNM estimates from several cities via multivariate meta-analysis (REML via mvmeta), producing a pooled curve, heterogeneity (Q, I²), and per-city BLUPs (Gasparrini et al., 2012).

Argument Description
fits Named list of climasus_dlnm objects (same cross-basis)
method "reml" (default), "ml", or "fixed"
blup Computes per-city BLUP predictions (default TRUE)
pred_at Percentiles of the pooled table

Returns: climasus_pool with $exposure_response, $lag_response, $heterogeneity, $city_table, and $blup_preds.

Statistical formulation

The DLNM models the daily count YtY_t as quasi-Poisson with mean μt\mu_t. The log-linear predictor combines the cross-basis cb(xt,ℓ;β)\mathrm{cb}(x_t, \ell; \beta) — the effect of the exposure history xt,xt−1,…,xt−Lx_{t}, x_{t-1}, \ldots, x_{t-L} — with control for seasonality (natural spline of time) and day of week (DOW):

log⁡𝔼(Yt)=α+cb(xt,ℓ;β)+ns(tempo)+DOW. \log\,\mathbb{E}(Y_t) \;=\; \alpha \;+\; \mathrm{cb}(x_t,\,\ell;\,\beta) \;+\; \mathrm{ns}(\text{tempo}) \;+\; \mathrm{DOW}.

The cross-basis is the outer product of a basis on exposure (argvar) and a basis on lag (arglag). The cumulative effect across all lags, at exposure level xx relative to the reference x0x_0, is

RRcum(x)=exp⁡{∑ℓ=0L[f(x,ℓ)−f(x0,ℓ)]}. RR_{\text{cum}}(x) \;=\; \exp\!\left\{ \sum_{\ell=0}^{L} \big[ f(x,\ell) - f(x_0,\ell) \big] \right\}.

The cumulative attributable fraction associated with the exposure history follows the form (Gasparrini & Armstrong, 2013):

AF=1−exp⁡(−∑ℓ=0Lβℓxt−ℓ),ANt=Yt(1−1RRt), AF \;=\; 1 - \exp\!\left( -\sum_{\ell=0}^{L} \beta_\ell\, x_{t-\ell} \right), \qquad AN_t \;=\; Y_t \left( 1 - \frac{1}{RR_t} \right),

with the daily attributable number ANtAN_t summed over the period and the CI obtained by Monte Carlo simulation of the cross-basis coefficients.

Modeling decisions

  • Overdispersion. Death counts typically have variance greater than the mean (ϕ̂>1\hat\phi > 1). The quasi-Poisson family inflates the standard error by ϕ̂\sqrt{\hat\phi} and is recommended when the dispersion ratio exceeds ~1.5 (Armstrong, 2006). The negative binomial (NB) alternative is incompatible with dlnm::crossbasis(), and climasus4r redirects it to quasi-Poisson.
  • Choice of knots and lags. For exposure, ns with 4 df captures the U/J shape (cold + heat) well; 3 df is more conservative, 5 df requires a large N. For lag, ns with 3 df models the smooth decay typical of the temperature-mortality relationship; a maximum lag of 21 days covers the delayed, prolonged effect of cold (Gasparrini, 2011).
  • Seasonality. The natural spline of time absorbs long-term trend and the annual cycle that confound the association. A value of ~dof_per_year = 4 is recommended for 1–2 years and 6–8 for subtropical Brazil and diseases with strong seasonality (Bhaskaran et al., 2013). Underfitting leaves residual autocorrelation (Ljung-Box test in $diagnostics); overfitting consumes the signal.

Illustrated results

Figure 1. DLNM response surface: log(RR) of deaths in older adults as a function of temperature (y-axis) and lag day (x-axis). Orange tones indicate increased risk (heat, short lags); green tones, reduced risk. The two-dimensional structure reveals when, after each exposure level, the effect manifests.

Figure 2. Cumulative exposure-response curve (RR and 95% CI) summed over all lags, centered on the median temperature (reference, RR = 1). The U/J shape — rising risk with both cold and heat — is the classic signature of the temperature-mortality relationship.

Figure 3. Lag-response curve at the 75th percentile of exposure: lag-specific RR between temperature and death. It allows distinguishing the acute effect (short lags) from the delayed, prolonged effect (long lags).

Figure 4. Attributable fraction of deaths (%) decomposed into total, heat, and cold, with 95% Monte Carlo CI. In subtropical climates the cold fraction tends to exceed the heat fraction, reflecting the prolonged effect of low temperatures on older adults.

Table 1. Summary of the DLNM fit (specification, RR at P75, and diagnostics).
Parametro Valor
Exposicao climatica tair_dry_bulb_c
Desfecho n_obitos
Familia GLM quasipoisson
Lag maximo (dias) 21
N observacoes (dias) 2170
Valor de referencia (°C) 18.8
RR no P75 (vs referencia) 0.922 [0.725; 1.172]
Lag de efeito pico (dias) 21
Razao de dispersao 1.05 (adequate)
AIC (Poisson) 11800.7
Table 2. Attributable fraction (AF) and attributable number (AN) by thermal component.
Componente FA (%) NA (obitos) N obitos
Total 5.78 [-7.25; 17.87] 1624 [-2037; 5018] 28081
Calor 2.90 [-3.93; 8.48] 814 [-1102; 2382] 28081
Frio 2.89 [-3.35; 9.49] 811 [-941; 2665] 28081

References

  • Armstrong, B. (2006). Models for the relationship between ambient temperature and daily mortality. Epidemiology, 17(6), 624–631.
  • Bernal, J. L., Cummins, S., & Gasparrini, A. (2017). Interrupted time series regression for the evaluation of public health interventions: a tutorial. International Journal of Epidemiology, 46(1), 348–355. https://doi.org/10.1093/ije/dyw098
  • Bhaskaran, K., Gasparrini, A., Hajat, S., Smeeth, L., & Armstrong, B. (2013). Time series regression studies in environmental epidemiology. International Journal of Epidemiology, 42(4), 1187–1195. https://doi.org/10.1093/ije/dyt092
  • Gasparrini, A. (2011). Distributed lag linear and non-linear models in R: the package dlnm. Journal of Statistical Software, 43(8), 1–20. https://doi.org/10.18637/jss.v043.i08
  • Gasparrini, A., Armstrong, B., & Kenward, M. G. (2010). Distributed lag non-linear models. Statistics in Medicine, 29(21), 2224–2234. https://doi.org/10.1002/sim.3940
  • Gasparrini, A., & Armstrong, B. (2013). Reducing and meta-analysing estimates from distributed lag non-linear models. BMC Medical Research Methodology, 14, 70. https://doi.org/10.1186/1471-2288-14-70
  • Maclure, M. (1991). The case-crossover design: a method for studying transient effects on the risk of acute events. American Journal of Epidemiology, 133(2), 144–153.