Skip to contents

Modeling 1: DLNM — Distributed Lag Non-Linear Model

climasus4r Team

17/10/2026

Source: vignettes/mod-01-dlnm.Rmd
mod-01-dlnm.Rmd

Objective

By the end of this tutorial you will be able to fit a DLNM (Distributed Lag Non-linear Model) to quantify the relationship between temperature and daily mortality, interpret the cumulative exposure-response curve, the lag-response curve, and the two-dimensional exposure × lag surface, and become familiar with all the visualizations available in sus_mod_plot_dlnm().

The DLNM is the cornerstone of climasus4r’s biostatistics hat. It connects the daily health series (produced in Tutorials 1–4) to the aligned climate exposure (Tutorials 7–8) and produces relative-risk estimates ready for attributable-burden quantification (Modeling Tutorial 2), multi-city meta-analysis (Modeling Tutorial 5), and vulnerability analysis (Modeling Tutorial 9).

Scientific context

Modern environmental epidemiology recognizes that the association between temperature and mortality is non-linear and lagged in time. Extreme heat kills within the first 1–2 days (acute effect); cold kills over 1 to 3 weeks (delayed, cumulative effect). Modeling only same-day temperature systematically underestimates the risk of cold and can mask the winter mortality peak in subtropical climates.

Gasparrini, Armstrong & Kenward (2010) proposed the DLNM to solve this problem with a cross-basis: the outer product of two spline bases — one on the exposure dimension (captures the U/J shape: risk rises with both cold and heat) and one on the lag dimension (captures how the effect is distributed over the following days). The R implementation is the dlnm package (Gasparrini, 2011), on top of which climasus4r builds sus_mod_dlnm().

The fitting procedure follows the best practices for time-series regression in environmental epidemiology systematized by Bhaskaran et al. (2013): a log-linear quasi-Poisson model, control of long-term trend and seasonality via a natural spline of time, and relative risk centered at the median exposure as the reference value (Armstrong, 2006).

The detailed formulas for the crossbasis and the cumulative RR are presented in the Formulas section below.

Driving case. We investigate temperature and health outcomes in older adults (60+) in the São Paulo Metropolitan Region, 2014–2019 (respiratory diseases J00–J99). We join the daily death series (dados/caso_serie.rds, Tutorial 4) with the temperature from INMET station A701 (dados/caso_clima_estacao.rds, Tutorial 7), build the distributed-lag structure, and fit the quasi-Poisson DLNM to estimate the exposure-response curve and the peak lag.

Formulas

The DLNM models the expected number of events YtY_t on date tt as:

log⁡{E(Yt)}=α+∑ℓ=0Ls(xt−ℓ,ℓ;𝛃)+ns(t)+γDOWt \log\{E(Y_t)\} = \alpha + \sum_{\ell=0}^{L} s(x_{t-\ell},\, \ell;\,\boldsymbol{\beta}) + \mathrm{ns}(t) + \gamma\,\mathrm{DOW}_t

where s(xt−ℓ,ℓ;𝛃)s(x_{t-\ell}, \ell;\,\boldsymbol{\beta}) is the cross-basis function — the outer product of a spline basis over the exposure xt−ℓx_{t-\ell} (temperature at lag ℓ\ell) and a spline basis over the lag ℓ∈{0,…,L}\ell \in \{0, \ldots, L\}, with coefficient vector 𝛃\boldsymbol{\beta}; ns(t)\mathrm{ns}(t) is a natural spline of time controlling for seasonality and long-term trend; and γDOWt\gamma\,\mathrm{DOW}_t controls for the day-of-week effect.

The cumulative effect of exposure xx relative to the reference value x0x_0 (lag-0 median) is:

RRcum(x)=exp⁡{∑ℓ=0L[s(x,ℓ;𝛃̂)−s(x0,ℓ;𝛃̂)]} RR_{\text{cum}}(x) = \exp\!\left\{ \sum_{\ell=0}^{L} \bigl[s(x,\, \ell;\,\boldsymbol{\hat\beta}) - s(x_0,\, \ell;\,\boldsymbol{\hat\beta})\bigr] \right\}

Key parameters implemented by sus_mod_dlnm():

Parameter Notation Controlled by
Exposure basis s(⋅;𝛃)s(\cdot; \boldsymbol{\beta}) over xx argvar
Lag basis s(⋅;𝛃)s(\cdot; \boldsymbol{\beta}) over ℓ\ell arglag
Lag window LL lag_max
Seasonality df df of ns(t)\mathrm{ns}(t) dof_per_year
Reference value x0x_0 ref_value

Data used

Input Origin 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 22 distributed-lag columns (tair_dry_bulb_c_lag0, …, tair_dry_bulb_c_lag21), a structure that in a real workflow is produced by sus_climate_aggregate(temporal_strategy = "distributed_lag", lag_days = 21). The script guards each step with file.exists() and, if the inputs are missing, generates a minimal, realistic synthetic series (with a V-shaped temperature-mortality signal) so the code still runs without the full pipeline.

Complete pipeline

The blocks below reproduce exactly what scripts/modelagem_01_dlnm.R executes.

library(climasus4r)

Block 0 — Load the chained inputs 

# Caminhos das entradas encadeadas
path_saude <- file.path("vignettes-pt", "dados", "caso_serie.rds")
path_clima <- file.path("vignettes-pt", "dados", "caso_clima_estacao.rds")

# Série diária de saúde
saude_diaria <- NULL
if (file.exists(path_saude)) {
  caso_serie   <- readRDS(path_saude)
  saude_df     <- dplyr::as_tibble(caso_serie)
  saude_diaria <- saude_df |>
    dplyr::summarise(n_obitos = sum(n_obitos, na.rm = TRUE), .by = date) |>
    dplyr::arrange(date)
}

# Série diária de temperatura
clima_diario <- NULL
if (file.exists(path_clima)) {
  caso_clima   <- readRDS(path_clima)
  clima_diario <- caso_clima$clima_diario |>
    dplyr::transmute(date = as.Date(date), tair_dry_bulb_c) |>
    dplyr::arrange(date)
}

# Fallback sintético realista (série em V, 6 anos)
if (is.null(clima_diario) || is.null(saude_diaria)) {
  # ... veja scripts/modelagem_01_dlnm.R para a implementação completa
}

# Juntar por data
serie <- dplyr::inner_join(saude_diaria, clima_diario, by = "date") |>
  dplyr::arrange(date) |>
  tidyr::drop_na(n_obitos, tair_dry_bulb_c)

Block 1 — Build the distributed-lag structure

In a complete climasus4r workflow, the _lagN columns are produced by sus_climate_aggregate(temporal_strategy = "distributed_lag", lag_days = 21). Here, with the series already aligned by date, we generate the lags directly:

VAR_TEMP <- "tair_dry_bulb_c"
LAG_MAX  <- 21L

df_lags <- serie
for (k in 0:LAG_MAX) {
  df_lags[[paste0(VAR_TEMP, "_lag", k)]] <- dplyr::lag(serie$tair_dry_bulb_c, n = k)
}

# Seleciona apenas: date, desfecho, colunas de lag
df_lags <- df_lags |>
  dplyr::select(
    date, n_obitos,
    dplyr::all_of(paste0(VAR_TEMP, "_lag", 0:LAG_MAX))
  )

# Promove a climasus_df com stage = "climate", type = "distributed_lag"
df_dl <- dplyr::as_tibble(df_lags)
attr(df_dl, "sus_meta") <- list(
  system  = "SIM",
  stage   = "climate",
  type    = "distributed_lag",
  backend = "tibble",
  history = "[derivado] serie diaria -> lags distribuidos tair_dry_bulb_c"
)
class(df_dl) <- c("climasus_df", class(df_dl))

Block 2 — sus_mod_dlnm(): fitting the quasi-Poisson DLNM 

fit <- sus_mod_dlnm(
  df           = df_dl,
  outcome_col  = "n_obitos",
  climate_col  = "tair_dry_bulb_c",
  lag_max      = 21L,
  covariates   = NULL,                       # confundidores adicionais (opcional)
  argvar       = list(fun = "ns", df = 4),   # spline natural, 4 gl exposição
  arglag       = list(fun = "ns", df = 3),   # spline natural, 3 gl lag
  family       = "quasipoisson",             # robusto à superdispersão
  ns_df        = NULL,                       # calculado automaticamente via dof_per_year
  dof_per_year = 4L,                         # gl/ano para sazonalidade (use 6 para Brasil subtropical)
  ref_value    = NULL,                       # mediana de lag-0 = valor de referência
  pred_at      = c(0.25, 0.50, 0.75, 0.90, 0.95, 0.99),
  alpha        = 0.05,
  lang         = "pt",
  verbose      = TRUE
)

# Inspecionar o resultado
print(fit)
fit$exposure_response   # RR cumulativo nos percentis pred_at
fit$lag_response        # RR por lag no percentil 75
fit$diagnostics         # razão de dispersão, Ljung-Box, AIC

Block 3 — sus_mod_plot_dlnm(): all visualizations 

# overall: curva exposição-resposta cumulativa (painel principal + histograma)
sus_mod_plot_dlnm(fit, type = "overall", lang = "pt")

# lag: curva lag-resposta no percentil 75 (RR específico e cumulativo)
sus_mod_plot_dlnm(fit, type = "lag", lang = "pt")

# surface: superfície 3-D bidimensional exposição × lag
# (interactive = FALSE usa contorno estático como fallback)
sus_mod_plot_dlnm(fit, type = "surface", interactive = FALSE, lang = "pt")

# contour: mapa de contorno 2-D log(RR) para publicação
sus_mod_plot_dlnm(fit, type = "contour", lang = "pt")

# slice: curvas de resposta por lag específico (valores > lag_max são descartados)
sus_mod_plot_dlnm(fit, type = "slice",
                  lags_at = c(0L, 3L, 7L, 14L, 21L), lang = "pt")

# distribution: histograma/densidade da exposição com percentis-chave
sus_mod_plot_dlnm(fit, type = "distribution", lang = "pt")

# series: série temporal desfecho (barras) + exposição lag-0 (eixo dual)
sus_mod_plot_dlnm(fit, type = "series", lang = "pt")

# output_type = "all": retorna lista nomeada $plot, $table, $data
out <- sus_mod_plot_dlnm(fit, type = "overall",
                         output_type = "all", lang = "pt")
out$table   # lista de 3 gt_tbl: especificação, curva RR, diagnósticos
out$data    # lista: $exposure_response, $lag_response, $data_daily, $diagnostics

Block 4 — Save the fit for reuse (M2 and M8) 

# O objeto climasus_dlnm é salvo para alimentar
# Tutorial M2 (fração atribuível) e M8 (meta-análise multi-cidade).
# save_tbl() (definida em scripts/_setup_tutoriais.R) persiste em
# vignettes-pt/dados/mod_dlnm_fit.rds; o script também copia para
# dados/mod_dlnm_fit.rds (raiz do pacote) para encadeamento com M2 e M8.
save_tbl(fit, "mod_dlnm_fit")
saveRDS(fit, file.path("dados", "mod_dlnm_fit.rds"))

Function characteristics

sus_mod_dlnm()

Fits a DLNM from a climasus_df at stage "climate" / type "distributed_lag", produced by sus_climate_aggregate(temporal_strategy = "distributed_lag"). Builds the crossbasis, fits a GLM, and precomputes the crosspred on a grid of 100 points between the 1st and 99th percentile of exposure.

Argument Default Description
df — climasus_df with date, outcome, and {var}_lagN columns
outcome_col "n_obitos" Name of the daily outcome count column
climate_col NULL Exposure name base (without _lagN); auto-detected by priority if NULL
lag_max NULL Maximum lag; auto-detected from the largest N among _lagN columns
covariates NULL Additional confounders to include linearly in the GLM
argvar list(fun="ns", df=4) Basis for the exposure dimension of the crossbasis
arglag list(fun="ns", df=3) Basis for the lag dimension of the crossbasis
family "quasipoisson" GLM family; "negbin" is redirected to "quasipoisson"
ns_df NULL df for the temporal spline; NULL uses dof_per_year × n_years
dof_per_year 4L df per year for seasonality (6 recommended for subtropical Brazil)
ref_value NULL Reference value (RR = 1); NULL uses the lag-0 median
pred_at c(.25,.5,.75,.9,.95,.99) Quantiles to report in the $exposure_response table
alpha 0.05 CI level (default: 95% CI)
lang "pt" Message language: "pt", "en", "es"

Return — climasus_dlnm object (S3 list):

Component Content
$model Fitted glm object (compatible with dlnm::crosspred())
$crossbasis dlnm::crossbasis object (exposure × lag basis)
$pred dlnm::crosspred object — allRRfit/low/high (overall effect) and matRRfit/low/high (lag-specific)
$exposure_response Tibble: pct, exposure, rr, lo, hi at the pred_at percentiles
$lag_response Tibble: lag, rr, lo, hi, rr_cum — RR per lag at P75
$models 1-row summary tibble: variable, n, family, lag_max, rr, lo, hi, lag_peak, disp_ratio, aic_poisson
$data_daily Tibble aggregated by date (outcome summed; climate averaged) used in the fit
$diagnostics List: disp_ratio, disp_category, autocorr_pval, has_autocorr, aic_poisson, deviance, df_residual
$meta List with all fit parameters plus the sus_meta history

Basis selection (argvar / arglag):

argvar When to use
list(fun="ns", df=3) Smooth non-linearity; conservative
list(fun="ns", df=4) U/J-shaped curve (heat + cold, default)
list(fun="ns", df=5) Complex dose-response; large N
list(fun="lin") Linear hypothesis test
arglag When to use
list(fun="ns", df=3) Smooth decay; temperature-related mortality (default)
list(fun="ns", df=4) Bimodal lag; cold effect
list(fun="poly", df=2) Conservative; small samples

Seasonality control (dof_per_year):

  • 4 df/year: 1–2 year series; annual seasonality
  • 6 df/year: subtropical Brazil; more pronounced seasonality
  • 8 df/year: arboviruses; sub-seasonal variation

Available S3 methods: print(), summary(), coef(), tidy() (returns a one-row tibble compatible with dplyr::bind_rows()).

sus_mod_plot_dlnm()

Produces scientific visualizations and statistical tables from a climasus_dlnm object. Implements all seven standard visualization types used in DLNM epidemiology.

Argument Default Description
fit — climasus_dlnm object from sus_mod_dlnm()
type "overall" Visualization type (see table below)
output_type "plot" "plot" → ggplot2/plotly; "table" → gt; "all" → list $plot, $table, $data
exposure_at NULL Focal exposure for the "lag" plot (default: P75)
lags_at c(0,3,7,14,21) Lags shown in the "slice" type (clipped at lag_max)
pred_at c(.05,.10,.25,.50,.75,.90,.95,.99) Quantiles in the statistical table
interactive FALSE TRUE → plotly::ggplotly() / 3-D plotly surface
color_palette "npg" ggsci palette (others: "lancet", "jama", "nejm")
lang "pt" Label language: "pt", "en", "es"

Visualization types (type =)

type What it shows Reference
"overall" Cumulative exposure-response curve (RR + 95% CI) with the exposure histogram in the bottom panel Gasparrini et al. (2010)
"lag" Lag-specific and cumulative RR at a given exposure level (default P75) Gasparrini et al. (2010)
"surface" Two-dimensional 3-D exposure × lag surface (interactive = TRUE → 3-D plotly; static → contour) Gasparrini (2011)
"contour" 2-D contour map (log(RR), divergent scale) for publication Gasparrini et al. (2014)
"slice" Overlaid response curves for specific lags (dose-response profiles by lag time) Armstrong (2006)
"distribution" Histogram + kernel density of exposure with markers at key percentiles Bhaskaran et al. (2013)
"series" Daily time series of the outcome (bars) and exposure (line, dual axis) Bhaskaran et al. (2013)

Return (output_type):

  • "plot" → ggplot or plotly object
  • "table" → list of gt_tbl (or tibble if gt is unavailable) with 3 sections: (1) model specification, (2) RR curve at percentiles,
    1. model diagnostics
  • "all" → named list $plot, $table, $data (containing $exposure_response, $lag_response, $data_daily, $diagnostics)

Illustrated results

Time series: outcome and exposure

Figure 1. Daily time series of deaths among older adults (green bars) and mean temperature (orange line, scaled secondary axis), São Paulo Metropolitan Region 2014–2019. The seasonal pattern is visible: mortality peaks coincide with the subtropical winter (June–August), when minimum daily temperatures drop below the thermal comfort threshold. Exposure to extreme heat (austral summer, December–February) produces a more acute effect.

Temperature distribution

Figure 2. Distribution of daily temperature (lag 0) with key percentiles marked (P10, P25, P50 — reference, P75, P90, P99). The median is the reference value (RR = 1) of the DLNM. The P25–P75 range represents “usual” thermal conditions; values above P90 and below P10 define the extreme heat and cold thresholds.

Cumulative exposure-response curve (type = "overall")

Figure 3. Cumulative exposure-response curve (RR and 95% CI via the delta method). The V-shaped pattern is characteristic of the temperature-mortality association in subtropical Brazil: elevated risk under moderate cold (delayed, cumulative effect) and, to a lesser degree, under extreme heat (acute effect). The horizontal dashed line marks RR = 1; the vertical dotted line marks the median exposure (reference value). The colored points indicate the quantiles reported in Table 1.

Lag-response curve (type = "lag")

Figure 4. Lag-response curve for the 75th percentile of exposure. The solid line shows the lag-specific RR (effect on each individual day); the dashed line shows the cumulative RR accumulated up to that lag. The peak lag — the day at which the lag-specific effect is maximal — is reported in the model diagnostics.

Contour map (type = "contour")

Figure 5. Contour map of log(RR) as a function of temperature (y-axis) and lag day (x-axis). Warm colors (orange/red) indicate excess risk; cool colors (green) indicate reduced risk relative to the median. The map simultaneously reveals the shape (non-linearity with respect to temperature) and the temporal distribution (how the effect persists across lags). White isolines mark the ±0.1 and ±0.2 log(RR) levels.

Curves by specific lag (type = "slice")

Figure 6. Exposure-response curves for selected lags (0, 3, 7, 14, and 21 days). Each curve shows how the dose-response profile changes over time. The extreme-heat effect is more pronounced at lags 0–3 (acute response); the moderate-cold effect persists up to lags 14–21 (delayed, distributed response).

Model statistical summary

Tabela 1. Especificação e diagnósticos do DLNM ajustado.
Parametro Valor
Exposicao climatica tair_dry_bulb_c
Desfecho n_obitos
Familia GLM quasipoisson
Lag maximo (dias) 21
Base exposicao (argvar) ns(df=4)
Base lag (arglag) ns(df=3)
gl/ano (sazonalidade) 4
N observacoes (dias) 2170
Valor de referencia (degC) 18.8
RR cumulativo no P75 0.7355
IC 95% RR (P75) [0.4858; 1.1134]
Lag de efeito pico (dias) 0
Razao de dispersao (phi) 1.132
Categoria dispersao adequate
AIC (Poisson re-fit) 8051.1
Ljung-Box p-valor 0.0138
Autocorrelacao residual Sim (p < 0.05)
Tabela 2. Curva exposição-resposta: RR cumulativo nos percentis da temperatura.
pct exposure rr lo hi
0.25 15.980 0.7695253 0.5363601 1.104052
0.50 18.805 1.0003781 0.9938598 1.006939
0.75 21.950 0.7354524 0.4857910 1.113422
0.90 23.711 0.8367105 0.4584446 1.527086
0.95 24.376 0.9651056 0.4831790 1.927710
0.99 25.673 1.3823675 0.5457206 3.501682

Table 2. Cumulative Relative Risk (95% CI) for each percentile of the temperature distribution, relative to the median (RR = 1). Values above 1.00 indicate excess risk; values below 1.00 indicate relative protection. Percentiles above P90 and below P10 indicate extreme exposures.

Interpretation and best practices

Dispersion ratio (ϕ̂\hat\phi). The quasi-Poisson is recommended when ϕ̂>1,5\hat\phi > 1{,}5 (Armstrong, 2006). The dispersion ratio is reported in fit$diagnostics$disp_ratio; category "adequate" indicates ϕ̂<1,5\hat\phi < 1{,}5.

Residual autocorrelation. The Ljung-Box test (10 lags of the deviance residuals) detects remaining autocorrelation after seasonality control. If positive (has_autocorr = TRUE), increase dof_per_year or add covariates = "dow" (day of week) (Bhaskaran et al., 2013).

Degrees of freedom of the temporal spline. The general rule is dof_per_year = 4 for 1–2 years of data and dof_per_year = 6 for subtropical Brazil (more pronounced seasonality). For arboviruses with sub-seasonal peaks, dof_per_year = 8. Visualize the seasonal residual with fit$data_daily to check the fit.

Reference value. The default is the median exposure at lag-0. In multi-city comparative studies, use a fixed common value (ref_value = 20.0) to ensure comparability across exposure-response curves.

Chaining. The mod_dlnm_fit.rds object (a climasus_dlnm component) is the input for the tutorials:

  • M2 (Attributable Fraction): sus_mod_af(fit = fit, ...)
  • M5 (Meta-analysis): sus_mod_pool(fits = list(rmsp = fit, ...), ...)
  • M8 (Space-Time Bayes): model with an offset based on the estimated RR

References

  • Armstrong, B. (2006). Models for the relationship between ambient temperature and daily deaths. Epidemiology, 17(6), 624–631. https://doi.org/10.1097/01.ede.0000239732.50999.58

  • 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., 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. (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. (2014). 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

  • Saldanha, R. F., Bastos, R. R., & Barcellos, C. (2019). Microdatasus: pacote para download e pré-processamento de microdados do DATASUS. Cadernos de Saúde Pública, 35(9), e00032419.

How to reproduce the figures. Run once, from the package root: Rscript scripts/modelagem_01_dlnm.R The script saves the figures to vignettes-pt/figuras/ (prefix moddlnm_), the tables to vignettes-pt/dados/, and the fitted object mod_dlnm_fit.rds (input for M2 and M8).