Distributed Lag Non-linear Model (DLNM) for Climate-Health Analyses

Fits a Distributed Lag Non-linear Model (DLNM) to quantify the association between a climate exposure and a daily health outcome count. The function accepts a climasus_df at stage "climate" / type "distributed_lag" produced by sus_climate_aggregate(temporal_strategy = "distributed_lag"), constructs the bidimensional crossbasis, fits a GLM, and returns a climasus_dlnm object with the fitted model, pre-computed crosspred, exposure-response and lag-response tables, and diagnostic statistics.

Usage

sus_mod_dlnm(
  df,
  outcome_col = "n_obitos",
  climate_col = NULL,
  lag_max = NULL,
  covariates = NULL,
  argvar = list(fun = "ns", df = 4),
  arglag = list(fun = "ns", df = 3),
  family = "quasipoisson",
  ns_df = NULL,
  dof_per_year = 4L,
  ref_value = NULL,
  pred_at = c(0.25, 0.5, 0.75, 0.9, 0.95, 0.99),
  alpha = 0.05,
  lang = "pt",
  verbose = TRUE
)

Arguments

df: A climasus_df at stage "climate" and type "distributed_lag", produced by sus_climate_aggregate(temporal_strategy = "distributed_lag"). Must contain a date column (Date), the outcome column, and lag columns {climate_col}_lag0, ..., {climate_col}_lag{L}. Geometry columns (sf) are automatically dropped before modelling.
outcome_col: Character. Name of the daily health count column (e.g., "n_obitos", "n_internacoes", "n_dengue"). Default: "n_obitos".
climate_col: Character or NULL. Base name of the climate exposure (without _lagN suffix, e.g., "tair_dry_bulb_c"). NULL auto-detects by priority: air temperature > rainfall > humidity > other known variables.
lag_max: Integer or NULL. Maximum lag to include in the crossbasis. NULL auto-detects from the highest N present in {climate_col}_lagN columns.
covariates: Character vector or NULL. Names of additional columns in df to include as linear confounders in the GLM formula (e.g., c("rainfall_mm", "rh_mean_porc")). These are averaged across municipalities before fitting. Default: NULL.
argvar: Named list. Arguments for the exposure dimension of the crossbasis (passed to dlnm::crossbasis(argvar = ...)). Default: list(fun = "ns", df = 4). See Basis selection section.
arglag: Named list. Arguments for the lag dimension of the crossbasis (passed to dlnm::crossbasis(arglag = ...)). Default: list(fun = "ns", df = 3). See Basis selection section.
family: Character. GLM family: "quasipoisson" (default, recommended), "poisson", or "negbin" (redirected to "quasipoisson" with warning).
ns_df: Integer or NULL. Degrees of freedom for the natural spline time trend (splines::ns(date, df = ns_df)). NULL (default) triggers automatic computation based on dof_per_year. Set to 0 to suppress time control entirely.
dof_per_year: Integer. Degrees of freedom per year of data used for automatic ns_df computation. Default: 4L. Ignored when ns_df is explicitly set.
ref_value: Numeric or NULL. Exposure value where RR = 1 (the centring value for dlnm::crosspred()). NULL (default) uses the sample median of lag-0 exposure values.
pred_at: Numeric vector of quantile probabilities (0-1) at which to report the cumulative exposure-response curve in $exposure_response. Default: c(0.25, 0.50, 0.75, 0.90, 0.95, 0.99).
alpha: Numeric (0-1). Confidence level for intervals: CI = $(1 - \alpha)$. Default: 0.05 (95% CI).
lang: Character. Language for messages: "pt" (default), "en", "es".
verbose: Logical. Print progress steps and warnings. Default TRUE.

Value

A climasus_dlnm list with the following components:

$model: The fitted glm object. Includes the crossbasis in its model frame — compatible with dlnm::crosspred().
$crossbasis: The dlnm::crossbasis object encoding both the exposure and lag dimensions.
$pred: A dlnm::crosspred object computed on a 100-point grid from the 1st to 99th percentile of exposure. Contains allRRfit, allRRlow, allRRhigh (cumulative overall effect) and matRRfit, matRRlow, matRRhigh (lag-specific effects).
$exposure_response: Tibble with cumulative RR at the percentiles specified in pred_at: columns pct, exposure, rr, lo, hi.
$lag_response: Tibble with lag-specific RR at the 75th percentile of exposure (for dlnm_lag plots): columns lag, rr, lo, hi, rr_cum.
$models: One-row summary tibble for integration with plot_climate_health(fit = ...): variable, n, family, lag_max, ref_value, exposure_p75, rr, lo, hi, lag_peak, disp_ratio, aic_poisson.
$data_daily: Aggregated daily tibble used for model fitting (outcome summed; climate averaged across municipalities).
$diagnostics: List: disp_ratio, disp_category, autocorr_pval, has_autocorr, aic_poisson, deviance, null_deviance, df_residual.
$meta: List of modelling parameters plus the input sus_meta history for pipeline tracing.

The object also carries attr(result, "sus_meta") with an updated history entry.

Statistical framework

The DLNM models the log-linear association between an exposure history $X_{t-0}, X_{t-1}, \ldots, X_{t-L}$ and the daily count outcome $Y_t \sim \text{Quasi-Poisson}(\mu_t)$:

$$\log(\mu_t) = \alpha + \mathbf{w}_t^\top \boldsymbol{\theta} + \sum_{k=1}^{K} \gamma_k s_k(t)$$

where $\mathbf{w}_t$ is the cross-basis vector (outer product of the exposure and lag basis matrices), and $s_k(t)$ are natural splines controlling seasonal confounding. The cumulative overall effect across all lags at exposure level $x$ is:

$$RR_{\text{cum}}(x) = \exp\!\left\{\sum_{l=0}^{L} [f(x,l) - f(x_0, l)]\right\}$$

where $x_0$ is the reference value (default: sample median).

Family selection

"quasipoisson" (default, recommended): Robust to overdispersion ($\hat\phi > 1$). SE inflated by $\sqrt{\hat\phi}$. Recommended when disp_ratio > 1.5 (Armstrong, 2006).
"poisson": Valid only when $\hat\phi \approx 1$. Provides AIC/BIC but underestimates SE under overdispersion.
"negbin": Redirects to "quasipoisson" with a warning — negative binomial is incompatible with dlnm::crossbasis().

Basis selection (`argvar`, `arglag`)

`argvar`	When to use
`list(fun = "ns", df = 3)`	Mild non-linearity; conservative
`list(fun = "ns", df = 4)`	U/J shape (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 mortality
`list(fun = "ns", df = 4)`	Bimodal lag; cold effect
`list(fun = "poly", df = 2)`	Conservative; small samples

Temporal confounding (`ns_df`, `dof_per_year`)

A natural spline $\text{ns}(\text{date}, df)$ controls long-term trends and seasonality that confound the climate-health association (Bhaskaran et al., 2013). When ns_df = NULL the function auto-computes $df = \texttt{dof\_per\_year} \times n_{\text{years}}$. Recommended:

dof_per_year = 4: 1-2 year series; captures annual seasonality.
dof_per_year = 6: Subtropical Brazil; stronger seasonality.
dof_per_year = 8: Arboviral diseases; sub-seasonal variation.
ns_df = NULL with dof_per_year = 0: suppress time control (use only when exposure is already an anomaly from climatology).

References

Gasparrini, A., Armstrong, B., & Kenward, M.G. (2010). Distributed lag non-linear models. Statistics in Medicine, 29(21), 2224-2234. doi: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. doi: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. doi:10.1186/1471-2288-14-70

Armstrong, B. (2006). Models for the relationship between ambient temperature and daily mortality. Epidemiology, 17(6), 624-631.

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. doi:10.1093/ije/dyt092

Perkins, S.E. (2015). A review on the scientific understanding of heatwaves. Atmospheric Research, 164-165, 242-267.

Examples

if (FALSE) { # \dontrun{
# Step 1: produce distributed_lag exposure data
df_dl <- sus_climate_aggregate(
  health_data       = sf_sim_spatial,
  climate_data      = df_inmet_filled,
  climate_var       = "tair_dry_bulb_c",
  temporal_strategy = "distributed_lag",
  lag_days          = 21,
  lang              = "pt"
)

# Step 2: fit DLNM
fit <- sus_mod_dlnm(
  df          = df_dl,
  outcome_col = "n_obitos",
  lag_max     = 21,
  argvar      = list(fun = "ns", df = 4),
  arglag      = list(fun = "ns", df = 3),
  family      = "quasipoisson",
  dof_per_year = 4L,
  lang        = "pt"
)

print(fit)
fit$exposure_response   # RR cumulative at percentile grid
fit$lag_response        # RR by lag at p75
fit$diagnostics$disp_ratio

# Step 3: visualise (requires plot_climate_health)
# plots <- plot_climate_health(
#   data = df_dl, fit = fit,
#   plot_type = c("dlnm_surface", "dlnm_overall", "dlnm_lag")
# )
} # }