---
title: "bdlnm"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{bdlnm}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
inla_available <- bdlnm::check_inla()
```

```{r setup, message=FALSE}
library(bdlnm)
library(dlnm)
library(splines)
library(utils)
```

# london dataset

As an example for this vignette, we will use the built-in dataset `london` which contains daily temperature and mortality data from 2000 to 2012. 

```{r}
head(london)
```

Using Bayesian Distributed Non-Linear Models (B-DLNM) we will estimate the lagged effect of the temperature exposure to mortality in the population with more than 75 years old, taking 21 days of lags, in the city of London (2000-2012).

# Set DLNM parameters

We will start by building the crossbasis with the `dlnm` package.

We will use a natural spline with three knots for the exposure–response function, placed at the 10th, 75th, and 90th percentiles of daily mean temperature, and a natural spline for the lag–response function, with knots equally spaced on the log scale up to a maximum lag of 21 days:

```{r}
# Exposure-response and lag-response spline parameters
dlnm_var <- list(
    var_prc = c(10, 75, 90),
    var_fun = "ns",
    lag_fun = "ns",
    max_lag = 21,
    lagnk = 3
)


# Cross-basis parameters
argvar <- list(fun = dlnm_var$var_fun, 
               knots = quantile(london$tmean, 
                                dlnm_var$var_prc/100, na.rm = TRUE),
               Bound = range(london$tmean, na.rm = TRUE))

arglag <- list(fun = dlnm_var$lag_fun, 
               knots = logknots(dlnm_var$max_lag, nk = dlnm_var$lagnk))

# Create crossbasis
cb <- crossbasis(london$tmean, lag = dlnm_var$max_lag, argvar, arglag)
```

Let’s define the seasonality of the mortality time series using a natural spline with 8 degrees of freedom per year, which flexibly controls for long-term and seasonal trends in mortality:

```{r}
# Seasonality of mortality time series
seas <- ns(london$date, df = round(8 * length(london$date) / 365.25))
```

Finally, we also define the temperature values for which predictions will be generated. These correspond to an evenly spaced grid with a 0.1ºC step covering the full range of observed temperatures in the data:

```{r}
# Prediction values (equidistant points)
temp <- round(seq(min(london$tmean), max(london$tmean), by = 0.1), 1)
```

# Fit the model and posterior samples

Fit the Bayesian DLNM using the function `bdlnm()`. In this example, we specify a Poisson model including the cross-basis term, a seasonal component, and an indicator for the day of the week. We draw 1000 samples from the posterior distribution and set a seed for reproducibility.

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

The returned object is of class bdlnm and contains the fitted INLA model, the crossbasis objects included in the formula, and posterior samples of the estimated coefficients together with their summaries (mean, standard deviation, and quantiles).

```{r, eval = inla_available}
str(mod, max.level = 1)
```

# Minimum mortality temperature

We estimate the minimum mortality temperature (MMT), defined as the temperature at which the overall mortality risk is minimized. This optimal value will later be used to center the estimated relative risks.

```{r, eval = inla_available}
mmt <- optimal_exposure(mod, exp_at = temp)

str(mmt)
```

We can visualize the posterior distribution of the MMT estimates. It is useful to inspect this distribution to assess whether multiple candidate optimal exposure values exist and to verify that the median provides a reasonable centering value:

```{r fig.width = 7, eval = inla_available}
plot(mmt, xlab = "Temperature (ºC)", main = paste0("MMT (Median = ", round(mmt$summary[["0.5quant"]], 1), "ºC)"))
```

To make the predictions, we will center the risk at the median of these values:

```{r, eval = inla_available}
cen <- mmt$summary[["0.5quant"]] 
cen
```

# Predict exposure-lag-response effects

Predict the exposure–lag–response association between temperature and mortality from the fitted model at the supplied temperature grid, centering the predictions at the MMT value:

```{r, eval = inla_available}
cpred <- bcrosspred(mod, exp_at = temp, cen = cen)
```

A `bcrosspred` object is returned containing posterior samples of the estimated exposure–lag–response association for the supplied temperature values and lags, together with their summaries.

```{r, eval = inla_available}
str(cpred)
```

For instance, the estimated `crossbasis` coefficients are stored as:

```{r, eval = inla_available}
cpred$coefficients |> 
  head(c(5, 5))
```

The relative risks (RR) for each temperature–lag combination are also stored in an array for each posterior sample. For example, for the first sample:

```{r, eval = inla_available}
cpred$matRRfit[, , "sample1"] |> 
  head()
```

The overall cumulative effects for each temperature (summing across all lags) are stored in:

```{r, eval = inla_available}
cpred$allRRfit |> 
  head(c(5, 5))
```

Summaries of these effects across posterior samples are also available:

```{r, eval = inla_available}
cpred$matRRfit.summary |> 
  head(5)

cpred$allRRfit.summary |> 
  head(5)
```

# Plot

Several visualizations can be produced from these predictions using the `plot()` method.

For example, we can plot the overall cumulative effect for each temperature value (suming across all lags):

```{r fig.width = 7, fig.height = 4, eval = inla_available}
plot(cpred, "overall", xlab = "Temperature (ºC)", ylab = "Relative Risk", col = 4, main="Overall", log = "y")
```

Alternatively, we can display the posterior sampling variability by plotting the curves from all posterior samples instead of the credible interval:

```{r fig.width = 7, fig.height = 4, eval = inla_available}
plot(cpred, "overall", xlab = "Temperature (ºC)", ylab = "Relative Risk", col = 4, main="Overall", log = "y", ci = "sampling")
```

We can also visualize the full exposure–lag–response association using a 3-D surface:

```{r fig.width = 7, fig.height = 6, eval = inla_available}
plot(cpred, "3d", zlab = "Relative risk", col = 4, lphi = 60, cex.axis = 0.6, xlab = "Temperature (ºC)", main = "3D graph of temperature effect")
```

A contour plot provides a 2-D projection of the same relationship:

```{r fig.width = 7, fig.height = 6, eval = inla_available}
plot(cpred, "contour", xlab = "Temperature (ºC)", ylab = "Lag", main = "Contour plot")
```

We can also examine the lag–response association for a high temperature (e.g., the 99th percentile):

```{r fig.width = 7, fig.height = 6, eval = inla_available}
htemp <- round(quantile(london$tmean, 0.99), 1)
plot(cpred , "slices", ci = "bars", type = "p", pch = 19, exp_at = htemp, ylab="RR", main=paste0("Association for a high temperature (", htemp, "ºC)"))
```

Or the exposure–response association at lag 0:

```{r fig.width = 7, fig.height = 6, eval = inla_available}
plot(cpred , "slices", lag_at = 0, col=4, ylab="RR", main=paste0("Association at lag 0"))
```