---
title: "Advanced Topics"
author: "Ricardo G. Tomáz"
description: "Advanced model specification, uncertainty, diagnostics, and computational considerations"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Advanced Topics}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  message = FALSE,
  warning = FALSE
)
```

# Introduction

Most `EpiExposure` workflows can be completed using the default model
specification introduced in previous sections

However, advanced applications often require additional flexibility
regarding:

- model engines;
- uncertainty estimation;
- spline specification;
- identifiability diagnostics;
- computational performance.

This section discusses these advanced topics and provides
recommendations for robust DLNM analyses.

# Alternative model engines

`EpiExposure` separates exposure-lag specification from model
estimation.

Consequently, the same DLNM structure can be fitted using different
regression engines.

# Supported families by modeling engine

The modeling engines available in `fit_epidlnm()` and `find_bestfit()`
differ in the response families exposed through the `EpiExposure`
interface. The table below summarizes the current `EpiExposure` support.
It describes the package interface rather than every family accepted
natively by the underlying modeling packages.

| Engine    | Beta | Binomial       | Poisson | Gamma | Gaussian | Negative binomial |
|:----------|:----:|:---------------|:-------:|:-----:|:--------:|:-----------------:|
| `glm`     |  No  | Yes            |   Yes   |  Yes  |   Yes    |        No         |
| `glmmTMB` | Yes  | Yes            |   Yes   |  Yes  |   Yes    |     Yes, NB2      |
| `gam`     | Yes  | Yes            |   Yes   |  Yes  |   Yes    |     Yes, NB2      |
| `gamm`    |  No  | Yes            |   Yes   |  Yes  |   Yes    |        No         |
| `gls`     |  No  | No             |   No    |  No   |   Yes    |        No         |
| `spaMM`   | Yes  | Yes            |   Yes   |  Yes  |   Yes    |     Yes, NB2      |
| `brms`    | Yes  | Yes, Bernoulli |   Yes   |  Yes  |   Yes    |     Yes, NB2      |
| `INLA`    | Yes  | Yes            |   Yes   |  Yes  |   Yes    |     Yes, NB2      |
| `bdlnm`   | Yes  | Yes            |   Yes   |  Yes  |   Yes    |     Yes, NB2      |

**Note:** Negative-binomial support refers to the quadratic-variance NB2
parameterization. The table describes the families currently exposed
through the `EpiExposure` interface, not every family that may be
available natively in the underlying modeling packages. For INLA-backed
engines, availability also depends on the likelihoods provided by the
installed INLA version.

Because family specification follows the same `fit_epidlnm()` interface
across engines, the sections below do not repeat every supported
engine-family combination. Instead, each example highlights a capability
that differs meaningfully across engines, such as conventional random
effects, spatial autocorrelation, or Bayesian inference.

# Reference DLNM design

The examples below use the same cross-basis specification whenever
possible. This allows the modeling engines to be compared without
repeatedly rebuilding the exposure-lag design.

```{r}
library(EpiExposure)
```


```{r}
data("poisson_data")

head(poisson_data)
```

```{r}
cb <- define_exposures(
  data = poisson_data,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  max_lag = 85,
  df_var = 3,
  df_lag = 2
)

dat <- build_design(
  data = poisson_data,
  cb_templates = cb,
  groups = c(
    "epi_id",
    "block"
  )
)
```

```{r}
dat_poisson <- prepare_response(
  data = dat,
  response = "y",
  family = "poisson"
)
```

# Fixed-effect reference model

A generalized linear model provides the simplest fixed-effect
specification. It is useful as a reference when random or spatial
dependence is not required.

```{r}
fit_glm <- fit_epidlnm(
  data = dat_poisson,
  model_engine = "glm",
  family = "poisson",
  epiexposure_spec = attr(
    cb,
    "spec"
  )
)
```

```{r}
summary(fit_glm)
```

# Conventional random effects

When epidemics are grouped within experimental blocks, observations from
the same block may share unmeasured conditions. A random intercept can
account for this hierarchical dependence.

```{r}
fit_glmm <- fit_epidlnm(
  data = dat_poisson,
  model_engine = "glmmTMB",
  family = "poisson",
  random_effect = "block",
  epiexposure_spec = attr(
    cb,
    "spec"
  )
)
```

```{r}
summary(fit_glmm)
```

The term generated internally is equivalent to:

`(1 | block)`

This effect represents discrete group-level heterogeneity. It does not
model spatial correlation as a function of distance.

## GAM estimation

The `gam` engine can fit the DLNM design with `mgcv`. Beta and
negative-binomial models use `mgcv` extended families and are fitted
through the supported REML route by default.

```{r}
data("nb2_data")

head(nb2_data)
```

```{r}
dat_nb2 <- build_design(
  data = nb2_data,
  cb_templates = cb,
  groups = c(
    "epi_id",
    "block"
  )
)

dat_nb2 <- prepare_response(
  data = dat_nb2,
  response = "y",
  family = "negative_binomial"
)
```

```{r}
fit_gam <- fit_epidlnm(
  data = dat_nb2,
  model_engine = "gam",
  family = "negative_binomial",
  epiexposure_spec = attr(
    cb,
    "spec"
  )
)
```

```{r}
summary(fit_gam)
```

## GAMM estimation

The `gamm` engine combines the fixed DLNM design with a conventional
random intercept. In the current `EpiExposure` interface,
`random_effect` is required.

```{r}
fit_gamm <- fit_epidlnm(
  data = dat_poisson,
  model_engine = "gamm",
  family = "poisson",
  random_effect = "block",
  epiexposure_spec = attr(
    cb,
    "spec"
  )
)
```

The fitted object contains two native components:

```{r}
summary(fit_gamm)
```

```{r}
summary(
  fit_gamm,
  component = "lme"
)
```

The default summary describes the fixed/population GAM component. The
`component = "lme"` option displays the mixed-model component.

For non-Gaussian responses, `mgcv::gamm()` uses a PQL-based fitting
route. Binary responses should therefore be treated cautiously.

# Spatially autocorrelated effects with `spaMM`

Consider a regional study of a foliar rice disease conducted at
spatially distributed field locations. Daily mean temperature, leaf
wetness duration, and rainfall are recorded during the 86 days preceding
the final assessment. The outcome is the final number of lesions per
plant.

The cross-basis terms estimate non-linear and delayed environmental effects. A
Matérn random effect accounts for residual geographic similarity among
nearby fields after those measured environmental effects have been considered.

## Spatially dependent fields

```{r}
data("st_poisson")

head(st_poisson)
```

```{r, eval=FALSE}

cb_spatial <- define_exposures(
  data = st_poisson,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  max_lag = 85,
  df_var = 3,
  df_lag = 2
)

dat_spatial <- build_design(
  data = st_poisson,
  cb_templates = cb_spatial,
  groups = c(
    "epi_id",
    "x_coord",
    "y_coord"
  )
)

dat_spatial <- prepare_response(
  data = dat_spatial,
  response = "y",
  family = "poisson"
)
```

```{r, eval=FALSE}
fit_spatial <- fit_epidlnm(
  data = dat_spatial,
  model_engine = "spamm",
  family = "poisson",
  random_effect = NULL,
  spatial_effect = c(
    "x_coord",
    "y_coord"
  ),
  spatial_structure = "matern",
  spatial_group = NULL,
  epiexposure_spec = attr(
    cb_spatial,
    "spec"
  )
)

```

```{r, eval=FALSE}
summary(fit_spatial)
```

This specification generates:

`Matern(1 | x_coord + y_coord)`

The Matérn term is distinct from a conventional random intercept. It
represents autocorrelation that changes continuously with the distance
between locations.

## Independent spatial fields by year

If the disease process is expected to have a different spatial
realization in each growing season, `year` can define independent Matérn
fields:

```{r}
data("spatial_year")

head(spatial_year)
```

```{r, eval=FALSE}
dat_spatial_year <- build_design(
  data = spatial_year,
  cb_templates = cb,
  groups = c(
    "epi_id",
    "x_coord",
    "y_coord",
    "year",
    "location_id",
    "block_id"
  )
)

dat_spatial_year <- prepare_response(
  data = dat_spatial_year,
  response = "y",
  family = "poisson"
)
```

```{r, eval=FALSE}
fit_spatial_year <- fit_epidlnm(
  data = dat_spatial_year,
  model_engine = "spamm",
  family = "poisson",
  random_effect = "block_id",
  spatial_effect = c(
    "x_coord",
    "y_coord"
  ),
  spatial_structure = "matern",
  spatial_group = "year",
  epiexposure_spec = attr(
    cb_spatial,
    "spec"
  )
)

```

```{r, eval=FALSE}
summary(fit_spatial_year)
```

The resulting spatial term is:

`Matern(1 | x_coord + y_coord %in% year)`

In this model, `(1 | block_id)` represents conventional experimental
heterogeneity, whereas the Matérn term represents spatially
autocorrelated residual variation within each year.

## Independent spatial fields across epidemiological replicates

The grouping variable in a Matérn term can identify independent
repetitions of an entire epidemiological experiment. Consider a field
study conducted at 130 fixed locations. The same host, pathogen,
experimental layout, and observation period are used in four independent
repetitions.

Each repetition contains an 86-day exposure profile followed by one
final disease assessment at each location. Temperature, rainfall, and
leaf wetness can differ among repetitions, but the temporal positions
remain comparable: the most recent observation corresponds to lag zero,
and the oldest observation corresponds to lag 85.

Each combination of location and repetition represents a separate
epidemic:

`L001_Replicate_1`

`L001_Replicate_2`

`L001_Replicate_3`

`L001_Replicate_4`

The cross-basis terms use the exposure profile from each epidemic to
estimate the shared non-linear and delayed environmental associations.
However, the residual spatial pattern of disease is not required to
remain in the same locations across repetitions.

Independent residual spatial fields can be specified as:

`Matern(1 | x_coord + y_coord %in% epidemic_replicate)`

The model therefore allows each repetition to produce a different
realized map of residual disease intensity. For example, areas with high
residual disease in the first repetition do not need to remain high in
the second repetition.

The Matérn parameters remain shared among repetitions. Consequently, all
repetitions contribute jointly to estimation of the general magnitude,
smoothness, and spatial scale of residual dependence, while retaining
separate realized spatial fields.

```{r}
data("spatial_replicate")

head(spatial_replicate)
```

```{r, eval=FALSE}

cb_replicate <- define_exposures(
  data = spatial_replicate,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  max_lag = 85,
  df_var = 3,
  df_lag = 2
)

dat_spatial_replicate <- build_design(
  data = spatial_replicate,
  cb_templates = cb_replicate,
  groups = c(
    "epi_id",
    "x_coord",
    "y_coord",
    "epidemic_replicate",
    "location_id",
    "block_id"
    
  )
)

dat_spatial_replicate <- prepare_response(
  data = dat_spatial_replicate,
  response = "y",
  family = "poisson"
)
```

```{r, eval=FALSE}
fit_spatial_replicate <- fit_epidlnm(
  data = dat_spatial_replicate,
  model_engine = "spamm",
  family = "poisson",
  random_effect = "block_id",
  spatial_effect = c(
    "x_coord",
    "y_coord"
  ),
  spatial_structure = "matern",
  spatial_group = "epidemic_replicate",
  epiexposure_spec = attr(
    cb_replicate,
    "spec"
  )
)
```

```{r, eval=FALSE}
summary(fit_spatial_replicate)
```

The conventional random intercept and the spatial term represent
different sources of dependence. The `block_id` effect accounts for
discrete experimental heterogeneity, whereas the Matérn term accounts
for continuous residual similarity among nearby locations within each
epidemiological repetition.

An epidemiological replicate must represent an independent repetition of
the complete disease process. Repeated assessments made during the
progression of a single epidemic should not be treated as independent
replicates.

# Gaussian generalized least squares

The `gls` engine is restricted to Gaussian outcomes with the identity
link and does not support the `EpiExposure` random-intercept interface.

```{r}
data("gaussian_data")

head(gaussian_data)
```

```{r, eval=FALSE}
cb_gaussian <- define_exposures(
  data = gaussian_data,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  max_lag = 85,
  df_var = 3,
  df_lag = 2
)

dat_gaussian <- build_design(
  data = gaussian_data,
  cb_templates = cb_gaussian,
  groups = c(
    "epi_id"))

dat_gaussian <- prepare_response(
  data = dat_gaussian,
  response = "y",
  family = "gaussian"
)


```

```{r, eval=FALSE}
fit_gls <- fit_epidlnm(
  data = dat_gaussian,
  model_engine = "gls",
  family = "gaussian",
  epiexposure_spec = attr(
    cb_gaussian,
    "spec"
  ),
  method = "REML"
)

```

```{r, eval=FALSE}
summary(fit_gls)
```

 **Note:** The INLA and B-DLNM commands in the following sections are
 shown but not evaluated during vignette construction because these
 models may require substantial computational time. Users can execute
 the code interactively after installing the optional packages `INLA`,
 `bdlnm`, and `sn`.

# MCMC estimation with `brms`

The `brms` engine provides posterior sampling for the fixed DLNM
coefficients and any requested hierarchical components.

```{r, eval=FALSE}
fit_brms <- fit_epidlnm(
  data = dat_poisson,
  model_engine = "brms",
  family = "poisson",
  random_effect = "block",
  epiexposure_spec = attr(
    cb,
    "spec"
  ),
  chains = 4,
  iter = 2000,
  seed = 123
)


```

```{r, include = FALSE, eval=FALSE}
fit_brms = readRDS("data/fit_brms.rds")
```

```{r, eval=FALSE}
summary(fit_brms, eval=FALSE)
```


# Bayesian DLNM with `bdlnm`

The `bdlnm` engine provides a Bayesian implementation specifically
designed for distributed lag linear and non-linear models. Internally,
the native `bdlnm` package fits the model using Integrated Nested
Laplace Approximation through `INLA::inla()` and then draws samples from
the approximate posterior distribution using
`INLA::inla.posterior.sample()`. 【1-fe3db1】

Unlike the direct `inla` engine, which receives the epidemic-level
`cb_*` columns as conventional fixed predictors, the `bdlnm` engine
receives the original named `crossbasis` objects in the model formula.
The fitted object stores the underlying INLA model, the basis objects,
posterior coefficient draws, and summaries of those draws. 【1-fe3db1】

In `EpiExposure`, this engine is integrated into the same prediction and
uncertainty contract used by the other engines. Posterior fixed-effect
draws are retained for draw-by-draw propagation through
`summarise_effects()`, `predict_outcomes()`, and other downstream
functions.


```{r, eval=FALSE}
cb_bdlnm <- define_exposures(
  data = poisson_data,
  vars = c(
    "tmean",
    "rain",
    "wetness"
  ),
  max_lag = 85,
  df_var = 3,
  df_lag = 2
)


dat_bdlnm <- build_design(
  data = poisson_data,
  cb_templates = cb_bdlnm,
  groups = c(
    "epi_id",
    "block"
  )
)

dat_bdlnm <- prepare_response(
  data = dat_bdlnm,
  response = "y",
  family = "poisson"
)

bdlnm_bases <- attr(
  dat_bdlnm,
  "epiexposure_bdlnm_basis_objects",
  exact = TRUE
)
```

The model can then be fitted through the harmonized `EpiExposure`
interface:

```{r, eval=FALSE}
fit_bdlnm <- fit_epidlnm(
  data = dat_bdlnm,
  model_engine = "bdlnm",
  family = "poisson",
  random_effect = "block",
  random_effect_prior = list(
    prec = list(prior = "pc.prec",param = c(1,0.01))),
  epiexposure_spec = attr(dat_bdlnm, "epiexposure_spec"),
  sample.arg = list(n = 1000, seed = 123))
```

```{r, eval=FALSE}
summary(fit_bdlnm)
```

```{r, eval=FALSE}
fit_bdlnm$model
#fit_bdlnm$coefficients
fit_bdlnm$coefficients.summary


fit_bdlnm$model$summary.random
fit_bdlnm$model$summary.hyperpar

fit_bdlnm$model$summary.random$.epiexposure_re_index

head(fit_bdlnm$model$marginals.hyperpar)

```

The `random_effect = "block"` argument is translated internally to an
independent latent effect of the form

`f(block, model = "iid")`

for INLA-backed engines. When supplied, `random_effect_prior` is passed
directly to the `hyper` argument of the latent effect specification,
allowing custom hyperprior definitions for the random-effect precision.
For example,

`f(block, model = "iid", hyper = random_effect_prior)`

The `sample.arg` list controls posterior coefficient sampling performed
by the native `bdlnm` engine. A nonzero seed provides reproducible
posterior sampling. The native package uses 1,000 posterior samples by
default, although the number of samples can be changed according to the
required precision and available computational resources.

For deterministic `EpiExposure` predictions, posterior means of the
fixed/population coefficients define the central expected response. When
uncertainty is requested, the stored joint posterior coefficient draws
are propagated through the inverse link one draw at a time.

The resulting intervals describe uncertainty in the population-level
expected response. They do not include observation-level noise or
group-specific random effects, even when random effects are included in
the fitted model.

For both `model_engine = "inla"` and `model_engine = "bdlnm"`, random-
effect hyperparameters can be inspected through the underlying INLA fit
(e.g., `summary.hyperpar`), allowing users to evaluate posterior
distributions of random-effect precision parameters.

# Direct INLA versus B-DLNM

Both engines rely on INLA for Bayesian inference and support the same
likelihood families, random effects, and hyperprior specifications.
The main difference lies in how the DLNM basis is represented and what
is returned to the user.

- `model_engine = "inla"` fits the epidemic-level `cb_*` columns
  directly as fixed predictors through `INLA::inla()`. The resulting
  object contains the standard INLA outputs, including posterior
  summaries, marginal distributions, hyperparameter estimates, and
  random-effect summaries.

- `model_engine = "bdlnm"` passes named `crossbasis` objects to the
  native `bdlnm` interface. Internally, the model is fitted through
  INLA, but the original DLNM basis representation is preserved and
  posterior coefficient samples are generated automatically for
  simulation-based uncertainty propagation.
  
  
```{r, eval=FALSE}
fit_inla <- fit_epidlnm(
  data = dat_poisson,
  model_engine = "inla",
  family = "poisson",
  random_effect = "block",
  random_effect_prior = list(
      prec = list(
      prior = "pc.prec",
      param = c(1,0.01))),
  epiexposure_spec = attr(cb, "spec"))
```

### Outputs shared by INLA and B-DLNM

Fixed-effect posterior summaries:

```{r, eval=FALSE}
fit_inla$summary.fixed

fit_bdlnm$model$summary.fixed
```

Random-effect posterior summaries:

```{r, eval=FALSE}
fit_inla$summary.random

fit_bdlnm$model$summary.random
```

Posterior distributions of hyperparameters:

```{r, eval=FALSE}
fit_inla$marginals.hyperpar

fit_bdlnm$model$marginals.hyperpar
```

Marginal likelihood:

```{r, eval=FALSE}
fit_inla$mlik

fit_bdlnm$model$mlik
```

For both engines, random-effect hyperparameters can be inspected through
the INLA summaries. For example:

```{r, eval=FALSE}
fit_inla$summary.hyperpar

fit_bdlnm$model$summary.hyperpar
```

## Outputs specific to B-DLNM

```{r, eval=FALSE}
names(fit_bdlnm)
```

```{r, eval=FALSE}
names(fit_inla)
```

Original DLNM basis objects:

```{r, eval=FALSE}
fit_bdlnm$basis
```

Posterior coefficient draws:

```{r, eval=FALSE}

dim(fit_bdlnm$coefficients)

head(fit_bdlnm$coefficients)

```

Posterior coefficient summaries:

```{r, eval=FALSE}

fit_bdlnm$coefficients.summary
```

### Model objects

```{r, eval=FALSE}
class(fit_inla)

class(fit_bdlnm)

class(fit_bdlnm$model)
```

The direct INLA engine is useful when users want a conventional INLA
representation of the epidemic-level design matrix and direct access
to standard INLA outputs. The B-DLNM engine is useful when users want
to retain the native DLNM representation and work directly with
posterior coefficient samples, fitted basis objects, and Bayesian
DLNM workflows.